[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),
     fTriggerMultiplicity(0),
     fTriggerMultiplicityEta(0),
     fTriggerMultiplicityPtMin(0)
{
// 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),
     fTriggerMultiplicity(0),
     fTriggerMultiplicityEta(0),
     fTriggerMultiplicityPtMin(0)
{
// 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),
     fTriggerMultiplicity(0),
     fTriggerMultiplicityEta(0),
     fTriggerMultiplicityPtMin(0)
{
    // 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();
      //
      // Multiplicity Trigger
      if (fTriggerMultiplicity > 0) {
	Int_t multiplicity = 0;
	for (Int_t i = 0; i < np; i++) {
	  TParticle *  iparticle = (TParticle *) fParticles.At(i);
	
	  Int_t statusCode = iparticle->GetStatusCode();
	
	  // Initial state particle
	  if (statusCode != 1)
	    continue;
	  // eta cut
	  if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
	    continue;
	  // pt cut
	  if (iparticle->Pt() < fTriggerMultiplicityPtMin) 
	    continue;
	  
	  TParticlePDG* pdgPart = iparticle->GetPDG();
	  if (pdgPart && pdgPart->Charge() == 0)
	    continue;
	  ++multiplicity;
	}
	//
	//
	if (multiplicity < fTriggerMultiplicity) continue;
	Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
      }    

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