[u/mrichter/AliRoot.git] / TDPMjet / TDPMjet.cxx

//*KEEP,CopyRight,T=C.
/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/
//*KEND.

//*    +-------------------------------------------------------------+
//*    |                                                             |
//*    |                                                             |
//*    |                        DPMJET 3.0                           |
//*    |                                                             |
//*    |							     |
//*    |	 S. Roesler+), R. Engel#), J. Ranft*)		     |
//*    |							     |
//*    |	 +) CERN, TIS-RP				     |
//*    |	    CH-1211 Geneva 23, Switzerland		     |
//*    |	    Email: Stefan.Roesler@cern.ch		     |
//*    |							     |
//*    |	 #) University of Delaware, BRI 		     |
//*    |	    Newark, DE 19716, USA			     |
//*    |							     |
//*    |	 *) University of Siegen, Dept. of Physics	     |
//*    |	    D-57068 Siegen, Germany			     |
//*    |							     |
//*    |							     |
//*    |       http://home.cern.ch/sroesler/dpmjet3.html	     |
//*    |							     |
//*    |							     |
//*    |       Monte Carlo models used for event generation:	     |
//*    |	  PHOJET 1.12, JETSET 7.4 and LEPTO 6.5.1	     |
//*    |							     |
//*    +-------------------------------------------------------------+

//*KEEP,TDPMjet.
#include "TDPMjet.h"
//*KEEP,DPMCOMMON.
#include "DPMcommon.h"
//*KEEP,TParticle,T=C++
#include "TParticle.h"
#include "TClonesArray.h"
//*KEND.

//*KEEP,TROOT.
#include "TROOT.h"
//*KEND.

#ifndef WIN32
# define dt_dtuini dt_dtuini_
# define dt_getemu de_getemu_
# define dt_kkinc  dt_kkinc_
# define pho_phist pho_phist_
# define dt_dtuout dt_dtuout_
# define dt_rndm   dt_rndm_
# define dt_rndmst dt_rndmst_
# define dt_rndmin dt_rndmin_
# define dt_rndmou dt_rndmou_
# define dpmjet_openinp dpmjet_openinp_
# define dt_evtout dt_evtout_
# define type_of_call
#else
# define dt_dtuini DT_DTUINI
# define dt_getemu DT_GETEMU
# define dt_kkinc  DT_KKINC
# define pho_phist PHO_PHIST
# define dt_dtuout DT_DTUOUT
# define dt_rndm   DT_RNDM
# define dt_rndmst DT_RNDMST
# define dt_rndmin DT_RNDMIN
# define dt_rndmou DT_RNDMOU
# define dt_evtout DT_EVTOUT_
# define dpmjet_openinp DPMJET_OPENINP
# define type_of_call _stdcall
#endif


#ifndef WIN32
extern "C" void   type_of_call dt_dtuini(Int_t & , Double_t &, Int_t & , Int_t &, 
			Int_t &, Int_t &, Int_t &, Int_t &);
extern "C" double type_of_call dt_getemu(Int_t &, Int_t &, Int_t &, Int_t &);
extern "C" void   type_of_call dt_kkinc(Int_t &, Int_t &, Int_t &, Int_t &,
				    Int_t &, Double_t &, Int_t &, Int_t &);
extern "C" void   type_of_call pho_phist(Int_t &, Double_t &);
extern "C" void   type_of_call dt_dtuout();
extern "C" int    type_of_call dt_evtout(const Int_t &);
extern "C" void   type_of_call dt_rndm(Int_t &);
extern "C" void   type_of_call dt_rndmst(Int_t &, Int_t &, Int_t &, Int_t &);
extern "C" void   type_of_call dt_rndmin(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
extern "C" void   type_of_call dt_rndmou(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
extern "C" void   type_of_call dpmjet_openinp();

#else

#endif

ClassImp(TDPMjet)


//______________________________________________________________________________
    TDPMjet::TDPMjet() : 
	TGenerator("dpmjet","dpmjet"),
	fNEvent(0),
	fIp(0),
	fIpz(0),
	fIt(0),
	fItz(0),
	fEpn(0.),
	fPpn(0.),
	fCMEn(0.),
	fIdp(0),
	fBmin(0.),
	fBmax(0.),
	fFCentr(0),
	fPi0Decay(0),
	fDecayAll(0),
	fProcess(kDpmMb),
	fFragmentation(kFALSE)
{
// Default Constructor
}

//______________________________________________________________________________
TDPMjet::TDPMjet(DpmProcess_t  iproc, Int_t Ip=208, Int_t Ipz=82, Int_t It=208, Int_t Itz=82, 
		 Double_t Epn=2700., Double_t CMEn=5400.) 
    : TGenerator("dpmjet","dpmjet"),
      fNEvent(0),
      fIp(Ip),
      fIpz(Ipz),
      fIt(It),
      fItz(Itz),
      fEpn(Epn),
      fPpn(0.),
      fCMEn(CMEn),
      fIdp(0),
      fBmin(0.),
      fBmax(0.),
      fFCentr(0),
      fPi0Decay(0),
      fDecayAll(0),
      fProcess(iproc),
      fFragmentation(kFALSE)
{  
    printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
}


//______________________________________________________________________________
Int_t TDPMjet::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 
//  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 &Particles = *particles;
  Particles.Clear();
  Int_t numpart = 0; 	 // Total number of produced particles
  Int_t numStabpart = 0; // Total number of produced stable particles
  Double_t entot = 0; 	 // Total energy in final state (from stable particles)
  
  numpart = DTEVT1.nhkk;
  for(Int_t i=0; i<numpart; i++){
     if(DTEVT1.isthkk[i]==1 || DTEVT1.isthkk[i]==-1 || DTEVT1.isthkk[i]==1001){
        numStabpart++;
	entot += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
     } 
  }
  Int_t nump = 0;
  
  if(!strcmp(option,"") || !strcmp(option,"Final")){
      for (Int_t i=0; i < numpart; i++) {
	  
	  if (DTEVT1.isthkk[i] == 1) {
	   //
	   //  Use the common block values for the TParticle constructor
	   //
	      new(Particles[nump]) TParticle(
		  DTEVT1.idhkk[i],
		  DTEVT1.isthkk[i],
		  -1,
		  -1,
		  -1,
		  -1,
		  DTEVT1.phkk[i][0],
		  DTEVT1.phkk[i][1],
		  DTEVT1.phkk[i][2],
		  DTEVT1.phkk[i][3],
		  
		  DTEVT1.vhkk[i][0],
		  DTEVT1.vhkk[i][1],
		  DTEVT1.vhkk[i][2],
		  DTEVT1.vhkk[i][3]);
	      nump++;
	  }
      }
  }
  else if(!strcmp(option,"All")){
      nump = numpart; 
      for (Int_t i=0; i <= numpart; i++){
        
	  // DTEVT1.JMOHKK[i][0] pointer to the entry of the 1st mother of entry i
	  Int_t iParent = DTEVT1.jmohkk[i][0] - 1;
	     
	  if(iParent >= 0){
	      TParticle *mother = (TParticle*) (Particles.UncheckedAt(iParent));
	      mother->SetLastDaughter(i);
	      if(mother->GetFirstDaughter() == -1) mother->SetFirstDaughter(i);
	  } 
	  // --- PDGcode for residual nuclei (idhkk=80000) 
	  // --- 	10000*Z + 10*A
	  // --- DPMJET -> idres = mass #, idxres = charge
	  if(DTEVT1.idhkk[i] == 80000) 
	      DTEVT1.idhkk[i] = 10000*DTEVT2.idxres[i]+10*DTEVT2.idres[i];
/*
	  if(DTEVT2.idxres[i] != 0) 
	      printf("\n	pc#%d -> A = %d, Z = %d -> PDGcode = %d\n",
		     i,DTEVT2.idres[i],DTEVT2.idxres[i],DTEVT1.idhkk[i]);
*/	  
/*
	  printf("%5d %5d %5d %5d %13.3f %13.3f %13.3f %13.3f \n", i,
		 DTEVT1.idhkk[i],
		 DTEVT1.isthkk[i],
		 iParent,  
		 DTEVT1.phkk[i][0],
		 DTEVT1.phkk[i][1],
		 DTEVT1.phkk[i][2],
		 DTEVT1.phkk[i][3]
		 );
*/
	  new(Particles[i]) TParticle(
	      DTEVT1.idhkk[i],
	      DTEVT1.isthkk[i],
	      iParent,
	      -1,
	      -1,
	      -1,
	      
	      DTEVT1.phkk[i][0],
	      DTEVT1.phkk[i][1],
	      DTEVT1.phkk[i][2],
	      DTEVT1.phkk[i][3],
	      
	      DTEVT1.vhkk[i][0],
	      DTEVT1.vhkk[i][1],
	      DTEVT1.vhkk[i][2],
	      DTEVT1.vhkk[i][3]);
      } // Particle loop  
  }
  return nump;
  
}


//====================== access to dpmjet subroutines =========================
//______________________________________________________________________________
void TDPMjet::Initialize()
{
//
//  Write standard DPMJET input cards 
//
    if(fFragmentation) printf("\tTDPMJet fragmentation/evaporation applied\n");

    FILE* out = fopen("dpmjet.inp","w");
//  Projectile and Target definition 
    if (fIp == 1 && fIpz ==1) {
	fprintf(out, "PROJPAR                                                               PROTON\n");
    } else if (fIp == 1 && fIpz == -1) {
	fprintf(out, "PROJPAR                                                               APROTON\n");
    } else {
	fprintf(out, "PROJPAR   %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", (Float_t) fIp, (Float_t) fIpz,  0., 0., 0., 0.);
    }
    
    if (fIt == 1 && fItz ==1) {
	fprintf(out, "TARPAR                                                                PROTON\n");
    } else if (fIt == 1 && fItz == -1) {
	fprintf(out, "TARPAR                                                                APROTON\n");
    } else {
	fprintf(out, "TARPAR    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", (Float_t) fIt, (Float_t) fItz,  0., 0., 0., 0.);
    }

//  Beam energy and crossing-angle
    fprintf(out, "CMENERGY      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fCMEn, 0., 0., 0., 0., 0.);      
    if(fIt == 1 && fIp ==1){ 
      fprintf(out, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.); //p-p
    }
    else if(fIp > 1 || fIt > 1){ 
      if(fIp>1 && fIt>1) fprintf(out, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.);//A-A
      else if(fIp==1 && fIt>1){ // proton towards A side (directed z>0)
        fprintf(out, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", fEpn,fEpn*fItz/fIt, 0., 0., 0., 0.);//pA
	printf("\n  TDPMjet::Initialize() -> p-A: projectile (p) energy =  %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn);
      }
      else if(fIt==1 && fIp>1){ // proton towards C side (directed z<0)
        fprintf(out, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", fEpn, fEpn*fIp/fIpz, 0., 0., 0., 0.);//A-p
	printf("\n  TDPMjet::Initialize() -> A-p: projectile (A) energy =  %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn);
      }
    }
//  Centrality
    if((fIp > 1 || fIt > 1) && fFragmentation)
	fprintf(out, "CENTRAL   %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-2., fBmin, fBmax, 0., 0., 0.);
    else if((fIp > 1 || fIt > 1) && !fFragmentation)
	fprintf(out, "CENTRAL   %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., fBmin, fBmax, 0., 0., 0.);
//  Particle decays
    if (fPi0Decay) 
	fprintf(out, "PARDECAY  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 2., 0., 0., 0., 0., 0.);    

    
//
//  PHOJET specific
    fprintf(out, "PHOINPUT\n");
    fprintf(out, "DEBUG      0 0 0 \n");

    if (fProcess == kDpmMb) {
	fprintf(out, "PROCESS           1 0 1 1 1 1 1 1\n");
    } else if (fProcess == kDpmMbNonDiffr) {
	fprintf(out, "PROCESS           1 0 1 0 0 0 0 1\n");
    } else if (fProcess == kDpmDiffr) {
	fprintf(out, "PROCESS           0 0 0 1 1 1 1 0\n");
    }else if (fProcess == kDpmSingleDiffr) {
        fprintf(out, "PROCESS           0 0 0 0 1 1 0 0\n");
    }else if (fProcess == kDpmDoubleDiffr) {
        fprintf(out, "PROCESS           0 0 0 0 0 0 1 0\n");
    } else if (fProcess == kDpmCentralDiffr){
        fprintf(out, "PROCESS           0 0 0 1 0 0 0 0\n");
    }
    
    Int_t iPDG[18] = 
	{
//          K0s   pi0  lam   sig+  sig-  tet0
	310,  111, 3122, 3222, 3112, 3322,
//          tet- om-    D+      D0     Ds+
	3312, 3334,  411,  421,  431,
//          etac lamc+ sigc++ sigc+ sigc0 Ksic+
        441, 4122, 4222, 4212, 4112, 4232,
//         Ksic0 sig0 
	4132
	};
    
    
    Int_t iON = (fDecayAll) ? 1:0;
    for (Int_t i = 0; i < 8; i++) {
	fprintf(out, "LUND-DECAY%5d %5d\n",  iPDG[i], iON);    
    }
	
    fprintf(out, "LUND-MSTJ %5d %5d\n",   22, 1);    

    fprintf(out, "ENDINPUT\n");
//
//  START card
    fprintf(out, "START            1.0       0.0\n");
    fprintf(out, "STOP\n");
    fclose(out);
    dpmjet_openinp();

//
//  Call DPMJET initialisation
    Int_t iemu = 0; // No emulsion (default)
    Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);

}


//______________________________________________________________________________
void TDPMjet::GenerateEvent()
{
   // Generates one event;
   fNEvent++;
   DTEVNO.nevent=fNEvent;
   Int_t kkmat=-1;
   Float_t Elab = fEpn;
   Int_t irej=0;
   Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
}
//______________________________________________________________________________
void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar, 
   			int ntmass, int ntchar, int idp, int iemu)
{
  // Call dmpjet routine DT_DTUINI passing the parameters 
  // in a way accepted by Fortran routines				   
     

   printf("\n-------------------------------------------\n");
   printf("\n		Dt_Dtuini called with:\n\n");
   printf(" Projectile	-> A = %d, Z = %d \n",npmass, npchar);
   printf(" Target    	-> A = %d, Z = %d \n",ntmass, ntchar);
   printf(" Proj. LAB E	-> E = %f GeV \n",epn);
   printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
   printf("\n-------------------------------------------\n");

   dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
    
}

//______________________________________________________________________________
void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar, 
   		       int idp, double elab, int kkmat, int irej)
{
  // Call dmpjet routine DT_KKINC passing the parameters 
  // in a way accepted by Fortran routines				   
  dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
}

//______________________________________________________________________________
void TDPMjet::Pho_Phist(int imode, double weight)
{
  // Call dmpjet routine PHO_PHIST passing the parameters 
  // in a way accepted by Fortran routines
  
  pho_phist(imode,weight);				   

}

//______________________________________________________________________________
void TDPMjet::Dt_Dtuout()
{
  // Call dmpjet routine DT_DTUOT passing the parameters 
  // in a way accepted by Fortran routines				   
  
  dt_dtuout();

}

//______________________________________________________________________________
Int_t TDPMjet::GetEvNum() const
{
	return DTEVT1.nevhkk;
}
//______________________________________________________________________________
Int_t TDPMjet::GetEntriesNum() const
{
	return DTEVT1.nhkk;
}
//______________________________________________________________________________
Int_t TDPMjet::GetNumStablePc() const
{
	Int_t NumStablePc = 0;
	for(Int_t i=0; i<DTEVT1.nhkk; i++){
	   if(DTEVT1.isthkk[i] == 1) NumStablePc++;
	}
	return NumStablePc;
}

//______________________________________________________________________________
Float_t TDPMjet::GetTotEnergy() const
{
  	Float_t TotEnergy = 0.;
  	for(Int_t i=0; i<DTEVT1.nhkk; i++){
     	  if(DTEVT1.isthkk[i] == 1)
	    TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
  	}
	return TotEnergy;
}

//______________________________________________________________________________
Int_t TDPMjet::GetStatusCode(Int_t evnum) const 
{
	return DTEVT1.isthkk[evnum];	
}
//______________________________________________________________________________
Int_t TDPMjet::GetPDGCode(Int_t evnum) const   
{
	return DTEVT1.idhkk[evnum];
}
//______________________________________________________________________________
Double_t TDPMjet::Getpx(Int_t evnum) const       
{
	return DTEVT1.phkk[evnum][0];
}
//______________________________________________________________________________
Double_t TDPMjet::Getpy(Int_t evnum) const      
{
	return DTEVT1.phkk[evnum][1];
}
//______________________________________________________________________________
Double_t TDPMjet::Getpz(Int_t evnum) const       
{
	return DTEVT1.phkk[evnum][2];
}
//______________________________________________________________________________
Double_t TDPMjet::GetEnergy(Int_t evnum) const	      
{
	return DTEVT1.phkk[evnum][3];
}
//______________________________________________________________________________
Double_t TDPMjet::GetMass(Int_t evnum) const 	      
{
	return DTEVT1.phkk[evnum][4];
}
//______________________________________________________________________________
Int_t    TDPMjet::GetFragmentA(Int_t evnum) const	
{
	return DTEVT2.idres[evnum];
}
//______________________________________________________________________________
Int_t    TDPMjet::GetFragmentZ(Int_t evnum) const	
{
	return DTEVT2.idxres[evnum];
}
//______________________________________________________________________________
Double_t TDPMjet::GetXSFrac() const 	      
{
	return DTIMPA.xsfrac;
}
//______________________________________________________________________________
Double_t TDPMjet::GetBImpac() const 	      
{
	return DTGLCP.bimpac;
}
//______________________________________________________________________________
Double_t TDPMjet::GetProjRadius() const 	      
{
	return DTGLCP.rproj;
}
//______________________________________________________________________________
Double_t TDPMjet::GetTargRadius() const 	      
{
	return DTGLCP.rtarg;
}
//______________________________________________________________________________
Int_t TDPMjet::GetProjWounded() const
{
	return DTGLCP.nwasam;
}
//______________________________________________________________________________
Int_t TDPMjet::GetTargWounded() const
{
	return DTGLCP.nwbsam;
}

//______________________________________________________________________________
Int_t TDPMjet::GetProjParticipants() const
{
	return DTGLCP.nwtaac;
}
//______________________________________________________________________________
Int_t TDPMjet::GetTargParticipants() const
{
	return DTGLCP.nwtbac;
}
//______________________________________________________________________________
Int_t TDPMjet::GetProcessCode() const
{
		return POPRCS.iproce;
}
//______________________________________________________________________________
void TDPMjet::Dt_Rndm(int idummy)
{
	dt_rndm(idummy);
}

//______________________________________________________________________________
void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
{
	dt_rndmst(na1, na2, na3, nb1);
}

//______________________________________________________________________________
void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
{
	dt_rndmin(u, c, cd, cm, i, j);
}

//______________________________________________________________________________
void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
{
	dt_rndmou(u, c, cd, cm, i, j);
}


Int_t TDPMjet::NHEP()                      const {return POEVT1.nhep;}
Int_t TDPMjet::ISTHEP(Int_t i)             const {return POEVT1.isthep[i];}
Int_t TDPMjet::IDHEP(Int_t i)              const {return POEVT1.idhep[i];}
Double_t TDPMjet::PHEP(Int_t i, Int_t j)   const {return POEVT1.phep[i][j];}
Commit	Line	Data
	1	//*KEEP,CopyRight,T=C.
	2	/*************************************************************************
	3	* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
	4	* All rights reserved. *
	5	* *
	6	* For the licensing terms see $ROOTSYS/LICENSE. *
	7	* For the list of contributors see $ROOTSYS/README/CREDITS. *
	8	*************************************************************************/
	9	//*KEND.
	10
	11	//* +-------------------------------------------------------------+
	12	//* \| \|
	13	//* \| \|
	14	//* \| DPMJET 3.0 \|
	15	//* \| \|
	16	//* \| \|
	17	//* \| S. Roesler+), R. Engel#), J. Ranft*) \|
	18	//* \| \|
	19	//* \| +) CERN, TIS-RP \|
	20	//* \| CH-1211 Geneva 23, Switzerland \|
	21	//* \| Email: Stefan.Roesler@cern.ch \|
	22	//* \| \|
	23	//* \| #) University of Delaware, BRI \|
	24	//* \| Newark, DE 19716, USA \|
	25	//* \| \|
	26	//* \| *) University of Siegen, Dept. of Physics \|
	27	//* \| D-57068 Siegen, Germany \|
	28	//* \| \|
	29	//* \| \|
	30	//* \| http://home.cern.ch/sroesler/dpmjet3.html \|
	31	//* \| \|
	32	//* \| \|
	33	//* \| Monte Carlo models used for event generation: \|
	34	//* \| PHOJET 1.12, JETSET 7.4 and LEPTO 6.5.1 \|
	35	//* \| \|
	36	//* +-------------------------------------------------------------+
	37
	38	//*KEEP,TDPMjet.
	39	#include "TDPMjet.h"
	40	//*KEEP,DPMCOMMON.
	41	#include "DPMcommon.h"
	42	//*KEEP,TParticle,T=C++
	43	#include "TParticle.h"
	44	#include "TClonesArray.h"
	45	//*KEND.
	46
	47	//*KEEP,TROOT.
	48	#include "TROOT.h"
	49	//*KEND.
	50
	51	#ifndef WIN32
	52	# define dt_dtuini dt_dtuini_
	53	# define dt_getemu de_getemu_
	54	# define dt_kkinc dt_kkinc_
	55	# define pho_phist pho_phist_
	56	# define dt_dtuout dt_dtuout_
	57	# define dt_rndm dt_rndm_
	58	# define dt_rndmst dt_rndmst_
	59	# define dt_rndmin dt_rndmin_
	60	# define dt_rndmou dt_rndmou_
	61	# define dpmjet_openinp dpmjet_openinp_
	62	# define dt_evtout dt_evtout_
	63	# define type_of_call
	64	#else
	65	# define dt_dtuini DT_DTUINI
	66	# define dt_getemu DT_GETEMU
	67	# define dt_kkinc DT_KKINC
	68	# define pho_phist PHO_PHIST
	69	# define dt_dtuout DT_DTUOUT
	70	# define dt_rndm DT_RNDM
	71	# define dt_rndmst DT_RNDMST
	72	# define dt_rndmin DT_RNDMIN
	73	# define dt_rndmou DT_RNDMOU
	74	# define dt_evtout DT_EVTOUT_
	75	# define dpmjet_openinp DPMJET_OPENINP
	76	# define type_of_call _stdcall
	77	#endif
	78
	79
	80	#ifndef WIN32
	81	extern "C" void type_of_call dt_dtuini(Int_t & , Double_t &, Int_t & , Int_t &,
	82	Int_t &, Int_t &, Int_t &, Int_t &);
	83	extern "C" double type_of_call dt_getemu(Int_t &, Int_t &, Int_t &, Int_t &);
	84	extern "C" void type_of_call dt_kkinc(Int_t &, Int_t &, Int_t &, Int_t &,
	85	Int_t &, Double_t &, Int_t &, Int_t &);
	86	extern "C" void type_of_call pho_phist(Int_t &, Double_t &);
	87	extern "C" void type_of_call dt_dtuout();
	88	extern "C" int type_of_call dt_evtout(const Int_t &);
	89	extern "C" void type_of_call dt_rndm(Int_t &);
	90	extern "C" void type_of_call dt_rndmst(Int_t &, Int_t &, Int_t &, Int_t &);
	91	extern "C" void type_of_call dt_rndmin(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
	92	extern "C" void type_of_call dt_rndmou(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
	93	extern "C" void type_of_call dpmjet_openinp();
	94
	95	#else
	96
	97	#endif
	98
	99	ClassImp(TDPMjet)
	100
	101
	102	//______________________________________________________________________________
	103	TDPMjet::TDPMjet() :
	104	TGenerator("dpmjet","dpmjet"),
	105	fNEvent(0),
	106	fIp(0),
	107	fIpz(0),
	108	fIt(0),
	109	fItz(0),
	110	fEpn(0.),
	111	fPpn(0.),
	112	fCMEn(0.),
	113	fIdp(0),
	114	fBmin(0.),
	115	fBmax(0.),
	116	fFCentr(0),
	117	fPi0Decay(0),
	118	fDecayAll(0),
	119	fProcess(kDpmMb),
	120	fFragmentation(kFALSE)
	121	{
	122	// Default Constructor
	123	}
	124
	125	//______________________________________________________________________________
	126	TDPMjet::TDPMjet(DpmProcess_t iproc, Int_t Ip=208, Int_t Ipz=82, Int_t It=208, Int_t Itz=82,
	127	Double_t Epn=2700., Double_t CMEn=5400.)
	128	: TGenerator("dpmjet","dpmjet"),
	129	fNEvent(0),
	130	fIp(Ip),
	131	fIpz(Ipz),
	132	fIt(It),
	133	fItz(Itz),
	134	fEpn(Epn),
	135	fPpn(0.),
	136	fCMEn(CMEn),
	137	fIdp(0),
	138	fBmin(0.),
	139	fBmax(0.),
	140	fFCentr(0),
	141	fPi0Decay(0),
	142	fDecayAll(0),
	143	fProcess(iproc),
	144	fFragmentation(kFALSE)
	145	{
	146	printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
	147	}
	148
	149
	150	//______________________________________________________________________________
	151	Int_t TDPMjet::ImportParticles(TClonesArray particles, Option_t option)
	152	{
	153	//
	154	// Default primary creation method. It reads the /HEPEVT/ common block which
	155	// has been filled by the GenerateEvent method. If the event generator does
	156	// not use the HEPEVT common block, This routine has to be overloaded by
	157	// the subclasses.
	158	// The function loops on the generated particles and store them in
	159	// the TClonesArray pointed by the argument particles.
	160	// The default action is to store only the stable particles
	161	// This can be demanded explicitly by setting the option = "Final"
	162	// If the option = "All", all the particles are stored.
	163	//
	164	if(particles==0) return 0;
	165	TClonesArray &Particles = *particles;
	166	Particles.Clear();
	167	Int_t numpart = 0; // Total number of produced particles
	168	Int_t numStabpart = 0; // Total number of produced stable particles
	169	Double_t entot = 0; // Total energy in final state (from stable particles)
	170
	171	numpart = DTEVT1.nhkk;
	172	for(Int_t i=0; i<numpart; i++){
	173	if(DTEVT1.isthkk[i]==1 \|\| DTEVT1.isthkk[i]==-1 \|\| DTEVT1.isthkk[i]==1001){
	174	numStabpart++;
	175	entot += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
	176	}
	177	}
	178	Int_t nump = 0;
	179
	180	if(!strcmp(option,"") \|\| !strcmp(option,"Final")){
	181	for (Int_t i=0; i < numpart; i++) {
	182
	183	if (DTEVT1.isthkk[i] == 1) {
	184	//
	185	// Use the common block values for the TParticle constructor
	186	//
	187	new(Particles[nump]) TParticle(
	188	DTEVT1.idhkk[i],
	189	DTEVT1.isthkk[i],
	190	-1,
	191	-1,
	192	-1,
	193	-1,
	194	DTEVT1.phkk[i][0],
	195	DTEVT1.phkk[i][1],
	196	DTEVT1.phkk[i][2],
	197	DTEVT1.phkk[i][3],
	198
	199	DTEVT1.vhkk[i][0],
	200	DTEVT1.vhkk[i][1],
	201	DTEVT1.vhkk[i][2],
	202	DTEVT1.vhkk[i][3]);
	203	nump++;
	204	}
	205	}
	206	}
	207	else if(!strcmp(option,"All")){
	208	nump = numpart;
	209	for (Int_t i=0; i <= numpart; i++){
	210
	211	// DTEVT1.JMOHKK[i][0] pointer to the entry of the 1st mother of entry i
	212	Int_t iParent = DTEVT1.jmohkk[i][0] - 1;
	213
	214	if(iParent >= 0){
	215	TParticle mother = (TParticle) (Particles.UncheckedAt(iParent));
	216	mother->SetLastDaughter(i);
	217	if(mother->GetFirstDaughter() == -1) mother->SetFirstDaughter(i);
	218	}
	219	// --- PDGcode for residual nuclei (idhkk=80000)
	220	// --- 10000Z + 10A
	221	// --- DPMJET -> idres = mass #, idxres = charge
	222	if(DTEVT1.idhkk[i] == 80000)
	223	DTEVT1.idhkk[i] = 10000DTEVT2.idxres[i]+10DTEVT2.idres[i];
	224	/*
	225	if(DTEVT2.idxres[i] != 0)
	226	printf("\n pc#%d -> A = %d, Z = %d -> PDGcode = %d\n",
	227	i,DTEVT2.idres[i],DTEVT2.idxres[i],DTEVT1.idhkk[i]);
	228	*/
	229	/*
	230	printf("%5d %5d %5d %5d %13.3f %13.3f %13.3f %13.3f \n", i,
	231	DTEVT1.idhkk[i],
	232	DTEVT1.isthkk[i],
	233	iParent,
	234	DTEVT1.phkk[i][0],
	235	DTEVT1.phkk[i][1],
	236	DTEVT1.phkk[i][2],
	237	DTEVT1.phkk[i][3]
	238	);
	239	*/
	240	new(Particles[i]) TParticle(
	241	DTEVT1.idhkk[i],
	242	DTEVT1.isthkk[i],
	243	iParent,
	244	-1,
	245	-1,
	246	-1,
	247
	248	DTEVT1.phkk[i][0],
	249	DTEVT1.phkk[i][1],
	250	DTEVT1.phkk[i][2],
	251	DTEVT1.phkk[i][3],
	252
	253	DTEVT1.vhkk[i][0],
	254	DTEVT1.vhkk[i][1],
	255	DTEVT1.vhkk[i][2],
	256	DTEVT1.vhkk[i][3]);
	257	} // Particle loop
	258	}
	259	return nump;
	260
	261	}
	262
	263
	264	//====================== access to dpmjet subroutines =========================
	265	//______________________________________________________________________________
	266	void TDPMjet::Initialize()
	267	{
	268	//
	269	// Write standard DPMJET input cards
	270	//
	271	if(fFragmentation) printf("\tTDPMJet fragmentation/evaporation applied\n");
	272
	273	FILE* out = fopen("dpmjet.inp","w");
	274	// Projectile and Target definition
	275	if (fIp == 1 && fIpz ==1) {
	276	fprintf(out, "PROJPAR PROTON\n");
	277	} else if (fIp == 1 && fIpz == -1) {
	278	fprintf(out, "PROJPAR APROTON\n");
	279	} else {
	280	fprintf(out, "PROJPAR %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", (Float_t) fIp, (Float_t) fIpz, 0., 0., 0., 0.);
	281	}
	282
	283	if (fIt == 1 && fItz ==1) {
	284	fprintf(out, "TARPAR PROTON\n");
	285	} else if (fIt == 1 && fItz == -1) {
	286	fprintf(out, "TARPAR APROTON\n");
	287	} else {
	288	fprintf(out, "TARPAR %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", (Float_t) fIt, (Float_t) fItz, 0., 0., 0., 0.);
	289	}
	290
	291	// Beam energy and crossing-angle
	292	fprintf(out, "CMENERGY %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fCMEn, 0., 0., 0., 0., 0.);
	293	if(fIt == 1 && fIp ==1){
	294	fprintf(out, "BEAM %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.); //p-p
	295	}
	296	else if(fIp > 1 \|\| fIt > 1){
	297	if(fIp>1 && fIt>1) fprintf(out, "BEAM %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.);//A-A
	298	else if(fIp==1 && fIt>1){ // proton towards A side (directed z>0)
	299	fprintf(out, "BEAM %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", fEpn,fEpn*fItz/fIt, 0., 0., 0., 0.);//pA
	300	printf("\n TDPMjet::Initialize() -> p-A: projectile (p) energy = %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn);
	301	}
	302	else if(fIt==1 && fIp>1){ // proton towards C side (directed z<0)
	303	fprintf(out, "BEAM %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", fEpn, fEpn*fIp/fIpz, 0., 0., 0., 0.);//A-p
	304	printf("\n TDPMjet::Initialize() -> A-p: projectile (A) energy = %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn);
	305	}
	306	}
	307	// Centrality
	308	if((fIp > 1 \|\| fIt > 1) && fFragmentation)
	309	fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-2., fBmin, fBmax, 0., 0., 0.);
	310	else if((fIp > 1 \|\| fIt > 1) && !fFragmentation)
	311	fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., fBmin, fBmax, 0., 0., 0.);
	312	// Particle decays
	313	if (fPi0Decay)
	314	fprintf(out, "PARDECAY %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 2., 0., 0., 0., 0., 0.);
	315
	316
	317	//
	318	// PHOJET specific
	319	fprintf(out, "PHOINPUT\n");
	320	fprintf(out, "DEBUG 0 0 0 \n");
	321
	322	if (fProcess == kDpmMb) {
	323	fprintf(out, "PROCESS 1 0 1 1 1 1 1 1\n");
	324	} else if (fProcess == kDpmMbNonDiffr) {
	325	fprintf(out, "PROCESS 1 0 1 0 0 0 0 1\n");
	326	} else if (fProcess == kDpmDiffr) {
	327	fprintf(out, "PROCESS 0 0 0 1 1 1 1 0\n");
	328	}else if (fProcess == kDpmSingleDiffr) {
	329	fprintf(out, "PROCESS 0 0 0 0 1 1 0 0\n");
	330	}else if (fProcess == kDpmDoubleDiffr) {
	331	fprintf(out, "PROCESS 0 0 0 0 0 0 1 0\n");
	332	} else if (fProcess == kDpmCentralDiffr){
	333	fprintf(out, "PROCESS 0 0 0 1 0 0 0 0\n");
	334	}
	335
	336	Int_t iPDG[18] =
	337	{
	338	// K0s pi0 lam sig+ sig- tet0
	339	310, 111, 3122, 3222, 3112, 3322,
	340	// tet- om- D+ D0 Ds+
	341	3312, 3334, 411, 421, 431,
	342	// etac lamc+ sigc++ sigc+ sigc0 Ksic+
	343	441, 4122, 4222, 4212, 4112, 4232,
	344	// Ksic0 sig0
	345	4132
	346	};
	347
	348
	349	Int_t iON = (fDecayAll) ? 1:0;
	350	for (Int_t i = 0; i < 8; i++) {
	351	fprintf(out, "LUND-DECAY%5d %5d\n", iPDG[i], iON);
	352	}
	353
	354	fprintf(out, "LUND-MSTJ %5d %5d\n", 22, 1);
	355
	356	fprintf(out, "ENDINPUT\n");
	357	//
	358	// START card
	359	fprintf(out, "START 1.0 0.0\n");
	360	fprintf(out, "STOP\n");
	361	fclose(out);
	362	dpmjet_openinp();
	363
	364	//
	365	// Call DPMJET initialisation
	366	Int_t iemu = 0; // No emulsion (default)
	367	Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
	368
	369	}
	370
	371
	372	//______________________________________________________________________________
	373	void TDPMjet::GenerateEvent()
	374	{
	375	// Generates one event;
	376	fNEvent++;
	377	DTEVNO.nevent=fNEvent;
	378	Int_t kkmat=-1;
	379	Float_t Elab = fEpn;
	380	Int_t irej=0;
	381	Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
	382	}
	383	//______________________________________________________________________________
	384	void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar,
	385	int ntmass, int ntchar, int idp, int iemu)
	386	{
	387	// Call dmpjet routine DT_DTUINI passing the parameters
	388	// in a way accepted by Fortran routines
	389
	390
	391	printf("\n-------------------------------------------\n");
	392	printf("\n Dt_Dtuini called with:\n\n");
	393	printf(" Projectile -> A = %d, Z = %d \n",npmass, npchar);
	394	printf(" Target -> A = %d, Z = %d \n",ntmass, ntchar);
	395	printf(" Proj. LAB E -> E = %f GeV \n",epn);
	396	printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
	397	printf("\n-------------------------------------------\n");
	398
	399	dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
	400
	401	}
	402
	403	//______________________________________________________________________________
	404	void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar,
	405	int idp, double elab, int kkmat, int irej)
	406	{
	407	// Call dmpjet routine DT_KKINC passing the parameters
	408	// in a way accepted by Fortran routines
	409	dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
	410	}
	411
	412	//______________________________________________________________________________
	413	void TDPMjet::Pho_Phist(int imode, double weight)
	414	{
	415	// Call dmpjet routine PHO_PHIST passing the parameters
	416	// in a way accepted by Fortran routines
	417
	418	pho_phist(imode,weight);
	419
	420	}
	421
	422	//______________________________________________________________________________
	423	void TDPMjet::Dt_Dtuout()
	424	{
	425	// Call dmpjet routine DT_DTUOT passing the parameters
	426	// in a way accepted by Fortran routines
	427
	428	dt_dtuout();
	429
	430	}
	431
	432	//______________________________________________________________________________
	433	Int_t TDPMjet::GetEvNum() const
	434	{
	435	return DTEVT1.nevhkk;
	436	}
	437	//______________________________________________________________________________
	438	Int_t TDPMjet::GetEntriesNum() const
	439	{
	440	return DTEVT1.nhkk;
	441	}
	442	//______________________________________________________________________________
	443	Int_t TDPMjet::GetNumStablePc() const
	444	{
	445	Int_t NumStablePc = 0;
	446	for(Int_t i=0; i<DTEVT1.nhkk; i++){
	447	if(DTEVT1.isthkk[i] == 1) NumStablePc++;
	448	}
	449	return NumStablePc;
	450	}
	451
	452	//______________________________________________________________________________
	453	Float_t TDPMjet::GetTotEnergy() const
	454	{
	455	Float_t TotEnergy = 0.;
	456	for(Int_t i=0; i<DTEVT1.nhkk; i++){
	457	if(DTEVT1.isthkk[i] == 1)
	458	TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
	459	}
	460	return TotEnergy;
	461	}
	462
	463	//______________________________________________________________________________
	464	Int_t TDPMjet::GetStatusCode(Int_t evnum) const
	465	{
	466	return DTEVT1.isthkk[evnum];
	467	}
	468	//______________________________________________________________________________
	469	Int_t TDPMjet::GetPDGCode(Int_t evnum) const
	470	{
	471	return DTEVT1.idhkk[evnum];
	472	}
	473	//______________________________________________________________________________
	474	Double_t TDPMjet::Getpx(Int_t evnum) const
	475	{
	476	return DTEVT1.phkk[evnum][0];
	477	}
	478	//______________________________________________________________________________
	479	Double_t TDPMjet::Getpy(Int_t evnum) const
	480	{
	481	return DTEVT1.phkk[evnum][1];
	482	}
	483	//______________________________________________________________________________
	484	Double_t TDPMjet::Getpz(Int_t evnum) const
	485	{
	486	return DTEVT1.phkk[evnum][2];
	487	}
	488	//______________________________________________________________________________
	489	Double_t TDPMjet::GetEnergy(Int_t evnum) const
	490	{
	491	return DTEVT1.phkk[evnum][3];
	492	}
	493	//______________________________________________________________________________
	494	Double_t TDPMjet::GetMass(Int_t evnum) const
	495	{
	496	return DTEVT1.phkk[evnum][4];
	497	}
	498	//______________________________________________________________________________
	499	Int_t TDPMjet::GetFragmentA(Int_t evnum) const
	500	{
	501	return DTEVT2.idres[evnum];
	502	}
	503	//______________________________________________________________________________
	504	Int_t TDPMjet::GetFragmentZ(Int_t evnum) const
	505	{
	506	return DTEVT2.idxres[evnum];
	507	}
	508	//______________________________________________________________________________
	509	Double_t TDPMjet::GetXSFrac() const
	510	{
	511	return DTIMPA.xsfrac;
	512	}
	513	//______________________________________________________________________________
	514	Double_t TDPMjet::GetBImpac() const
	515	{
	516	return DTGLCP.bimpac;
	517	}
	518	//______________________________________________________________________________
	519	Double_t TDPMjet::GetProjRadius() const
	520	{
	521	return DTGLCP.rproj;
	522	}
	523	//______________________________________________________________________________
	524	Double_t TDPMjet::GetTargRadius() const
	525	{
	526	return DTGLCP.rtarg;
	527	}
	528	//______________________________________________________________________________
	529	Int_t TDPMjet::GetProjWounded() const
	530	{
	531	return DTGLCP.nwasam;
	532	}
	533	//______________________________________________________________________________
	534	Int_t TDPMjet::GetTargWounded() const
	535	{
	536	return DTGLCP.nwbsam;
	537	}
	538
	539	//______________________________________________________________________________
	540	Int_t TDPMjet::GetProjParticipants() const
	541	{
	542	return DTGLCP.nwtaac;
	543	}
	544	//______________________________________________________________________________
	545	Int_t TDPMjet::GetTargParticipants() const
	546	{
	547	return DTGLCP.nwtbac;
	548	}
	549	//______________________________________________________________________________
	550	Int_t TDPMjet::GetProcessCode() const
	551	{
	552	return POPRCS.iproce;
	553	}
	554	//______________________________________________________________________________
	555	void TDPMjet::Dt_Rndm(int idummy)
	556	{
	557	dt_rndm(idummy);
	558	}
	559
	560	//______________________________________________________________________________
	561	void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
	562	{
	563	dt_rndmst(na1, na2, na3, nb1);
	564	}
	565
	566	//______________________________________________________________________________
	567	void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
	568	{
	569	dt_rndmin(u, c, cd, cm, i, j);
	570	}
	571
	572	//______________________________________________________________________________
	573	void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
	574	{
	575	dt_rndmou(u, c, cd, cm, i, j);
	576	}
	577
	578
	579	Int_t TDPMjet::NHEP() const {return POEVT1.nhep;}
	580	Int_t TDPMjet::ISTHEP(Int_t i) const {return POEVT1.isthep[i];}
	581	Int_t TDPMjet::IDHEP(Int_t i) const {return POEVT1.idhep[i];}
	582	Double_t TDPMjet::PHEP(Int_t i, Int_t j) const {return POEVT1.phep[i][j];}
	583