[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
1a52e0ed	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"
07c71c2e	42	//*KEEP,TParticle,T=C++
1a52e0ed	43	#include "TParticle.h"
07c71c2e	44	#include "TClonesArray.h"
1a52e0ed	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_
db274076	61	# define dpmjet_openinp dpmjet_openinp_
06f923c8	62	# define dt_evtout dt_evtout_
1a52e0ed	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
06f923c8	74	# define dt_evtout DT_EVTOUT_
db274076	75	# define dpmjet_openinp DPMJET_OPENINP
1a52e0ed	76	# define type_of_call _stdcall
	77	#endif
	78
06f923c8	79
1a52e0ed	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();
06f923c8	88	extern "C" int type_of_call dt_evtout(const Int_t &);
1a52e0ed	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 &);
db274076	93	extern "C" void type_of_call dpmjet_openinp();
1a52e0ed	94
	95	#else
	96
	97	#endif
	98
	99	ClassImp(TDPMjet)
	100
f97d4a8e	101
1a52e0ed	102	//______________________________________________________________________________
f97d4a8e	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),
717765ce	118	fDecayAll(0),
5fed2e29	119	fProcess(kDpmMb),
5fed2e29	120	fFragmentation(kFALSE)
1a52e0ed	121	{
f97d4a8e	122	// Default Constructor
1a52e0ed	123	}
	124
	125	//______________________________________________________________________________
d38751f5	126	TDPMjet::TDPMjet(DpmProcess_t iproc, Int_t Ip=208, Int_t Ipz=82, Int_t It=208, Int_t Itz=82,
f97d4a8e	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),
2d22f1fe	135	fPpn(0.),
f97d4a8e	136	fCMEn(CMEn),
	137	fIdp(0),
	138	fBmin(0.),
	139	fBmax(0.),
	140	fFCentr(0),
	141	fPi0Decay(0),
717765ce	142	fDecayAll(0),
5fed2e29	143	fProcess(iproc),
7cbda79e	144	fFragmentation(kFALSE)
1a52e0ed	145	{
d38751f5	146	printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
1a52e0ed	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	}
1a52e0ed	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	//
1a52e0ed	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));
1a52e0ed	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];
d38751f5	224	/*
1a52e0ed	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]);
d38751f5	228	*/
06f923c8	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	*/
1a52e0ed	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]);
1a52e0ed	257	} // Particle loop
	258	}
	259	return nump;
fb941ba2	260
1a52e0ed	261	}
	262
	263
	264	//====================== access to dpmjet subroutines =========================
	265	//______________________________________________________________________________
	266	void TDPMjet::Initialize()
	267	{
d38751f5	268	//
	269	// Write standard DPMJET input cards
	270	//
7cbda79e	271	if(fFragmentation) printf("\tTDPMJet fragmentation/evaporation applied\n");
7cbda79e	272
d38751f5	273	FILE* out = fopen("dpmjet.inp","w");
d38751f5	274	// Projectile and Target definition
cd19c617	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
d38751f5	291	// Beam energy and crossing-angle
7cbda79e	292	fprintf(out, "CMENERGY %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fCMEn, 0., 0., 0., 0., 0.);
e2b2a1dd	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
5fed2e29	298	else if(fIp==1 && fIt>1){ // proton towards A side (directed z>0)
e2b2a1dd	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
7cbda79e	300	printf("\n TDPMjet::Initialize() -> p-A: projectile (p) energy = %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn);
e2b2a1dd	301	}
e2b2a1dd	302	else if(fIt==1 && fIp>1){ // proton towards C side (directed z<0)
7cbda79e	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
7cbda79e	304	printf("\n TDPMjet::Initialize() -> A-p: projectile (A) energy = %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn);
e2b2a1dd	305	}
e2b2a1dd	306	}
d38751f5	307	// Centrality
5fed2e29	308	if((fIp > 1 \|\| fIt > 1) && fFragmentation)
e2b2a1dd	309	fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-2., fBmin, fBmax, 0., 0., 0.);
5fed2e29	310	else if((fIp > 1 \|\| fIt > 1) && !fFragmentation)
5fed2e29	311	fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., fBmin, fBmax, 0., 0., 0.);
f97d4a8e	312	// Particle decays
f97d4a8e	313	if (fPi0Decay)
717765ce	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
d38751f5	317	//
	318	// PHOJET specific
	319	fprintf(out, "PHOINPUT\n");
	320	fprintf(out, "DEBUG 0 0 0 \n");
6a45631b	321
d38751f5	322	if (fProcess == kDpmMb) {
	323	fprintf(out, "PROCESS 1 0 1 1 1 1 1 1\n");
	324	} else if (fProcess == kDpmMbNonDiffr) {
cab41d67	325	fprintf(out, "PROCESS 1 0 1 0 0 0 0 1\n");
6a45631b	326	} else if (fProcess == kDpmDiffr) {
cab41d67	327	fprintf(out, "PROCESS 0 0 0 1 1 1 1 0\n");
6a45631b	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");
cab41d67	332	} else if (fProcess == kDpmCentralDiffr){
cab41d67	333	fprintf(out, "PROCESS 0 0 0 1 0 0 0 0\n");
d38751f5	334	}
d38751f5	335
06f923c8	336	Int_t iPDG[18] =
717765ce	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
06f923c8	345	4132
717765ce	346	};
	347
	348
	349	Int_t iON = (fDecayAll) ? 1:0;
06f923c8	350	for (Int_t i = 0; i < 8; i++) {
717765ce	351	fprintf(out, "LUND-DECAY%5d %5d\n", iPDG[i], iON);
	352	}
	353
	354	fprintf(out, "LUND-MSTJ %5d %5d\n", 22, 1);
	355
d38751f5	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);
db274076	362	dpmjet_openinp();
db274076	363
d38751f5	364	//
	365	// Call DPMJET initialisation
	366	Int_t iemu = 0; // No emulsion (default)
	367	Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
1a52e0ed	368
	369	}
	370
	371
	372	//______________________________________________________________________________
	373	void TDPMjet::GenerateEvent()
	374	{
	375	// Generates one event;
	376	fNEvent++;
	377	DTEVNO.nevent=fNEvent;
1a52e0ed	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);
1a52e0ed	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
d38751f5	390
d38751f5	391	printf("\n-------------------------------------------\n");
1a52e0ed	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);
d38751f5	397	printf("\n-------------------------------------------\n");
1a52e0ed	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
1a52e0ed	409	dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
1a52e0ed	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	}
3c349853	538
1a52e0ed	539	//______________________________________________________________________________
3c349853	540	Int_t TDPMjet::GetProjParticipants() const
1a52e0ed	541	{
	542	return DTGLCP.nwtaac;
	543	}
	544	//______________________________________________________________________________
3c349853	545	Int_t TDPMjet::GetTargParticipants() const
1a52e0ed	546	{
	547	return DTGLCP.nwtbac;
	548	}
9dcc0beb	549	//______________________________________________________________________________
	550	Int_t TDPMjet::GetProcessCode() const
	551	{
	552	return POPRCS.iproce;
	553	}
1a52e0ed	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
fb941ba2	578
492d4eba	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];}
fb941ba2	583