[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()
{
// 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.),
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
	174	Float_t polar[3] = {0,0,0};
	175	Float_t origin[3] = {0,0,0};
e67e0fb7	176	Float_t p[3];
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();
e67e0fb7	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
1a52e0ed	337	tof = kconv*iparticle->T();
	338
	339	imo = -1;
	340	TParticle* mother = 0;
	341	if (hasMother) {
	342	imo = iparticle->GetFirstMother();
8507138f	343	mother = (TParticle *) fParticles.At(imo);
cd19c617	344	imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
1a52e0ed	345	} // if has mother
1a52e0ed	346
9c97cedd	347
9c97cedd	348
1a52e0ed	349	Bool_t tFlag = (fTrackIt && (ks == 1));
642f15cf	350	PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
1a52e0ed	351	KeepTrack(nt);
	352	newPos[i] = nt;
	353	} // if selected
	354	} // particle loop
	355	delete[] newPos;
	356	delete[] pSelected;
1a52e0ed	357	if (nc>0) {
	358	jev += nc;
	359	if (jev >= fNpart \|\| fNpart == -1) {
1a52e0ed	360	break;
1a52e0ed	361	}
1a52e0ed	362	}
	363	} // event loop
	364	MakeHeader();
	365	SetHighWaterMark(nt);
	366	}
	367
1a52e0ed	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++){
8507138f	384	TParticle * jparticle = (TParticle *) fParticles.At(i);
1a52e0ed	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
1a52e0ed	435	//______________________________________________________________________________
	436	void AliGenDPMjet::MakeHeader()
	437	{
	438	// Builds the event header, to be called after each event
	439	AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
	440	((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
	441	((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
	442	((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
1a52e0ed	443	((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
1a52e0ed	444	fDPMjet->GetfIt());
9dcc0beb	445	((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
1a52e0ed	446	// Bookkeeping for kinematic bias
	447	((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
	448	// Event Vertex
f2ba02a4	449	header->SetPrimaryVertex(fVertex);
1a52e0ed	450	gAlice->SetGenEventHeader(header);
d2593e7a	451	AddHeader(header);
1a52e0ed	452	}
1a52e0ed	453
9dcc0beb	454	void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
	455	{
	456	// Add header to container or runloader
	457	if (fContainer) {
	458	fContainer->AddHeader(header);
	459	} else {
33c3c91a	460	AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
9dcc0beb	461	}
9dcc0beb	462	}
1a52e0ed	463
	464
	465	//______________________________________________________________________________
ba758f5a	466	AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /rhs/)
1a52e0ed	467	{
	468	// Assignment operator
	469	return *this;
	470	}
	471
	472
d2593e7a	473	void AliGenDPMjet::FinishRun()
	474	{
	475	// Print run statistics
	476	fDPMjet->Dt_Dtuout();
	477	}
	478
	479
1a52e0ed	480
1a52e0ed	481	//______________________________________________________________________________