[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.),
     fMinImpactParam(0.),
     fMaxImpactParam(5.),
     fICentr(0),
     fSelectAll(0),
     fFlavor(0),
     fTrials(0),
     fSpectators(1),
     fSpecn(0),
     fSpecp(0),
     fDPMjet(0),
     fNoGammas(0),
     fLHC(0),
     fPi0Decay(1),
     fDecayAll(0),
     fGenImpPar(0.),
     fProcess(kDpmMb)
{
// Constructor
    fEnergyCMS = 5500.;
    AliDpmJetRndm::SetDpmJetRandom(GetRandom());
}


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

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

//______________________________________________________________________________
AliGenDPMjet::~AliGenDPMjet()
{
// Destructor
}
//______________________________________________________________________________
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 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
    
    fDPMjet->SetfFCentr(fICentr);  
    fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
    fDPMjet->SetPi0Decay(fPi0Decay);
    fDPMjet->SetDecayAll(fDecayAll);
//
//  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 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;
      }
      
//      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 || kf == 99999) 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 && mother->GetPdgCode() != 99999) ? 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::Instance()->GetHeader()->SetGenEventHeader(header);
    }
}


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


void AliGenDPMjet::FinishRun()
{
    // Print run statistics
    fDPMjet->Dt_Dtuout();
}

    
//______________________________________________________________________________
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(),
f97d4a8e	46	fBeamEn(2750.),
f97d4a8e	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),
f97d4a8e	57	fNoGammas(0),
f97d4a8e	58	fLHC(0),
9c97cedd	59	fPi0Decay(1),
717765ce	60	fDecayAll(0),
f97d4a8e	61	fGenImpPar(0.),
f97d4a8e	62	fProcess(kDpmMb)
1a52e0ed	63	{
1a52e0ed	64	// Constructor
e7c989e4	65	fEnergyCMS = 5500.;
f2ba02a4	66	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed	67	}
	68
	69
	70	//______________________________________________________________________________
	71	AliGenDPMjet::AliGenDPMjet(Int_t npart)
f97d4a8e	72	:AliGenMC(npart),
f97d4a8e	73	fBeamEn(2750.),
f97d4a8e	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),
f97d4a8e	84	fNoGammas(0),
f97d4a8e	85	fLHC(0),
9c97cedd	86	fPi0Decay(1),
717765ce	87	fDecayAll(0),
f97d4a8e	88	fGenImpPar(0.),
f97d4a8e	89	fProcess(kDpmMb)
1a52e0ed	90	{
	91	// Default PbPb collisions at 5. 5 TeV
	92	//
e7c989e4	93	fEnergyCMS = 5500.;
1a52e0ed	94	fName = "DPMJET";
1a52e0ed	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(),
f97d4a8e	104	fBeamEn(2750.),
f97d4a8e	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),
f97d4a8e	115	fNoGammas(0),
f97d4a8e	116	fLHC(0),
9c97cedd	117	fPi0Decay(1),
717765ce	118	fDecayAll(0),
f97d4a8e	119	fGenImpPar(0.),
f97d4a8e	120	fProcess(kDpmMb)
ba758f5a	121	{
f97d4a8e	122	// Dummy copy constructor
e7c989e4	123	fEnergyCMS = 5500.;
ba758f5a	124	}
1a52e0ed	125
	126	//______________________________________________________________________________
	127	AliGenDPMjet::~AliGenDPMjet()
	128	{
	129	// Destructor
1a52e0ed	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
cd19c617	147	if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
ba758f5a	148
1a52e0ed	149	fDPMjet->SetfFCentr(fICentr);
e67e0fb7	150	fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
f97d4a8e	151	fDPMjet->SetPi0Decay(fPi0Decay);
717765ce	152	fDPMjet->SetDecayAll(fDecayAll);
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};
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();
9c97cedd	188
1a52e0ed	189	fTrials++;
1a52e0ed	190
8507138f	191	fDPMjet->ImportParticles(&fParticles,"All");
1a52e0ed	192	if (fLHC) Boost();
	193
	194	// Temporaneo
	195	fGenImpPar = fDPMjet->GetBImpac();
	196
8507138f	197	Int_t np = fParticles.GetEntriesFast();
1a52e0ed	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
1a52e0ed	210	// First select parent particles
	211
	212	for (i = 0; i<np; i++) {
8507138f	213	TParticle iparticle = (TParticle ) fParticles.At(i);
1a52e0ed	214
	215	// Is this a parent particle ?
	216
	217	if (Stable(iparticle)) continue;
	218
	219	Bool_t selected = kTRUE;
	220	Bool_t hasSelectedDaughters = kFALSE;
	221
	222	kf = iparticle->GetPdgCode();
cd19c617	223	if (kf == 92 \|\| kf == 99999) continue;
1a52e0ed	224	ks = iparticle->GetStatusCode();
	225	// No initial state partons
	226	if (ks==21) continue;
	227
	228	if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
	229	SelectFlavor(kf);
9c97cedd	230
9c97cedd	231
1a52e0ed	232	hasSelectedDaughters = DaughtersSelection(iparticle);
1a52e0ed	233
9c97cedd	234
1a52e0ed	235	// Put particle on the stack if it is either selected or
	236	// it is the mother of at least one seleted particle
	237
	238	if (selected \|\| hasSelectedDaughters) {
	239	nc++;
	240	pSelected[i] = 1;
	241	} // selected
	242	} // particle loop parents
	243
	244	// Now select the final state particles
	245
	246
	247	for (i=0; i<np; i++) {
8507138f	248	TParticle iparticle = (TParticle ) fParticles.At(i);
1a52e0ed	249
	250	// Is this a final state particle ?
	251
	252	if (!Stable(iparticle)) continue;
	253
	254	Bool_t selected = kTRUE;
	255	kf = iparticle->GetPdgCode();
	256	ks = iparticle->GetStatusCode();
	257
	258	// --------------------------------------------------------------------------
	259	// Count spectator neutrons and protons (ks == 13, 14)
	260	if(ks == 13 \|\| ks == 14){
	261	if(kf == kNeutron) fSpecn += 1;
	262	if(kf == kProton) fSpecp += 1;
	263	}
	264	// --------------------------------------------------------------------------
	265
	266	if (!fSelectAll) {
	267	selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
	268	if (!fSpectators && selected) selected = (ks == 13 \|\| ks == 14);
	269	}
	270
	271	// Put particle on the stack if selected
	272
	273	if (selected) {
	274	nc++;
	275	pSelected[i] = 1;
	276	} // selected
	277	} // particle loop final state
	278
	279	// Write particles to stack
	280
	281	for (i = 0; i<np; i++) {
8507138f	282	TParticle * iparticle = (TParticle *) fParticles.At(i);
1a52e0ed	283	Bool_t hasMother = (iparticle->GetFirstMother()>=0);
1a52e0ed	284	if (pSelected[i]) {
	285
	286	kf = iparticle->GetPdgCode();
1a52e0ed	287	ks = iparticle->GetStatusCode();
	288
	289	p[0] = iparticle->Px();
	290	p[1] = iparticle->Py();
	291	p[2] = iparticle->Pz();
e67e0fb7	292	origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
	293	origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
	294	origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
	295
1a52e0ed	296	tof = kconv*iparticle->T();
	297
	298	imo = -1;
	299	TParticle* mother = 0;
	300	if (hasMother) {
	301	imo = iparticle->GetFirstMother();
8507138f	302	mother = (TParticle *) fParticles.At(imo);
cd19c617	303	imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
1a52e0ed	304	} // if has mother
1a52e0ed	305
9c97cedd	306
9c97cedd	307
1a52e0ed	308	Bool_t tFlag = (fTrackIt && (ks == 1));
642f15cf	309	PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
1a52e0ed	310	KeepTrack(nt);
	311	newPos[i] = nt;
	312	} // if selected
	313	} // particle loop
	314	delete[] newPos;
	315	delete[] pSelected;
1a52e0ed	316	if (nc>0) {
	317	jev += nc;
	318	if (jev >= fNpart \|\| fNpart == -1) {
1a52e0ed	319	break;
1a52e0ed	320	}
1a52e0ed	321	}
	322	} // event loop
	323	MakeHeader();
	324	SetHighWaterMark(nt);
	325	}
	326
1a52e0ed	327	//______________________________________________________________________________
	328	Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
	329	{
	330	//
	331	// Looks recursively if one of the daughters has been selected
	332	//
	333	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	334	Int_t imin = -1;
	335	Int_t imax = -1;
	336	Int_t i;
	337	Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
	338	Bool_t selected = kFALSE;
	339	if (hasDaughters) {
	340	imin = iparticle->GetFirstDaughter();
	341	imax = iparticle->GetLastDaughter();
	342	for (i = imin; i <= imax; i++){
8507138f	343	TParticle * jparticle = (TParticle *) fParticles.At(i);
1a52e0ed	344	Int_t ip = jparticle->GetPdgCode();
	345	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	346	selected=kTRUE; break;
	347	}
	348	if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
	349	}
	350	} else {
	351	return kFALSE;
	352	}
	353	return selected;
	354	}
	355
	356
	357
	358	//______________________________________________________________________________
	359	Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
	360	{
	361	// Select flavor of particle
	362	// 0: all
	363	// 4: charm and beauty
	364	// 5: beauty
	365	Bool_t res = 0;
	366
	367	if (fFlavor == 0) {
	368	res = kTRUE;
	369	} else {
	370	Int_t ifl = TMath::Abs(pid/100);
	371	if (ifl > 10) ifl/=10;
	372	res = (fFlavor == ifl);
	373	}
	374	//
	375	// This part if gamma writing is inhibited
	376	if (fNoGammas)
	377	res = res && (pid != kGamma && pid != kPi0);
	378	//
	379	return res;
	380	}
	381
	382	//______________________________________________________________________________
	383	Bool_t AliGenDPMjet::Stable(TParticle* particle)
	384	{
	385	// Return true for a stable particle
	386	//
	387
	388	// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
	389	if (particle->GetStatusCode() == 1) return kTRUE;
	390	else return kFALSE;
	391
	392	}
	393
1a52e0ed	394	//______________________________________________________________________________
	395	void AliGenDPMjet::MakeHeader()
	396	{
	397	// Builds the event header, to be called after each event
	398	AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
	399	((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
	400	((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
	401	((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
1a52e0ed	402	((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
1a52e0ed	403	fDPMjet->GetfIt());
9dcc0beb	404	((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
1a52e0ed	405	// Bookkeeping for kinematic bias
	406	((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
	407	// Event Vertex
f2ba02a4	408	header->SetPrimaryVertex(fVertex);
1a52e0ed	409	gAlice->SetGenEventHeader(header);
d2593e7a	410	AddHeader(header);
1a52e0ed	411	}
1a52e0ed	412
9dcc0beb	413	void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
	414	{
	415	// Add header to container or runloader
	416	if (fContainer) {
	417	fContainer->AddHeader(header);
	418	} else {
33c3c91a	419	AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
9dcc0beb	420	}
9dcc0beb	421	}
1a52e0ed	422
	423
	424	//______________________________________________________________________________
ba758f5a	425	AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /rhs/)
1a52e0ed	426	{
	427	// Assignment operator
	428	return *this;
	429	}
	430
	431
d2593e7a	432	void AliGenDPMjet::FinishRun()
	433	{
	434	// Print run statistics
	435	fDPMjet->Dt_Dtuout();
	436	}
	437
	438
1a52e0ed	439
1a52e0ed	440	//______________________________________________________________________________