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

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


// Generator using DPMJET as an external generator
// The main DPMJET options are accessable for the user through this interface.
// Uses the TDPMjet implementation of TGenerator.

#include <TDPMjet.h>
#include <TRandom.h>
#include <TArrayI.h>
#include <TParticle.h>
#include <TGraph.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
#include <TParticleClassPDG.h>
#include <TPDGCode.h>
#include <TLorentzVector.h>
#include "AliRunLoader.h"
#include "AliGenDPMjet.h"
#include "AliGenDPMjetEventHeader.h"
#include "AliRun.h"
#include "AliDpmJetRndm.h"
#include "AliHeader.h"
#include "AliStack.h"
#include "AliMC.h"

ClassImp(AliGenDPMjet)

//______________________________________________________________________________
AliGenDPMjet::AliGenDPMjet()
    :AliGenMC(), 
     fBeamEn(2750.),
     fEnergyCMS(5500.),
     fMinImpactParam(0.),
     fMaxImpactParam(5.),
     fICentr(0),
     fSelectAll(0),
     fFlavor(0),
     fTrials(0),
     fSpectators(1),
     fSpecn(0),
     fSpecp(0),
     fDPMjet(0),
     fParticles(0),
     fNoGammas(0),
     fLHC(0),
     fPi0Decay(0),
     fGenImpPar(0.),
     fProcess(kDpmMb)
{
// Constructor
    AliDpmJetRndm::SetDpmJetRandom(GetRandom());
}


//______________________________________________________________________________
AliGenDPMjet::AliGenDPMjet(Int_t npart)
    :AliGenMC(npart),
     fBeamEn(2750.),
     fEnergyCMS(5500.),
     fMinImpactParam(0.),
     fMaxImpactParam(5.),
     fICentr(0),
     fSelectAll(0),
     fFlavor(0),
     fTrials(0),
     fSpectators(1),
     fSpecn(0),
     fSpecp(0),
     fDPMjet(0),
     fParticles(new TClonesArray("TParticle",10000)),
     fNoGammas(0),
     fLHC(0),
     fPi0Decay(0),
     fGenImpPar(0.),
     fProcess(kDpmMb)
{
// Default PbPb collisions at 5. 5 TeV
//
    fName = "DPMJET";
    fTitle= "Particle Generator using DPMJET";
    SetTarget();
    SetProjectile();
    fVertex.Set(3);
    AliDpmJetRndm::SetDpmJetRandom(GetRandom());
}

AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
    :AliGenMC(),
     fBeamEn(2750.),
     fEnergyCMS(5500.),
     fMinImpactParam(0.),
     fMaxImpactParam(5.),
     fICentr(0),
     fSelectAll(0),
     fFlavor(0),
     fTrials(0),
     fSpectators(1),
     fSpecn(0),
     fSpecp(0),
     fDPMjet(0),
     fParticles(0),
     fNoGammas(0),
     fLHC(0),
     fPi0Decay(0),
     fGenImpPar(0.),
     fProcess(kDpmMb)
{
    // Dummy copy constructor
}

//______________________________________________________________________________
AliGenDPMjet::~AliGenDPMjet()
{
// Destructor
    delete fParticles;
}
//______________________________________________________________________________
void AliGenDPMjet::Init()
{
// Initialization
    
    SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget, 
		      fBeamEn,fEnergyCMS));

    fDPMjet=(TDPMjet*) fMCEvGen;
    //
    // **** Flag to force central production
    // fICentr=1. central production forced 
    // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam	  
    // fICentr<-99 -> fraction of x-sec. = XSFRAC		  
    // fICentr=-1. -> evaporation/fzc suppressed		  
    // fICentr<-1. -> evaporation/fzc suppressed		  
    if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
    
    fDPMjet->SetfFCentr(fICentr);  
    fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
    fDPMjet->SetPi0Decay(fPi0Decay);
//
//  Initialize DPMjet  
//    
    fDPMjet->Initialize();
}


//______________________________________________________________________________
void AliGenDPMjet::Generate()
{
// Generate one event

  Float_t polar[3]    =   {0,0,0};
  Float_t origin[3]   =   {0,0,0};
  Float_t origin0[3]  =   {0,0,0};
  Float_t p[3];
  Float_t tof;

//  converts from mm/c to s
  const Float_t kconv = 0.001/2.999792458e8;
  Int_t nt  = 0;
  Int_t jev = 0;
  Int_t kf, ks, imo;
  kf = 0;
  fTrials = 0;
  //  Set collision vertex position 
  if (fVertexSmear == kPerEvent) Vertex();
  
  while(1)
  {
//    Generate one event
// --------------------------------------------------------------------------
      fSpecn = 0;  
      fSpecp = 0;
// --------------------------------------------------------------------------
      fDPMjet->GenerateEvent();
      fTrials++;

      fDPMjet->ImportParticles(fParticles,"All");      
      if (fLHC) Boost();

      // Temporaneo
      fGenImpPar = fDPMjet->GetBImpac();
      
      Int_t np = fParticles->GetEntriesFast();
      printf("\n **************************************************%d\n",np);
      Int_t nc=0;
      if (np==0) continue;
      Int_t i;
      Int_t* newPos     = new Int_t[np];
      Int_t* pSelected  = new Int_t[np];

      for (i = 0; i<np; i++) {
	  newPos[i]    = i;
	  pSelected[i] = 0;
      }
      
//      Get event vertex
      
      fVertex[0] = origin0[0];
      fVertex[1] = origin0[1];	
      fVertex[2] = origin0[2];
      
//      First select parent particles

      for (i = 0; i<np; i++) {
	  TParticle *iparticle = (TParticle *) fParticles->At(i);

// Is this a parent particle ?

	  if (Stable(iparticle)) continue;

	  Bool_t  selected             =  kTRUE;
	  Bool_t  hasSelectedDaughters =  kFALSE;
	  
	  kf = iparticle->GetPdgCode();
	  if (kf == 92) continue;
	  ks = iparticle->GetStatusCode();
// No initial state partons
          if (ks==21) continue;
	    
	  if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
			       SelectFlavor(kf);
	  hasSelectedDaughters = DaughtersSelection(iparticle);

// Put particle on the stack if it is either selected or 
// it is the mother of at least one seleted particle

	  if (selected || hasSelectedDaughters) {
	      nc++;
	      pSelected[i] = 1;
	  } // selected
      } // particle loop parents

// Now select the final state particles


      for (i=0; i<np; i++) {
	  TParticle *iparticle = (TParticle *) fParticles->At(i);

// Is this a final state particle ?

	  if (!Stable(iparticle)) continue;
      
	  Bool_t  selected =  kTRUE;
	  kf = iparticle->GetPdgCode();
	  ks = iparticle->GetStatusCode();

// --------------------------------------------------------------------------
// Count spectator neutrons and protons (ks == 13, 14)
	  if(ks == 13 || ks == 14){
	      if(kf == kNeutron) fSpecn += 1;
	      if(kf == kProton)  fSpecp += 1;
	  }
// --------------------------------------------------------------------------

	  if (!fSelectAll) {
	      selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
	      if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
	  }

// Put particle on the stack if selected

	  if (selected) {
	      nc++;
	      pSelected[i] = 1;
	  } // selected
      } // particle loop final state

// Write particles to stack

      for (i = 0; i<np; i++) {
	  TParticle *  iparticle = (TParticle *) fParticles->At(i);
	  Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
	  if (pSelected[i]) {
	      
	      kf   = iparticle->GetPdgCode();	      
	      ks   = iparticle->GetStatusCode();	      
	      
	      p[0] = iparticle->Px();
	      p[1] = iparticle->Py();
	      p[2] = iparticle->Pz();
	      origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
	      origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
	      origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
		    
	      tof = kconv*iparticle->T();
	      
	      imo = -1;
	      TParticle* mother = 0;
	      if (hasMother) {
		  imo = iparticle->GetFirstMother();
		  mother = (TParticle *) fParticles->At(imo);
		  imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
	      } // if has mother   

	      Bool_t tFlag = (fTrackIt && (ks == 1));
	      
	      PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
	      KeepTrack(nt);
	      newPos[i] = nt;
	  } // if selected
      } // particle loop
      delete[] newPos;
      delete[] pSelected;
      if (nc>0) {
	  jev += nc;
	  if (jev >= fNpart || fNpart == -1) {
	      break;
	  }
      }
  } // event loop
  MakeHeader();
  SetHighWaterMark(nt);
}

//______________________________________________________________________________
Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
{
//
// Looks recursively if one of the daughters has been selected
//
//    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
    Int_t imin = -1;
    Int_t imax = -1;
    Int_t i;
    Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
    Bool_t selected = kFALSE;
    if (hasDaughters) {
	imin = iparticle->GetFirstDaughter();
	imax = iparticle->GetLastDaughter();       
	for (i = imin; i <= imax; i++){
	    TParticle *  jparticle = (TParticle *) fParticles->At(i);	
	    Int_t ip = jparticle->GetPdgCode();
	    if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
		selected=kTRUE; break;
	    }
	    if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
	}
    } else {
	return kFALSE;
    }
    return selected;
}


//______________________________________________________________________________
Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
{
// Select flavor of particle
// 0: all
// 4: charm and beauty
// 5: beauty
    Bool_t res = 0;
    
    if (fFlavor == 0) {
	res = kTRUE;
    } else {
	Int_t ifl = TMath::Abs(pid/100);
	if (ifl > 10) ifl/=10;
	res = (fFlavor == ifl);
    }
//
//  This part if gamma writing is inhibited
    if (fNoGammas) 
	res = res && (pid != kGamma && pid != kPi0);
//
    return res;
}

//______________________________________________________________________________
Bool_t AliGenDPMjet::Stable(TParticle*  particle)
{
// Return true for a stable particle
//
    
//    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
    if (particle->GetStatusCode() == 1) return kTRUE;
    else return kFALSE;

}

//______________________________________________________________________________
void AliGenDPMjet::MakeHeader()
{
// Builds the event header, to be called after each event
    AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
    ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
    ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
    ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
    ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
    							 fDPMjet->GetfIt());
 ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
// Bookkeeping for kinematic bias
    ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
// Event Vertex
    header->SetPrimaryVertex(fVertex);
    gAlice->SetGenEventHeader(header);    
 AddHeader(header);
}

void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
{
    // Add header to container or runloader
    if (fContainer) {
        fContainer->AddHeader(header);
    } else {
        AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
    }
}


//______________________________________________________________________________
AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
{
// Assignment operator
    return *this;
}


//______________________________________________________________________________
Commit	Line	Data
1a52e0ed	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16
	17	// Generator using DPMJET as an external generator
	18	// The main DPMJET options are accessable for the user through this interface.
	19	// Uses the TDPMjet implementation of TGenerator.
	20
	21	#include <TDPMjet.h>
	22	#include <TRandom.h>
	23	#include <TArrayI.h>
	24	#include <TParticle.h>
	25	#include <TGraph.h>
	26	#include <TDatabasePDG.h>
	27	#include <TParticlePDG.h>
	28	#include <TParticleClassPDG.h>
	29	#include <TPDGCode.h>
	30	#include <TLorentzVector.h>
9dcc0beb	31	#include "AliRunLoader.h"
1a52e0ed	32	#include "AliGenDPMjet.h"
1a52e0ed	33	#include "AliGenDPMjetEventHeader.h"
1a52e0ed	34	#include "AliRun.h"
f2ba02a4	35	#include "AliDpmJetRndm.h"
9dcc0beb	36	#include "AliHeader.h"
	37	#include "AliStack.h"
	38	#include "AliMC.h"
1a52e0ed	39
f97d4a8e	40	ClassImp(AliGenDPMjet)
1a52e0ed	41
	42	//______________________________________________________________________________
	43	AliGenDPMjet::AliGenDPMjet()
f97d4a8e	44	:AliGenMC(),
	45	fBeamEn(2750.),
	46	fEnergyCMS(5500.),
	47	fMinImpactParam(0.),
	48	fMaxImpactParam(5.),
	49	fICentr(0),
	50	fSelectAll(0),
	51	fFlavor(0),
	52	fTrials(0),
	53	fSpectators(1),
	54	fSpecn(0),
	55	fSpecp(0),
	56	fDPMjet(0),
	57	fParticles(0),
	58	fNoGammas(0),
	59	fLHC(0),
	60	fPi0Decay(0),
	61	fGenImpPar(0.),
	62	fProcess(kDpmMb)
1a52e0ed	63	{
1a52e0ed	64	// Constructor
f2ba02a4	65	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed	66	}
	67
	68
	69	//______________________________________________________________________________
	70	AliGenDPMjet::AliGenDPMjet(Int_t npart)
f97d4a8e	71	:AliGenMC(npart),
	72	fBeamEn(2750.),
	73	fEnergyCMS(5500.),
	74	fMinImpactParam(0.),
	75	fMaxImpactParam(5.),
	76	fICentr(0),
	77	fSelectAll(0),
	78	fFlavor(0),
	79	fTrials(0),
	80	fSpectators(1),
	81	fSpecn(0),
	82	fSpecp(0),
	83	fDPMjet(0),
	84	fParticles(new TClonesArray("TParticle",10000)),
	85	fNoGammas(0),
	86	fLHC(0),
	87	fPi0Decay(0),
	88	fGenImpPar(0.),
	89	fProcess(kDpmMb)
1a52e0ed	90	{
	91	// Default PbPb collisions at 5. 5 TeV
	92	//
	93	fName = "DPMJET";
	94	fTitle= "Particle Generator using DPMJET";
1a52e0ed	95	SetTarget();
1a52e0ed	96	SetProjectile();
f2ba02a4	97	fVertex.Set(3);
f2ba02a4	98	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed	99	}
1a52e0ed	100
ba758f5a	101	AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/Dpmjet/)
f97d4a8e	102	:AliGenMC(),
	103	fBeamEn(2750.),
	104	fEnergyCMS(5500.),
	105	fMinImpactParam(0.),
	106	fMaxImpactParam(5.),
	107	fICentr(0),
	108	fSelectAll(0),
	109	fFlavor(0),
	110	fTrials(0),
	111	fSpectators(1),
	112	fSpecn(0),
	113	fSpecp(0),
	114	fDPMjet(0),
	115	fParticles(0),
	116	fNoGammas(0),
	117	fLHC(0),
	118	fPi0Decay(0),
	119	fGenImpPar(0.),
	120	fProcess(kDpmMb)
ba758f5a	121	{
f97d4a8e	122	// Dummy copy constructor
ba758f5a	123	}
1a52e0ed	124
	125	//______________________________________________________________________________
	126	AliGenDPMjet::~AliGenDPMjet()
	127	{
	128	// Destructor
	129	delete fParticles;
	130	}
1a52e0ed	131	//______________________________________________________________________________
	132	void AliGenDPMjet::Init()
	133	{
	134	// Initialization
	135
e67e0fb7	136	SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
1a52e0ed	137	fBeamEn,fEnergyCMS));
1a52e0ed	138
35c87559	139	fDPMjet=(TDPMjet*) fMCEvGen;
1a52e0ed	140	//
	141	// **** Flag to force central production
	142	// fICentr=1. central production forced
	143	// fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
	144	// fICentr<-99 -> fraction of x-sec. = XSFRAC
	145	// fICentr=-1. -> evaporation/fzc suppressed
	146	// fICentr<-1. -> evaporation/fzc suppressed
ba758f5a	147	if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
ba758f5a	148
1a52e0ed	149	fDPMjet->SetfFCentr(fICentr);
e67e0fb7	150	fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
f97d4a8e	151	fDPMjet->SetPi0Decay(fPi0Decay);
1a52e0ed	152	//
	153	// Initialize DPMjet
	154	//
	155	fDPMjet->Initialize();
1a52e0ed	156	}
	157
	158
	159	//______________________________________________________________________________
	160	void AliGenDPMjet::Generate()
	161	{
	162	// Generate one event
	163
	164	Float_t polar[3] = {0,0,0};
	165	Float_t origin[3] = {0,0,0};
	166	Float_t origin0[3] = {0,0,0};
e67e0fb7	167	Float_t p[3];
1a52e0ed	168	Float_t tof;
	169
	170	// converts from mm/c to s
	171	const Float_t kconv = 0.001/2.999792458e8;
	172	Int_t nt = 0;
	173	Int_t jev = 0;
e67e0fb7	174	Int_t kf, ks, imo;
1a52e0ed	175	kf = 0;
1a52e0ed	176	fTrials = 0;
e67e0fb7	177	// Set collision vertex position
e67e0fb7	178	if (fVertexSmear == kPerEvent) Vertex();
1a52e0ed	179
	180	while(1)
	181	{
	182	// Generate one event
	183	// --------------------------------------------------------------------------
	184	fSpecn = 0;
	185	fSpecp = 0;
	186	// --------------------------------------------------------------------------
	187	fDPMjet->GenerateEvent();
	188	fTrials++;
	189
	190	fDPMjet->ImportParticles(fParticles,"All");
	191	if (fLHC) Boost();
	192
	193	// Temporaneo
	194	fGenImpPar = fDPMjet->GetBImpac();
	195
	196	Int_t np = fParticles->GetEntriesFast();
	197	printf("\n **************************************************%d\n",np);
	198	Int_t nc=0;
	199	if (np==0) continue;
	200	Int_t i;
	201	Int_t* newPos = new Int_t[np];
	202	Int_t* pSelected = new Int_t[np];
	203
	204	for (i = 0; i<np; i++) {
	205	newPos[i] = i;
	206	pSelected[i] = 0;
	207	}
	208
	209	// Get event vertex
	210
f2ba02a4	211	fVertex[0] = origin0[0];
	212	fVertex[1] = origin0[1];
	213	fVertex[2] = origin0[2];
1a52e0ed	214
	215	// First select parent particles
	216
	217	for (i = 0; i<np; i++) {
	218	TParticle iparticle = (TParticle ) fParticles->At(i);
1a52e0ed	219
	220	// Is this a parent particle ?
	221
	222	if (Stable(iparticle)) continue;
	223
	224	Bool_t selected = kTRUE;
	225	Bool_t hasSelectedDaughters = kFALSE;
	226
	227	kf = iparticle->GetPdgCode();
	228	if (kf == 92) continue;
	229	ks = iparticle->GetStatusCode();
	230	// No initial state partons
	231	if (ks==21) continue;
	232
	233	if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
	234	SelectFlavor(kf);
	235	hasSelectedDaughters = DaughtersSelection(iparticle);
	236
	237	// Put particle on the stack if it is either selected or
	238	// it is the mother of at least one seleted particle
	239
	240	if (selected \|\| hasSelectedDaughters) {
	241	nc++;
	242	pSelected[i] = 1;
	243	} // selected
	244	} // particle loop parents
	245
	246	// Now select the final state particles
	247
	248
	249	for (i=0; i<np; i++) {
	250	TParticle iparticle = (TParticle ) fParticles->At(i);
1a52e0ed	251
	252	// Is this a final state particle ?
	253
	254	if (!Stable(iparticle)) continue;
	255
	256	Bool_t selected = kTRUE;
	257	kf = iparticle->GetPdgCode();
	258	ks = iparticle->GetStatusCode();
	259
	260	// --------------------------------------------------------------------------
	261	// Count spectator neutrons and protons (ks == 13, 14)
	262	if(ks == 13 \|\| ks == 14){
	263	if(kf == kNeutron) fSpecn += 1;
	264	if(kf == kProton) fSpecp += 1;
	265	}
	266	// --------------------------------------------------------------------------
	267
	268	if (!fSelectAll) {
	269	selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
	270	if (!fSpectators && selected) selected = (ks == 13 \|\| ks == 14);
	271	}
	272
	273	// Put particle on the stack if selected
	274
	275	if (selected) {
	276	nc++;
	277	pSelected[i] = 1;
	278	} // selected
	279	} // particle loop final state
	280
	281	// Write particles to stack
	282
	283	for (i = 0; i<np; i++) {
	284	TParticle * iparticle = (TParticle *) fParticles->At(i);
	285	Bool_t hasMother = (iparticle->GetFirstMother()>=0);
1a52e0ed	286	if (pSelected[i]) {
	287
	288	kf = iparticle->GetPdgCode();
1a52e0ed	289	ks = iparticle->GetStatusCode();
	290
	291	p[0] = iparticle->Px();
	292	p[1] = iparticle->Py();
	293	p[2] = iparticle->Pz();
e67e0fb7	294	origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
	295	origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
	296	origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
	297
1a52e0ed	298	tof = kconv*iparticle->T();
	299
	300	imo = -1;
	301	TParticle* mother = 0;
	302	if (hasMother) {
	303	imo = iparticle->GetFirstMother();
	304	mother = (TParticle *) fParticles->At(imo);
	305	imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
	306	} // if has mother
	307
1a52e0ed	308	Bool_t tFlag = (fTrackIt && (ks == 1));
1a52e0ed	309
642f15cf	310	PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
1a52e0ed	311	KeepTrack(nt);
	312	newPos[i] = nt;
	313	} // if selected
	314	} // particle loop
	315	delete[] newPos;
	316	delete[] pSelected;
1a52e0ed	317	if (nc>0) {
	318	jev += nc;
	319	if (jev >= fNpart \|\| fNpart == -1) {
1a52e0ed	320	break;
1a52e0ed	321	}
1a52e0ed	322	}
	323	} // event loop
	324	MakeHeader();
	325	SetHighWaterMark(nt);
	326	}
	327
1a52e0ed	328	//______________________________________________________________________________
	329	Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
	330	{
	331	//
	332	// Looks recursively if one of the daughters has been selected
	333	//
	334	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	335	Int_t imin = -1;
	336	Int_t imax = -1;
	337	Int_t i;
	338	Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
	339	Bool_t selected = kFALSE;
	340	if (hasDaughters) {
	341	imin = iparticle->GetFirstDaughter();
	342	imax = iparticle->GetLastDaughter();
	343	for (i = imin; i <= imax; i++){
	344	TParticle * jparticle = (TParticle *) fParticles->At(i);
	345	Int_t ip = jparticle->GetPdgCode();
	346	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	347	selected=kTRUE; break;
	348	}
	349	if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
	350	}
	351	} else {
	352	return kFALSE;
	353	}
	354	return selected;
	355	}
	356
	357
	358
	359	//______________________________________________________________________________
	360	Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
	361	{
	362	// Select flavor of particle
	363	// 0: all
	364	// 4: charm and beauty
	365	// 5: beauty
	366	Bool_t res = 0;
	367
	368	if (fFlavor == 0) {
	369	res = kTRUE;
	370	} else {
	371	Int_t ifl = TMath::Abs(pid/100);
	372	if (ifl > 10) ifl/=10;
	373	res = (fFlavor == ifl);
	374	}
	375	//
	376	// This part if gamma writing is inhibited
	377	if (fNoGammas)
	378	res = res && (pid != kGamma && pid != kPi0);
	379	//
	380	return res;
	381	}
	382
	383	//______________________________________________________________________________
	384	Bool_t AliGenDPMjet::Stable(TParticle* particle)
	385	{
	386	// Return true for a stable particle
	387	//
	388
	389	// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
	390	if (particle->GetStatusCode() == 1) return kTRUE;
	391	else return kFALSE;
392
393	}
394
1a52e0ed	395	//______________________________________________________________________________
	396	void AliGenDPMjet::MakeHeader()
	397	{
	398	// Builds the event header, to be called after each event
	399	AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
	400	((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
	401	((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
	402	((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
1a52e0ed	403	((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
1a52e0ed	404	fDPMjet->GetfIt());
9dcc0beb	405	((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
1a52e0ed	406	// Bookkeeping for kinematic bias
	407	((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
	408	// Event Vertex
f2ba02a4	409	header->SetPrimaryVertex(fVertex);
1a52e0ed	410	gAlice->SetGenEventHeader(header);
9dcc0beb	411	AddHeader(header);
1a52e0ed	412	}
1a52e0ed	413
9dcc0beb	414	void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
	415	{
	416	// Add header to container or runloader
	417	if (fContainer) {
	418	fContainer->AddHeader(header);
	419	} else {
	420	AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
	421	}
	422	}
1a52e0ed	423
	424
	425	//______________________________________________________________________________
ba758f5a	426	AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /rhs/)
1a52e0ed	427	{
	428	// Assignment operator
	429	return *this;
	430	}
	431
	432
	433
	434	//______________________________________________________________________________