[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

  Double_t polar[3]    =   {0,0,0};
  Double_t origin[3]   =   {0,0,0};
  Double_t p[4]        =   {0};
  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();
	      p[3] = iparticle->Energy();
	      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 = fTime + 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[0], p[1], p[2], p[3], 
			origin[0], origin[1], origin[2], tof,
			polar[0], polar[1], polar[2],
			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);
    header->SetInteractionTime(fTime);
    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.),
3b50c39e	62	fProcess(kDpmMb),
	63	fTriggerMultiplicity(0),
	64	fTriggerMultiplicityEta(0),
	65	fTriggerMultiplicityPtMin(0)
1a52e0ed	66	{
1a52e0ed	67	// Constructor
e7c989e4	68	fEnergyCMS = 5500.;
f2ba02a4	69	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed	70	}
	71
	72
	73	//______________________________________________________________________________
	74	AliGenDPMjet::AliGenDPMjet(Int_t npart)
f97d4a8e	75	:AliGenMC(npart),
f97d4a8e	76	fBeamEn(2750.),
f97d4a8e	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),
f97d4a8e	87	fNoGammas(0),
f97d4a8e	88	fLHC(0),
9c97cedd	89	fPi0Decay(1),
717765ce	90	fDecayAll(0),
f97d4a8e	91	fGenImpPar(0.),
3b50c39e	92	fProcess(kDpmMb),
	93	fTriggerMultiplicity(0),
	94	fTriggerMultiplicityEta(0),
	95	fTriggerMultiplicityPtMin(0)
1a52e0ed	96	{
	97	// Default PbPb collisions at 5. 5 TeV
	98	//
e7c989e4	99	fEnergyCMS = 5500.;
1a52e0ed	100	fName = "DPMJET";
1a52e0ed	101	fTitle= "Particle Generator using DPMJET";
1a52e0ed	102	SetTarget();
1a52e0ed	103	SetProjectile();
f2ba02a4	104	fVertex.Set(3);
f2ba02a4	105	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed	106	}
1a52e0ed	107
ba758f5a	108	AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/Dpmjet/)
f97d4a8e	109	:AliGenMC(),
f97d4a8e	110	fBeamEn(2750.),
f97d4a8e	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),
f97d4a8e	121	fNoGammas(0),
f97d4a8e	122	fLHC(0),
9c97cedd	123	fPi0Decay(1),
717765ce	124	fDecayAll(0),
f97d4a8e	125	fGenImpPar(0.),
3b50c39e	126	fProcess(kDpmMb),
	127	fTriggerMultiplicity(0),
	128	fTriggerMultiplicityEta(0),
	129	fTriggerMultiplicityPtMin(0)
ba758f5a	130	{
f97d4a8e	131	// Dummy copy constructor
e7c989e4	132	fEnergyCMS = 5500.;
ba758f5a	133	}
1a52e0ed	134
	135	//______________________________________________________________________________
	136	AliGenDPMjet::~AliGenDPMjet()
	137	{
	138	// Destructor
1a52e0ed	139	}
1a52e0ed	140	//______________________________________________________________________________
	141	void AliGenDPMjet::Init()
	142	{
	143	// Initialization
	144
e67e0fb7	145	SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
1a52e0ed	146	fBeamEn,fEnergyCMS));
1a52e0ed	147
35c87559	148	fDPMjet=(TDPMjet*) fMCEvGen;
1a52e0ed	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
cd19c617	156	if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
ba758f5a	157
1a52e0ed	158	fDPMjet->SetfFCentr(fICentr);
e67e0fb7	159	fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
f97d4a8e	160	fDPMjet->SetPi0Decay(fPi0Decay);
717765ce	161	fDPMjet->SetDecayAll(fDecayAll);
1a52e0ed	162	//
	163	// Initialize DPMjet
	164	//
	165	fDPMjet->Initialize();
1a52e0ed	166	}
	167
	168
	169	//______________________________________________________________________________
	170	void AliGenDPMjet::Generate()
	171	{
	172	// Generate one event
	173
c0d429ef	174	Double_t polar[3] = {0,0,0};
	175	Double_t origin[3] = {0,0,0};
	176	Double_t p[4] = {0};
1a52e0ed	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;
e67e0fb7	183	Int_t kf, ks, imo;
1a52e0ed	184	kf = 0;
1a52e0ed	185	fTrials = 0;
e67e0fb7	186	// Set collision vertex position
e67e0fb7	187	if (fVertexSmear == kPerEvent) Vertex();
1a52e0ed	188
	189	while(1)
	190	{
	191	// Generate one event
	192	// --------------------------------------------------------------------------
	193	fSpecn = 0;
	194	fSpecp = 0;
	195	// --------------------------------------------------------------------------
	196	fDPMjet->GenerateEvent();
9c97cedd	197
1a52e0ed	198	fTrials++;
1a52e0ed	199
8507138f	200	fDPMjet->ImportParticles(&fParticles,"All");
1a52e0ed	201	if (fLHC) Boost();
	202
	203	// Temporaneo
	204	fGenImpPar = fDPMjet->GetBImpac();
41cc7c18	205
	206	if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
	207
8507138f	208	Int_t np = fParticles.GetEntriesFast();
3b50c39e	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
1a52e0ed	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
1a52e0ed	251	// First select parent particles
	252
	253	for (i = 0; i<np; i++) {
8507138f	254	TParticle iparticle = (TParticle ) fParticles.At(i);
1a52e0ed	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();
cd19c617	264	if (kf == 92 \|\| kf == 99999) continue;
1a52e0ed	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);
9c97cedd	271
9c97cedd	272
1a52e0ed	273	hasSelectedDaughters = DaughtersSelection(iparticle);
1a52e0ed	274
9c97cedd	275
1a52e0ed	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++) {
8507138f	289	TParticle iparticle = (TParticle ) fParticles.At(i);
1a52e0ed	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++) {
8507138f	323	TParticle * iparticle = (TParticle *) fParticles.At(i);
1a52e0ed	324	Bool_t hasMother = (iparticle->GetFirstMother()>=0);
1a52e0ed	325	if (pSelected[i]) {
	326
	327	kf = iparticle->GetPdgCode();
1a52e0ed	328	ks = iparticle->GetStatusCode();
	329
	330	p[0] = iparticle->Px();
	331	p[1] = iparticle->Py();
	332	p[2] = iparticle->Pz();
c0d429ef	333	p[3] = iparticle->Energy();
e67e0fb7	334	origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
	335	origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
	336	origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
	337
21391258	338	tof = fTime + kconv*iparticle->T();
1a52e0ed	339
	340	imo = -1;
	341	TParticle* mother = 0;
	342	if (hasMother) {
	343	imo = iparticle->GetFirstMother();
8507138f	344	mother = (TParticle *) fParticles.At(imo);
cd19c617	345	imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
1a52e0ed	346	} // if has mother
1a52e0ed	347
9c97cedd	348
9c97cedd	349
1a52e0ed	350	Bool_t tFlag = (fTrackIt && (ks == 1));
c0d429ef	351	PushTrack(tFlag, imo, kf,
	352	p[0], p[1], p[2], p[3],
	353	origin[0], origin[1], origin[2], tof,
	354	polar[0], polar[1], polar[2],
	355	kPNoProcess, nt, 1., ks);
1a52e0ed	356	KeepTrack(nt);
	357	newPos[i] = nt;
	358	} // if selected
	359	} // particle loop
	360	delete[] newPos;
	361	delete[] pSelected;
1a52e0ed	362	if (nc>0) {
	363	jev += nc;
	364	if (jev >= fNpart \|\| fNpart == -1) {
1a52e0ed	365	break;
1a52e0ed	366	}
1a52e0ed	367	}
	368	} // event loop
	369	MakeHeader();
	370	SetHighWaterMark(nt);
	371	}
	372
1a52e0ed	373	//______________________________________________________________________________
	374	Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
	375	{
	376	//
	377	// Looks recursively if one of the daughters has been selected
	378	//
	379	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	380	Int_t imin = -1;
	381	Int_t imax = -1;
	382	Int_t i;
	383	Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
	384	Bool_t selected = kFALSE;
	385	if (hasDaughters) {
	386	imin = iparticle->GetFirstDaughter();
	387	imax = iparticle->GetLastDaughter();
	388	for (i = imin; i <= imax; i++){
8507138f	389	TParticle * jparticle = (TParticle *) fParticles.At(i);
1a52e0ed	390	Int_t ip = jparticle->GetPdgCode();
	391	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	392	selected=kTRUE; break;
	393	}
	394	if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
	395	}
	396	} else {
	397	return kFALSE;
	398	}
	399	return selected;
	400	}
	401
	402
	403
	404	//______________________________________________________________________________
	405	Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
	406	{
	407	// Select flavor of particle
	408	// 0: all
	409	// 4: charm and beauty
	410	// 5: beauty
	411	Bool_t res = 0;
	412
	413	if (fFlavor == 0) {
	414	res = kTRUE;
	415	} else {
	416	Int_t ifl = TMath::Abs(pid/100);
	417	if (ifl > 10) ifl/=10;
	418	res = (fFlavor == ifl);
	419	}
	420	//
	421	// This part if gamma writing is inhibited
	422	if (fNoGammas)
	423	res = res && (pid != kGamma && pid != kPi0);
	424	//
	425	return res;
	426	}
	427
	428	//______________________________________________________________________________
	429	Bool_t AliGenDPMjet::Stable(TParticle* particle)
	430	{
	431	// Return true for a stable particle
	432	//
	433
	434	// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
	435	if (particle->GetStatusCode() == 1) return kTRUE;
	436	else return kFALSE;
	437
	438	}
	439
1a52e0ed	440	//______________________________________________________________________________
	441	void AliGenDPMjet::MakeHeader()
	442	{
4eb21a7b	443	printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
1a52e0ed	444	// Builds the event header, to be called after each event
	445	AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
	446	((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
	447	((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
	448	((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
3c349853	449	((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetProjParticipants(),
	450	fDPMjet->GetTargParticipants());
	451	((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
	452	// Bookkeeping for kinematic bias
1a52e0ed	453	((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
3c349853	454	// Event Vertex
f2ba02a4	455	header->SetPrimaryVertex(fVertex);
21391258	456	header->SetInteractionTime(fTime);
1a52e0ed	457	gAlice->SetGenEventHeader(header);
d2593e7a	458	AddHeader(header);
1a52e0ed	459	}
1a52e0ed	460
9dcc0beb	461	void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
	462	{
	463	// Add header to container or runloader
	464	if (fContainer) {
	465	fContainer->AddHeader(header);
	466	} else {
33c3c91a	467	AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
9dcc0beb	468	}
9dcc0beb	469	}
1a52e0ed	470
	471
	472	//______________________________________________________________________________
ba758f5a	473	AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /rhs/)
1a52e0ed	474	{
	475	// Assignment operator
	476	return *this;
	477	}
	478
	479
d2593e7a	480	void AliGenDPMjet::FinishRun()
	481	{
	482	// Print run statistics
	483	fDPMjet->Dt_Dtuout();
	484	}
	485
	486
1a52e0ed	487
1a52e0ed	488	//______________________________________________________________________________