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