[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(0),
     fGenImpPar(0.),
     fProcess(kDpmMb)
{
// Constructor
    fEnergyCMS = 5500.;
    AliDpmJetRndm::SetDpmJetRandom(GetRandom());
}


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

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

//______________________________________________________________________________
AliGenDPMjet::~AliGenDPMjet()
{
// Destructor
}
//______________________________________________________________________________
void AliGenDPMjet::Init()
{
// Initialization
    
    SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget, 
		      fBeamEn,fEnergyCMS));

    fDPMjet=(TDPMjet*) fMCEvGen;
    //
    // **** Flag to force central production
    // fICentr=1. central production forced 
    // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam	  
    // fICentr<-99 -> fraction of x-sec. = XSFRAC		  
    // fICentr=-1. -> evaporation/fzc suppressed		  
    // fICentr<-1. -> evaporation/fzc suppressed		  
    if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
    
    fDPMjet->SetfFCentr(fICentr);  
    fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
    fDPMjet->SetPi0Decay(fPi0Decay);
//
//  Initialize DPMjet  
//    
    fDPMjet->Initialize();
}


//______________________________________________________________________________
void AliGenDPMjet::Generate()
{
// Generate one event

  Float_t polar[3]    =   {0,0,0};
  Float_t origin[3]   =   {0,0,0};
  Float_t p[3];
  Float_t tof;

//  converts from mm/c to s
  const Float_t kconv = 0.001/2.999792458e8;
  Int_t nt  = 0;
  Int_t jev = 0;
  Int_t kf, ks, imo;
  kf = 0;
  fTrials = 0;
  //  Set collision vertex position 
  if (fVertexSmear == kPerEvent) Vertex();
  
  while(1)
  {
//    Generate one event
// --------------------------------------------------------------------------
      fSpecn = 0;  
      fSpecp = 0;
// --------------------------------------------------------------------------
      fDPMjet->GenerateEvent();
      fTrials++;

      fDPMjet->ImportParticles(&fParticles,"All");      
      if (fLHC) Boost();

      // Temporaneo
      fGenImpPar = fDPMjet->GetBImpac();
      
      Int_t np = fParticles.GetEntriesFast();
      printf("\n **************************************************%d\n",np);
      Int_t nc=0;
      if (np==0) continue;
      Int_t i;
      Int_t* newPos     = new Int_t[np];
      Int_t* pSelected  = new Int_t[np];

      for (i = 0; i<np; i++) {
	  newPos[i]    = i;
	  pSelected[i] = 0;
      }
      
//      First select parent particles

      for (i = 0; i<np; i++) {
	  TParticle *iparticle = (TParticle *) fParticles.At(i);

// Is this a parent particle ?

	  if (Stable(iparticle)) continue;

	  Bool_t  selected             =  kTRUE;
	  Bool_t  hasSelectedDaughters =  kFALSE;
	  
	  kf = iparticle->GetPdgCode();
	  if (kf == 92 || kf == 99999) continue;
	  ks = iparticle->GetStatusCode();
// No initial state partons
          if (ks==21) continue;
	    
	  if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
			       SelectFlavor(kf);
	  hasSelectedDaughters = DaughtersSelection(iparticle);

// Put particle on the stack if it is either selected or 
// it is the mother of at least one seleted particle

	  if (selected || hasSelectedDaughters) {
	      nc++;
	      pSelected[i] = 1;
	  } // selected
      } // particle loop parents

// Now select the final state particles


      for (i=0; i<np; i++) {
	  TParticle *iparticle = (TParticle *) fParticles.At(i);

// Is this a final state particle ?

	  if (!Stable(iparticle)) continue;
      
	  Bool_t  selected =  kTRUE;
	  kf = iparticle->GetPdgCode();
	  ks = iparticle->GetStatusCode();

// --------------------------------------------------------------------------
// Count spectator neutrons and protons (ks == 13, 14)
	  if(ks == 13 || ks == 14){
	      if(kf == kNeutron) fSpecn += 1;
	      if(kf == kProton)  fSpecp += 1;
	  }
// --------------------------------------------------------------------------

	  if (!fSelectAll) {
	      selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
	      if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
	  }

// Put particle on the stack if selected

	  if (selected) {
	      nc++;
	      pSelected[i] = 1;
	  } // selected
      } // particle loop final state

// Write particles to stack

      for (i = 0; i<np; i++) {
	  TParticle *  iparticle = (TParticle *) fParticles.At(i);
	  Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
	  if (pSelected[i]) {
	      
	      kf   = iparticle->GetPdgCode();	      
	      ks   = iparticle->GetStatusCode();	      
	      
	      p[0] = iparticle->Px();
	      p[1] = iparticle->Py();
	      p[2] = iparticle->Pz();
	      origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
	      origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
	      origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
		    
	      tof = kconv*iparticle->T();
	      
	      imo = -1;
	      TParticle* mother = 0;
	      if (hasMother) {
		  imo = iparticle->GetFirstMother();
		  mother = (TParticle *) fParticles.At(imo);
		  imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
	      } // if has mother   

	      Bool_t tFlag = (fTrackIt && (ks == 1));
	      PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
	      KeepTrack(nt);
	      newPos[i] = nt;
	  } // if selected
      } // particle loop
      delete[] newPos;
      delete[] pSelected;
      if (nc>0) {
	  jev += nc;
	  if (jev >= fNpart || fNpart == -1) {
	      break;
	  }
      }
  } // event loop
  MakeHeader();
  SetHighWaterMark(nt);
}

//______________________________________________________________________________
Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
{
//
// Looks recursively if one of the daughters has been selected
//
//    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
    Int_t imin = -1;
    Int_t imax = -1;
    Int_t i;
    Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
    Bool_t selected = kFALSE;
    if (hasDaughters) {
	imin = iparticle->GetFirstDaughter();
	imax = iparticle->GetLastDaughter();       
	for (i = imin; i <= imax; i++){
	    TParticle *  jparticle = (TParticle *) fParticles.At(i);	
	    Int_t ip = jparticle->GetPdgCode();
	    if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
		selected=kTRUE; break;
	    }
	    if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
	}
    } else {
	return kFALSE;
    }
    return selected;
}


//______________________________________________________________________________
Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
{
// Select flavor of particle
// 0: all
// 4: charm and beauty
// 5: beauty
    Bool_t res = 0;
    
    if (fFlavor == 0) {
	res = kTRUE;
    } else {
	Int_t ifl = TMath::Abs(pid/100);
	if (ifl > 10) ifl/=10;
	res = (fFlavor == ifl);
    }
//
//  This part if gamma writing is inhibited
    if (fNoGammas) 
	res = res && (pid != kGamma && pid != kPi0);
//
    return res;
}

//______________________________________________________________________________
Bool_t AliGenDPMjet::Stable(TParticle*  particle)
{
// Return true for a stable particle
//
    
//    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
    if (particle->GetStatusCode() == 1) return kTRUE;
    else return kFALSE;

}

//______________________________________________________________________________
void AliGenDPMjet::MakeHeader()
{
// Builds the event header, to be called after each event
    AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
    ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
    ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
    ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
    ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
    							 fDPMjet->GetfIt());
 ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
// Bookkeeping for kinematic bias
    ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
// Event Vertex
    header->SetPrimaryVertex(fVertex);
    gAlice->SetGenEventHeader(header);    
    AddHeader(header);
}

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