// @(#)root/pythia6:$Name:  $:$Id: TPythia6.cxx,v 1.19 2005/09/04 19:08:11 brun Exp $
// Author: Rene Brun   19/10/99
//
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// TPythia6                                                                   //
//                                                                            //
// TPythia is an interface class to F77 version of Pythia 6.2                 //
// To use this class you must install a version of pythia6.                   //
//   see the installation instructions at                                     //
//       http://root.cern.ch/root/Install.html                                //
//                                                                            // 
// CERNLIB event generators, written by T.Sjostrand.                          //
// For the details about these generators look at Pythia/Jetset manual:       //
//                                                                            //
//******************************************************************************
//**                                                                          **
//**                                                                          **
//**  PPP  Y   Y TTTTT H   H III   A        JJJJ EEEE TTTTT  SSS  EEEE TTTTT  **
//**  P  P  Y Y    T   H   H  I   A A          J E      T   S     E      T    **
//**  PPP    Y     T   HHHHH  I  AAAAA         J EEE    T    SSS  EEE    T    **
//**  P      Y     T   H   H  I  A   A      J  J E      T       S E      T    **
//**  P      Y     T   H   H III A   A       JJ  EEEE   T    SSS  EEEE   T    **
//**                                                                          **
//**                                                                          **
//**              *......*                  Welcome to the Lund Monte Carlo!  **
//**         *:::!!:::::::::::*                                               **
//**      *::::::!!::::::::::::::*            This is PYTHIA version 5.720    **
//**    *::::::::!!::::::::::::::::*        Last date of change: 29 Nov 1995  **
//**   *:::::::::!!:::::::::::::::::*                                         **
//**   *:::::::::!!:::::::::::::::::*         This is JETSET version 7.408    **
//**    *::::::::!!::::::::::::::::*!       Last date of change: 23 Aug 1995  **
//**      *::::::!!::::::::::::::* !!                                         **
//**      !! *:::!!:::::::::::*    !!                 Main author:            **
//**      !!     !* -><- *         !!              Torbjorn Sjostrand         **
//**      !!     !!                !!        Dept. of theoretical physics 2   **
//**      !!     !!                !!              University of Lund         **
//**      !!                       !!                Solvegatan 14A           **
//**      !!        ep             !!             S-223 62 Lund, Sweden       **
//**      !!                       !!          phone: +46 - 46 - 222 48 16    **
//**      !!                 pp    !!          E-mail: torbjorn@thep.lu.se    **
//**      !!   e+e-                !!                                         **
//**      !!                       !!         Copyright Torbjorn Sjostrand    **
//**      !!                                     and CERN, Geneva 1993        **
//**                                                                          **
//**                                                                          **
//** The latest program versions and documentation is found on WWW address    **
//** http://thep.lu.se/tf2/staff/torbjorn/Welcome.html                        **
//**                                                                          **
//** When you cite these programs, priority should always be given to the     **
//** latest published description. Currently this is                          **
//** T. Sjostrand, Computer Physics Commun. 82 (1994) 74.                     **
//** The most recent long description (unpublished) is                        **
//** T. Sjostrand, LU TP 95-20 and CERN-TH.7112/93 (revised August 1995).     **
//** Also remember that the programs, to a large extent, represent original   **
//** physics research. Other publications of special relevance to your        **
//** studies may therefore deserve separate mention.                          **
//**                                                                          **
//**                                                                          **
//******************************************************************************

#include "TPythia6.h"

#include "TClonesArray.h"
#include "TMCParticle6.h"
#include "TParticle.h"

TPythia6*  TPythia6::fgInstance = 0;

#ifndef WIN32
# define pyevnt pyevnt_
# define pyinit pyinit_
# define pychge pychge_
# define pycomp pycomp_
# define pyedit pyedit_
# define pyexec pyexec_
# define pyhepc pyhepc_
# define pylist pylist_
# define pymass pymass_
# define pyname pyname_
# define pyr    pyr_
# define pyrget pyrget_
# define pyrset pyrset_
# define pystat pystat_
# define pytest pytest_
# define pyupda pyupda_
# define tpythia6_open_fortran_file tpythia6_open_fortran_file_
# define tpythia6_close_fortran_file tpythia6_close_fortran_file_
# define type_of_call
#else
# define pyevnt PYEVNT
# define pyinit PYINIT
# define pychge PYCHGE
# define pycomp PYCOMP
# define pyedit PYEDIT
# define pyexec PYEXEC
# define pyhepc PYHEPC
# define pylist PYLIST
# define pymass PYMASS
# define pyname PYNAME
# define pyr    PYR
# define pyrget PYRGET
# define pyrset PYRSET
# define pystat PYSTAT
# define pytest PYTEST
# define pyupda PYUPDA
# define tpythia6_open_fortran_file TPYTHIA6_OPEN_FORTRAN_FILE
# define tpythia6_close_fortran_file TPYTHIA6_CLOSE_FORTRAN_FILE
# define type_of_call _stdcall
#endif


extern "C" void type_of_call pyevnt();
extern "C" void type_of_call pystat(int *key);
extern "C" void type_of_call pylist(int *key);
extern "C" void type_of_call pyedit(int *medit);
extern "C" void type_of_call pyexec();
extern "C" void type_of_call pyhepc(int *mconv);
extern "C" void type_of_call pylist(int *flag);
extern "C" int  type_of_call pychge(int *kf);
extern "C" int  type_of_call pycomp(int *kf);
extern "C" double  type_of_call pymass(int *kf);
extern "C" void type_of_call pyname(int *kf, char *name, Long_t l_name);
extern "C" int  type_of_call pyr(int *dummy);
extern "C" int  type_of_call pyrget(int *lun, int *move);
extern "C" int  type_of_call pyrset(int *lun, int *move);
extern "C" int  type_of_call pytest(int *flag);
extern "C" int  type_of_call pyupda(int *mupda, int *lun);

#ifndef WIN32
extern "C" void type_of_call pyinit(char *frame, char *beam, char *target,
                                     double *win, Long_t l_frame, Long_t l_beam,
                                     Long_t l_target);
#else
extern "C" void type_of_call pyinit(char *frame,  Long_t l_frame,
                                    char *beam,   Long_t l_beam,
                                    char *target, Long_t l_target,
                                    double *win
                                    );
#endif

extern "C" {
  void*  pythia6_common_address(const char*);
  void   type_of_call tpythia6_open_fortran_file(int* lun, char* name, int);
  void   type_of_call tpythia6_close_fortran_file(int* lun);
}

ClassImp(TPythia6)

//------------------------------------------------------------------------------
TPythia6::TPythia6Cleaner::TPythia6Cleaner() {
}

//------------------------------------------------------------------------------
TPythia6::TPythia6Cleaner::~TPythia6Cleaner() {

  if (TPythia6::fgInstance) {
    delete TPythia6::fgInstance;
    TPythia6::fgInstance = 0;
  }
}

//------------------------------------------------------------------------------
//  constructor is not supposed to be called from the outside - only
//  Initialize() method
//------------------------------------------------------------------------------
 TPythia6::TPythia6() : TGenerator("TPythia6","TPythia6") {
// TPythia6 constructor: creates a TClonesArray in which it will store all
// particles. Note that there may be only one functional TPythia6 object
// at a time, so it's not use to create more than one instance of it.

  delete fParticles; // was allocated as TObjArray in TGenerator

  fParticles = new TClonesArray("TMCParticle",50);

  // initialize common-blocks
  // the functions/subroutines referenced by TPythia6 can be found
  // at ftp://root.cern.ch/root/pythia6.tar.gz
  
  fPyjets = (Pyjets_t*) pythia6_common_address("PYJETS");
  fPydat1 = (Pydat1_t*) pythia6_common_address("PYDAT1");
  fPydat2 = (Pydat2_t*) pythia6_common_address("PYDAT2");
  fPydat3 = (Pydat3_t*) pythia6_common_address("PYDAT3");
  fPydat4 = (Pydat4_t*) pythia6_common_address("PYDAT4");
  fPydatr = (Pydatr_t*) pythia6_common_address("PYDATR");
  fPysubs = (Pysubs_t*) pythia6_common_address("PYSUBS");
  fPypars = (Pypars_t*) pythia6_common_address("PYPARS");
  fPyint1 = (Pyint1_t*) pythia6_common_address("PYINT1");
  fPyint2 = (Pyint2_t*) pythia6_common_address("PYINT2");
  fPyint3 = (Pyint3_t*) pythia6_common_address("PYINT3");
  fPyint4 = (Pyint4_t*) pythia6_common_address("PYINT4");
  fPyint5 = (Pyint5_t*) pythia6_common_address("PYINT5");
  fPyint6 = (Pyint6_t*) pythia6_common_address("PYINT6");
  fPyint7 = (Pyint7_t*) pythia6_common_address("PYINT7");
  fPyint8 = (Pyint8_t*) pythia6_common_address("PYINT8");
  fPyint9 = (Pyint9_t*) pythia6_common_address("PYINT9");
  fPyuppr = (Pyuppr_t*) pythia6_common_address("PYUPPR");
  fPymssm = (Pymssm_t*) pythia6_common_address("PYMSSM");
  fPyssmt = (Pyssmt_t*) pythia6_common_address("PYSSMT");
  fPyints = (Pyints_t*) pythia6_common_address("PYINTS");
  fPybins = (Pybins_t*) pythia6_common_address("PYBINS");
}

//------------------------------------------------------------------------------
  TPythia6::~TPythia6()
 {
// Destroys the object, deletes and disposes all TMCParticles currently on list.

    if (fParticles) {
      fParticles->Delete();
      delete fParticles;
      fParticles = 0;
   }
 }

//------------------------------------------------------------------------------
 TPythia6* TPythia6::Instance() {
  // model of automatic memory cleanup suggested by Jim Kowalkovski:
  // destructor for local static variable `cleaner' is  always called in the end
  // of the job thus deleting the only TPythia6 instance

  static TPythia6::TPythia6Cleaner cleaner;
  return fgInstance ? fgInstance : (fgInstance=new TPythia6()) ;
}





//______________________________________________________________________________
 void TPythia6::GenerateEvent() {

  //  generate event and copy the information from /HEPEVT/ to fPrimaries

  pyevnt();
  ImportParticles();
}

//______________________________________________________________________________
 void TPythia6::OpenFortranFile(int lun, char* name) {
  tpythia6_open_fortran_file(&lun, name, strlen(name));
}

//______________________________________________________________________________
 void TPythia6::CloseFortranFile(int lun) {
  tpythia6_close_fortran_file(&lun);
}


//______________________________________________________________________________
 TObjArray *TPythia6::ImportParticles(Option_t *)
{
// Fills TObjArray fParticles list with particles from common LUJETS
// Old contents of a list are cleared. This function should be called after
// any change in common LUJETS, however GetParticles() method  calls it
// automatically - user don't need to care about it. In case you make a call
// to LuExec() you must call this method yourself to transfer new data from
// common LUJETS to the fParticles list.

   fParticles->Clear();
   Int_t numpart   = fPyjets->N;
   TClonesArray &a = *((TClonesArray*)fParticles);
   for (Int_t i = 0; i<numpart; i++) {
      new(a[i]) TMCParticle(fPyjets->K[0][i] ,
                            fPyjets->K[1][i] ,
                            fPyjets->K[2][i] ,
                            fPyjets->K[3][i] ,
                            fPyjets->K[4][i] ,

                            fPyjets->P[0][i] ,
                            fPyjets->P[1][i] ,
                            fPyjets->P[2][i] ,
                            fPyjets->P[3][i] ,
                            fPyjets->P[4][i] ,

                            fPyjets->V[0][i] ,
                            fPyjets->V[1][i] ,
                            fPyjets->V[2][i] ,
                            fPyjets->V[3][i] ,
                            fPyjets->V[4][i]);
   }
   return fParticles;
}

//______________________________________________________________________________
 Int_t TPythia6::ImportParticles(TClonesArray *particles, Option_t *option)
{
//
//  Default primary creation method. It reads the /HEPEVT/ common block which
//  has been filled by the GenerateEvent method. If the event generator does
//  not use the HEPEVT common block, This routine has to be overloaded by
//  the subclasses.
//  The function loops on the generated particles and store them in
//  the TClonesArray pointed by the argument particles.
//  The default action is to store only the stable particles (ISTHEP = 1)
//  This can be demanded explicitly by setting the option = "Final"
//  If the option = "All", all the particles are stored.
//

  if (particles == 0) return 0;
  TClonesArray &clonesParticles = *particles;
  clonesParticles.Clear();
  Int_t numpart = fPyjets->N;
  Int_t nparts=0;
  if (!strcmp(option,"") || !strcmp(option,"Final")) {
    for (Int_t i = 0; i<numpart; i++) {

    if (fPyjets->K[0][i] == 1) {
//
//  Use the common block values for the TParticle constructor
//
        new(clonesParticles[nparts]) TParticle(
                            fPyjets->K[1][i] ,
                            fPyjets->K[0][i] ,
                            fPyjets->K[2][i] ,
                            -1,
                            fPyjets->K[3][i] ,
                            fPyjets->K[4][i] ,

                            fPyjets->P[0][i] ,
                            fPyjets->P[1][i] ,
                            fPyjets->P[2][i] ,
                            fPyjets->P[3][i] ,

                            fPyjets->V[0][i] ,
                            fPyjets->V[1][i] ,
                            fPyjets->V[2][i] ,
                            fPyjets->V[3][i]);

         //     if(gDebug) printf("%d %d %d! ",i,fPyjets->K[1][i],numpart);
         nparts++;
      }
    }
  }
  else if (!strcmp(option,"All")) {
    for (Int_t i = 0; i<numpart; i++) {
        new(clonesParticles[i]) TParticle(
                            fPyjets->K[1][i] ,
                            fPyjets->K[0][i] ,
                            fPyjets->K[2][i] ,
                            -1,
                            fPyjets->K[3][i] ,
                            fPyjets->K[4][i] ,

                            fPyjets->P[0][i] ,
                            fPyjets->P[1][i] ,
                            fPyjets->P[2][i] ,
                            fPyjets->P[3][i] ,

                            fPyjets->V[0][i] ,
                            fPyjets->V[1][i] ,
                            fPyjets->V[2][i] ,
                            fPyjets->V[3][i]);
    }
    nparts=numpart;
  }

  return nparts;
}

//______________________________________________________________________________
 void TPythia6::Initialize(const char *frame, const char *beam, const char *target, float win)
{
// Calls PyInit with the same parameters after performing some checking,
// sets correct title. This method should preferably be called instead of PyInit.
// PURPOSE: to initialize the generation procedure.
// ARGUMENTS: See documentation for details.
//    frame:  - specifies the frame of the experiment:
//                "CMS","FIXT","USER","FOUR","FIVE","NONE"
//    beam,
//    target: - beam and target particles (with additionaly cahrges, tildes or "bar":
//              e,nu_e,mu,nu_mu,tau,nu_tau,gamma,pi,n,p,Lambda,Sigma,Xi,Omega,
//              pomeron,reggeon
//    win:    - related to energy system:
//              for frame=="CMS" - total energy of system
//              for frame=="FIXT" - momentum of beam particle
//              for frame=="USER" - dummy - see documentation.
////////////////////////////////////////////////////////////////////////////////////

   char  cframe[4];
   strncpy(cframe,frame,4);
   char  cbeam[8];
   strncpy(cbeam,beam,8);
   char  ctarget[8];
   strncpy(ctarget,target,8);

   if ( (!strncmp(frame, "CMS"  ,3)) &&
        (!strncmp(frame, "FIXT" ,4)) &&
        (!strncmp(frame, "USER" ,4)) &&
        (!strncmp(frame, "FOUR" ,4)) &&
        (!strncmp(frame, "FIVE" ,4)) &&
        (!strncmp(frame, "NONE" ,4)) ) {
      printf("WARNING! In TPythia6:Initialize():\n");
      printf(" specified frame=%s is neither of CMS,FIXT,USER,FOUR,FIVE,NONE\n",frame);
      printf(" resetting to \"CMS\" .");
      sprintf(cframe,"CMS");
   }

   if ( (!strncmp(beam, "e"       ,1)) &&
        (!strncmp(beam, "nu_e"    ,4)) &&
        (!strncmp(beam, "mu"      ,2)) &&
        (!strncmp(beam, "nu_mu"   ,5)) &&
        (!strncmp(beam, "tau"     ,3)) &&
        (!strncmp(beam, "nu_tau"  ,6)) &&
        (!strncmp(beam, "gamma"   ,5)) &&
        (!strncmp(beam, "pi"      ,2)) &&
        (!strncmp(beam, "n"       ,1)) &&
        (!strncmp(beam, "p"       ,1)) &&
        (!strncmp(beam, "Lambda"  ,6)) &&
        (!strncmp(beam, "Sigma"   ,5)) &&
        (!strncmp(beam, "Xi"      ,2)) &&
        (!strncmp(beam, "Omega"   ,5)) &&
        (!strncmp(beam, "pomeron" ,7)) &&
        (!strncmp(beam, "reggeon" ,7)) ) {
      printf("WARNING! In TPythia6:Initialize():\n");
      printf(" specified beam=%s is unrecognized .\n",beam);
      printf(" resetting to \"p+\" .");
      sprintf(cbeam,"p+");
   }

   if ( (!strncmp(target, "e"       ,1)) &&
        (!strncmp(target, "nu_e"    ,4)) &&
        (!strncmp(target, "mu"      ,2)) &&
        (!strncmp(target, "nu_mu"   ,5)) &&
        (!strncmp(target, "tau"     ,3)) &&
        (!strncmp(target, "nu_tau"  ,6)) &&
        (!strncmp(target, "gamma"   ,5)) &&
        (!strncmp(target, "pi"      ,2)) &&
        (!strncmp(target, "n"       ,1)) &&
        (!strncmp(target, "p"       ,1)) &&
        (!strncmp(target, "Lambda"  ,6)) &&
        (!strncmp(target, "Sigma"   ,5)) &&
        (!strncmp(target, "Xi"      ,2)) &&
        (!strncmp(target, "Omega"   ,5)) &&
        (!strncmp(target, "pomeron" ,7)) &&
        (!strncmp(target, "reggeon" ,7)) ){
      printf("WARNING! In TPythia6:Initialize():\n");
      printf(" specified target=%s is unrecognized.\n",target);
      printf(" resetting to \"p+\" .");
      sprintf(ctarget,"p+");
   }



   Pyinit(cframe, cbeam ,ctarget, win);

   char atitle[32];
   sprintf(atitle," %s-%s at %g GeV",cbeam,ctarget,win);
   SetTitle(atitle);

}


 void TPythia6::Pyinit(char* frame, char* beam, char* target, double win) {
//------------------------------------------------------------------------------
// Calls Pyinit with the same parameters after performing some checking,
// sets correct title. This method should preferably be called instead of PyInit.
// PURPOSE: to initialize the generation procedure.
// ARGUMENTS: See documentation for details.
//    frame:  - specifies the frame of the experiment:
//                "CMS","FIXT","USER","FOUR","FIVE","NONE"
//    beam,
//    target: - beam and target particles (with additionaly charges,
//              tildes or "bar":
//              e,nu_e,mu,nu_mu,tau,nu_tau,gamma,pi,n,p,Lambda,Sigma,Xi,Omega,
//              pomeron,reggeon
//    win:    - related to energy system:
//              for frame=="CMS" - total energy of system
//              for frame=="FIXT" - momentum of beam particle
//              for frame=="USER" - dummy - see documentation.
//------------------------------------------------------------------------------

   Double_t lwin = win;
   Long_t  s1    = strlen(frame);
   Long_t  s2    = strlen(beam);
   Long_t  s3    = strlen(target);
#ifndef WIN32
   pyinit(frame,beam,target,&lwin,s1,s2,s3);
#else
   pyinit(frame, s1, beam , s2, target, s3, &lwin);
#endif
}


 int TPythia6::Pycomp(int kf) {
  return pycomp(&kf);
}

 void TPythia6::Pyedit(int medit) {
  pyedit(&medit);
  ImportParticles();
}

 void TPythia6::Pyevnt() {
  pyevnt();
}

 void TPythia6::Pyexec() {
  pyexec();
}

 void TPythia6::Pyhepc(int mconv) {
  pyhepc(&mconv);
}

 void TPythia6::Pylist(int flag) {
  pylist(&flag);
}

 void TPythia6::Pyname(int kf, char* name) {
  //Note that the array name must be dimensioned in the calling program
  //to at least name[16]
   
   pyname(&kf,name,15);
  // cut trailing blanks to get C string
 name[15] = 0;
  //for (int i=15; (i>=0) && (name[i] == ' '); i--) {
  //  name[i] = 0;
  // }
}

 double TPythia6::Pyr(int idummy) {
  return pyr(&idummy);
}

 void TPythia6::Pyrget(int lun, int move) {
  pyrget(&lun,&move);
}

 void TPythia6::Pyrset(int lun, int move) {
  pyrset(&lun,&move);
}

 void TPythia6::Pystat(int flag) {
  pystat(&flag);
}

 void TPythia6::Pytest(int flag) {
  pytest(&flag);
}

 void TPythia6::Pyupda(int mupda, int lun) {
  pyupda(&mupda,&lun);
}

 double TPythia6::Pymass(int kf) {
  return pymass(&kf);
}

 int TPythia6::Pychge(int kf) {
  return pychge(&kf);
}

//______________________________________________________________________________
 void TPythia6::SetupTest()
{
// Exemplary setup of Pythia parameters:
// Switches on processes 102,123,124 (Higgs generation) and switches off
// interactions, fragmentation, ISR, FSR...

   SetMSEL(0);            // full user controll;

   SetMSUB(102,1);        // g + g -> H0
   SetMSUB(123,1);        // f + f' -> f + f' + H0
   SetMSUB(124,1);        // f + f' -> f" + f"' + H0


   SetPMAS(6,1,175.0);   // mass of TOP
   SetPMAS(25,1,300);    // mass of Higgs


   SetCKIN(1,290.0);     // range of allowed mass
   SetCKIN(2,310.0);

   SetMSTP(61,  0);      // switch off ISR
   SetMSTP(71,  0);      // switch off FSR
   SetMSTP(81,  0);      // switch off multiple interactions
   SetMSTP(111, 0);      // switch off fragmentation/decay
}