[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.),
     fEnergyCMS(5500.),
     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
    AliDpmJetRndm::SetDpmJetRandom(GetRandom());
}


//______________________________________________________________________________
AliGenDPMjet::AliGenDPMjet(Int_t npart)
    :AliGenMC(npart),
     fBeamEn(2750.),
     fEnergyCMS(5500.),
     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
//
    fName = "DPMJET";
    fTitle= "Particle Generator using DPMJET";
    SetTarget();
    SetProjectile();
    fVertex.Set(3);
    AliDpmJetRndm::SetDpmJetRandom(GetRandom());
}

AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
    :AliGenMC(),
     fBeamEn(2750.),
     fEnergyCMS(5500.),
     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
}

//______________________________________________________________________________
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 && 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) 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) ? 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
	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	fEnergyCMS(5500.),
	48	fMinImpactParam(0.),
	49	fMaxImpactParam(5.),
	50	fICentr(0),
	51	fSelectAll(0),
	52	fFlavor(0),
	53	fTrials(0),
	54	fSpectators(1),
	55	fSpecn(0),
	56	fSpecp(0),
	57	fDPMjet(0),
	58	fNoGammas(0),
	59	fLHC(0),
	60	fPi0Decay(0),
	61	fGenImpPar(0.),
	62	fProcess(kDpmMb)
	63	{
	64	// Constructor
	65	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
	66	}
	67
	68
	69	//______________________________________________________________________________
	70	AliGenDPMjet::AliGenDPMjet(Int_t npart)
	71	:AliGenMC(npart),
	72	fBeamEn(2750.),
	73	fEnergyCMS(5500.),
	74	fMinImpactParam(0.),
	75	fMaxImpactParam(5.),
	76	fICentr(0),
	77	fSelectAll(0),
	78	fFlavor(0),
	79	fTrials(0),
	80	fSpectators(1),
	81	fSpecn(0),
	82	fSpecp(0),
	83	fDPMjet(0),
	84	fNoGammas(0),
	85	fLHC(0),
	86	fPi0Decay(0),
	87	fGenImpPar(0.),
	88	fProcess(kDpmMb)
	89	{
	90	// Default PbPb collisions at 5. 5 TeV
	91	//
	92	fName = "DPMJET";
	93	fTitle= "Particle Generator using DPMJET";
	94	SetTarget();
	95	SetProjectile();
	96	fVertex.Set(3);
	97	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
	98	}
	99
	100	AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/Dpmjet/)
	101	:AliGenMC(),
	102	fBeamEn(2750.),
	103	fEnergyCMS(5500.),
	104	fMinImpactParam(0.),
	105	fMaxImpactParam(5.),
	106	fICentr(0),
	107	fSelectAll(0),
	108	fFlavor(0),
	109	fTrials(0),
	110	fSpectators(1),
	111	fSpecn(0),
	112	fSpecp(0),
	113	fDPMjet(0),
	114	fNoGammas(0),
	115	fLHC(0),
	116	fPi0Decay(0),
	117	fGenImpPar(0.),
	118	fProcess(kDpmMb)
	119	{
	120	// Dummy copy constructor
	121	}
	122
	123	//______________________________________________________________________________
	124	AliGenDPMjet::~AliGenDPMjet()
	125	{
	126	// Destructor
	127	}
	128	//______________________________________________________________________________
	129	void AliGenDPMjet::Init()
	130	{
	131	// Initialization
	132
	133	SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
	134	fBeamEn,fEnergyCMS));
	135
	136	fDPMjet=(TDPMjet*) fMCEvGen;
	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
	144	if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
	145
	146	fDPMjet->SetfFCentr(fICentr);
	147	fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
	148	fDPMjet->SetPi0Decay(fPi0Decay);
	149	//
	150	// Initialize DPMjet
	151	//
	152	fDPMjet->Initialize();
	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};
	163	Float_t p[3];
	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;
	170	Int_t kf, ks, imo;
	171	kf = 0;
	172	fTrials = 0;
	173	// Set collision vertex position
	174	if (fVertexSmear == kPerEvent) Vertex();
	175
	176	while(1)
	177	{
	178	// Generate one event
	179	// --------------------------------------------------------------------------
	180	fSpecn = 0;
	181	fSpecp = 0;
	182	// --------------------------------------------------------------------------
	183	fDPMjet->GenerateEvent();
	184	fTrials++;
	185
	186	fDPMjet->ImportParticles(&fParticles,"All");
	187	if (fLHC) Boost();
	188
	189	// Temporaneo
	190	fGenImpPar = fDPMjet->GetBImpac();
	191
	192	Int_t np = fParticles.GetEntriesFast();
	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
	205	// First select parent particles
	206
	207	for (i = 0; i<np; i++) {
	208	TParticle iparticle = (TParticle ) fParticles.At(i);
	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();
	218	if (kf == 92) continue;
	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++) {
	240	TParticle iparticle = (TParticle ) fParticles.At(i);
	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++) {
	274	TParticle * iparticle = (TParticle *) fParticles.At(i);
	275	Bool_t hasMother = (iparticle->GetFirstMother()>=0);
	276	if (pSelected[i]) {
	277
	278	kf = iparticle->GetPdgCode();
	279	ks = iparticle->GetStatusCode();
	280
	281	p[0] = iparticle->Px();
	282	p[1] = iparticle->Py();
	283	p[2] = iparticle->Pz();
	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
	288	tof = kconv*iparticle->T();
	289
	290	imo = -1;
	291	TParticle* mother = 0;
	292	if (hasMother) {
	293	imo = iparticle->GetFirstMother();
	294	mother = (TParticle *) fParticles.At(imo);
	295	imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
	296	} // if has mother
	297
	298	Bool_t tFlag = (fTrackIt && (ks == 1));
	299
	300	PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
	301	KeepTrack(nt);
	302	newPos[i] = nt;
	303	} // if selected
	304	} // particle loop
	305	delete[] newPos;
	306	delete[] pSelected;
	307	if (nc>0) {
	308	jev += nc;
	309	if (jev >= fNpart \|\| fNpart == -1) {
	310	break;
	311	}
	312	}
	313	} // event loop
	314	MakeHeader();
	315	SetHighWaterMark(nt);
	316	}
	317
	318	//______________________________________________________________________________
	319	Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
	320	{
	321	//
	322	// Looks recursively if one of the daughters has been selected
	323	//
	324	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	325	Int_t imin = -1;
	326	Int_t imax = -1;
	327	Int_t i;
	328	Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
	329	Bool_t selected = kFALSE;
	330	if (hasDaughters) {
	331	imin = iparticle->GetFirstDaughter();
	332	imax = iparticle->GetLastDaughter();
	333	for (i = imin; i <= imax; i++){
	334	TParticle * jparticle = (TParticle *) fParticles.At(i);
	335	Int_t ip = jparticle->GetPdgCode();
	336	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	337	selected=kTRUE; break;
	338	}
	339	if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
	340	}
	341	} else {
	342	return kFALSE;
	343	}
	344	return selected;
	345	}
	346
	347
	348
	349	//______________________________________________________________________________
	350	Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
	351	{
	352	// Select flavor of particle
	353	// 0: all
	354	// 4: charm and beauty
	355	// 5: beauty
	356	Bool_t res = 0;
	357
	358	if (fFlavor == 0) {
	359	res = kTRUE;
	360	} else {
	361	Int_t ifl = TMath::Abs(pid/100);
	362	if (ifl > 10) ifl/=10;
	363	res = (fFlavor == ifl);
	364	}
	365	//
	366	// This part if gamma writing is inhibited
	367	if (fNoGammas)
	368	res = res && (pid != kGamma && pid != kPi0);
	369	//
	370	return res;
	371	}
	372
	373	//______________________________________________________________________________
	374	Bool_t AliGenDPMjet::Stable(TParticle* particle)
	375	{
	376	// Return true for a stable particle
	377	//
	378
	379	// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
	380	if (particle->GetStatusCode() == 1) return kTRUE;
	381	else return kFALSE;
	382
	383	}
	384
	385	//______________________________________________________________________________
	386	void AliGenDPMjet::MakeHeader()
	387	{
	388	// Builds the event header, to be called after each event
	389	AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
	390	((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
	391	((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
	392	((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
	393	((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
	394	fDPMjet->GetfIt());
	395	((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
	396	// Bookkeeping for kinematic bias
	397	((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
	398	// Event Vertex
	399	header->SetPrimaryVertex(fVertex);
	400	gAlice->SetGenEventHeader(header);
	401	AddHeader(header);
	402	}
	403
	404	void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
	405	{
	406	// Add header to container or runloader
	407	if (fContainer) {
	408	fContainer->AddHeader(header);
	409	} else {
	410	AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
	411	}
	412	}
	413
	414
	415	//______________________________________________________________________________
	416	AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /rhs/)
	417	{
	418	// Assignment operator
	419	return *this;
	420	}
	421
	422
	423	void AliGenDPMjet::FinishRun()
	424	{
	425	// Print run statistics
	426	fDPMjet->Dt_Dtuout();
	427	}
	428
	429
	430
	431	//______________________________________________________________________________