[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();

      if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();

      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()
{
  printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
// 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->GetProjParticipants(), 
    							 fDPMjet->GetTargParticipants());
     ((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	if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
	207
	208	Int_t np = fParticles.GetEntriesFast();
	209	//
	210	// Multiplicity Trigger
	211	if (fTriggerMultiplicity > 0) {
	212	Int_t multiplicity = 0;
	213	for (Int_t i = 0; i < np; i++) {
	214	TParticle * iparticle = (TParticle *) fParticles.At(i);
	215
	216	Int_t statusCode = iparticle->GetStatusCode();
	217
	218	// Initial state particle
	219	if (statusCode != 1)
	220	continue;
	221	// eta cut
	222	if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
	223	continue;
	224	// pt cut
	225	if (iparticle->Pt() < fTriggerMultiplicityPtMin)
	226	continue;
	227
	228	TParticlePDG* pdgPart = iparticle->GetPDG();
	229	if (pdgPart && pdgPart->Charge() == 0)
	230	continue;
	231	++multiplicity;
	232	}
	233	//
	234	//
	235	if (multiplicity < fTriggerMultiplicity) continue;
	236	Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
	237	}
	238
	239	Int_t nc = 0;
	240	if (np == 0) continue;
	241
	242	Int_t i;
	243	Int_t* newPos = new Int_t[np];
	244	Int_t* pSelected = new Int_t[np];
	245
	246	for (i = 0; i<np; i++) {
	247	newPos[i] = i;
	248	pSelected[i] = 0;
	249	}
	250
	251	// First select parent particles
	252
	253	for (i = 0; i<np; i++) {
	254	TParticle iparticle = (TParticle ) fParticles.At(i);
	255
	256	// Is this a parent particle ?
	257
	258	if (Stable(iparticle)) continue;
	259
	260	Bool_t selected = kTRUE;
	261	Bool_t hasSelectedDaughters = kFALSE;
	262
	263	kf = iparticle->GetPdgCode();
	264	if (kf == 92 \|\| kf == 99999) continue;
	265	ks = iparticle->GetStatusCode();
	266	// No initial state partons
	267	if (ks==21) continue;
	268
	269	if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
	270	SelectFlavor(kf);
	271
	272
	273	hasSelectedDaughters = DaughtersSelection(iparticle);
	274
	275
	276	// Put particle on the stack if it is either selected or
	277	// it is the mother of at least one seleted particle
	278
	279	if (selected \|\| hasSelectedDaughters) {
	280	nc++;
	281	pSelected[i] = 1;
	282	} // selected
	283	} // particle loop parents
	284
	285	// Now select the final state particles
	286
	287
	288	for (i=0; i<np; i++) {
	289	TParticle iparticle = (TParticle ) fParticles.At(i);
	290
	291	// Is this a final state particle ?
	292
	293	if (!Stable(iparticle)) continue;
	294
	295	Bool_t selected = kTRUE;
	296	kf = iparticle->GetPdgCode();
	297	ks = iparticle->GetStatusCode();
	298
	299	// --------------------------------------------------------------------------
	300	// Count spectator neutrons and protons (ks == 13, 14)
	301	if(ks == 13 \|\| ks == 14){
	302	if(kf == kNeutron) fSpecn += 1;
	303	if(kf == kProton) fSpecp += 1;
	304	}
	305	// --------------------------------------------------------------------------
	306
	307	if (!fSelectAll) {
	308	selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
	309	if (!fSpectators && selected) selected = (ks == 13 \|\| ks == 14);
	310	}
	311
	312	// Put particle on the stack if selected
	313
	314	if (selected) {
	315	nc++;
	316	pSelected[i] = 1;
	317	} // selected
	318	} // particle loop final state
	319
	320	// Write particles to stack
	321
	322	for (i = 0; i<np; i++) {
	323	TParticle * iparticle = (TParticle *) fParticles.At(i);
	324	Bool_t hasMother = (iparticle->GetFirstMother()>=0);
	325	if (pSelected[i]) {
	326
	327	kf = iparticle->GetPdgCode();
	328	ks = iparticle->GetStatusCode();
	329
	330	p[0] = iparticle->Px();
	331	p[1] = iparticle->Py();
	332	p[2] = iparticle->Pz();
	333	origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
	334	origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
	335	origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
	336
	337	tof = kconv*iparticle->T();
	338
	339	imo = -1;
	340	TParticle* mother = 0;
	341	if (hasMother) {
	342	imo = iparticle->GetFirstMother();
	343	mother = (TParticle *) fParticles.At(imo);
	344	imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
	345	} // if has mother
	346
	347
	348
	349	Bool_t tFlag = (fTrackIt && (ks == 1));
	350	PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
	351	KeepTrack(nt);
	352	newPos[i] = nt;
	353	} // if selected
	354	} // particle loop
	355	delete[] newPos;
	356	delete[] pSelected;
	357	if (nc>0) {
	358	jev += nc;
	359	if (jev >= fNpart \|\| fNpart == -1) {
	360	break;
	361	}
	362	}
	363	} // event loop
	364	MakeHeader();
	365	SetHighWaterMark(nt);
	366	}
	367
	368	//______________________________________________________________________________
	369	Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
	370	{
	371	//
	372	// Looks recursively if one of the daughters has been selected
	373	//
	374	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	375	Int_t imin = -1;
	376	Int_t imax = -1;
	377	Int_t i;
	378	Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
	379	Bool_t selected = kFALSE;
	380	if (hasDaughters) {
	381	imin = iparticle->GetFirstDaughter();
	382	imax = iparticle->GetLastDaughter();
	383	for (i = imin; i <= imax; i++){
	384	TParticle * jparticle = (TParticle *) fParticles.At(i);
	385	Int_t ip = jparticle->GetPdgCode();
	386	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	387	selected=kTRUE; break;
	388	}
	389	if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
	390	}
	391	} else {
	392	return kFALSE;
	393	}
	394	return selected;
	395	}
	396
	397
	398
	399	//______________________________________________________________________________
	400	Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
	401	{
	402	// Select flavor of particle
	403	// 0: all
	404	// 4: charm and beauty
	405	// 5: beauty
	406	Bool_t res = 0;
	407
	408	if (fFlavor == 0) {
	409	res = kTRUE;
	410	} else {
	411	Int_t ifl = TMath::Abs(pid/100);
	412	if (ifl > 10) ifl/=10;
	413	res = (fFlavor == ifl);
	414	}
	415	//
	416	// This part if gamma writing is inhibited
	417	if (fNoGammas)
	418	res = res && (pid != kGamma && pid != kPi0);
	419	//
	420	return res;
	421	}
	422
	423	//______________________________________________________________________________
	424	Bool_t AliGenDPMjet::Stable(TParticle* particle)
	425	{
	426	// Return true for a stable particle
	427	//
	428
	429	// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
	430	if (particle->GetStatusCode() == 1) return kTRUE;
	431	else return kFALSE;
	432
	433	}
	434
	435	//______________________________________________________________________________
	436	void AliGenDPMjet::MakeHeader()
	437	{
	438	printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
	439	// Builds the event header, to be called after each event
	440	AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
	441	((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
	442	((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
	443	((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
	444	((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetProjParticipants(),
	445	fDPMjet->GetTargParticipants());
	446	((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
	447	// Bookkeeping for kinematic bias
	448	((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
	449	// Event Vertex
	450	header->SetPrimaryVertex(fVertex);
	451	gAlice->SetGenEventHeader(header);
	452	AddHeader(header);
	453	}
	454
	455	void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
	456	{
	457	// Add header to container or runloader
	458	if (fContainer) {
	459	fContainer->AddHeader(header);
	460	} else {
	461	AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
	462	}
	463	}
	464
	465
	466	//______________________________________________________________________________
	467	AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /rhs/)
	468	{
	469	// Assignment operator
	470	return *this;
	471	}
	472
	473
	474	void AliGenDPMjet::FinishRun()
	475	{
	476	// Print run statistics
	477	fDPMjet->Dt_Dtuout();
	478	}
	479
	480
	481
	482	//______________________________________________________________________________