[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 "AliGenDPMjet.h"
#include "AliGenDPMjetEventHeader.h"
#include "AliPythia.h"
#include "AliRun.h"
#include "AliDpmJetRndm.h"

 ClassImp(AliGenDPMjet)


//______________________________________________________________________________
AliGenDPMjet::AliGenDPMjet()
                 :AliGenMC()
{
// Constructor
    fParticles = 0;
    fDPMjet = 0;
    AliDpmJetRndm::SetDpmJetRandom(GetRandom());
}


//______________________________________________________________________________
AliGenDPMjet::AliGenDPMjet(Int_t npart)
    :AliGenMC(npart)
{
// Default PbPb collisions at 5. 5 TeV
//
    fName = "DPMJET";
    fTitle= "Particle Generator using DPMJET";

    SetEnergyCMS();
    SetTarget();
    SetProjectile();
    SetCentral();
    SetImpactParameterRange();
    SetBoostLHC();
    
    fKeep       =  0;
    fDecaysOff  =  1;
    fEvaluate   =  0;
    fSelectAll  =  0;
    fFlavor     =  0;
    fSpectators =  1;
    fVertex.Set(3);
        
    fParticles = new TClonesArray("TParticle",10000);    

    // Set random number generator   
    fDPMjet = 0;
    // Instance AliPythia
    AliPythia::Instance(); 
    AliDpmJetRndm::SetDpmJetRandom(GetRandom());
}

AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
:AliGenMC()
{
}

//______________________________________________________________________________
AliGenDPMjet::~AliGenDPMjet()
{
// Destructor
    delete fParticles;
}


//______________________________________________________________________________
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); 
    
//
//  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 origin0[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;
      }
      
//      Get event vertex
      
      fVertex[0] = origin0[0];
      fVertex[1] = origin0[1];	
      fVertex[2] = origin0[2];
      
//      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) ? imo = 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);
}


//______________________________________________________________________________
void AliGenDPMjet::KeepFullEvent()
{
    fKeep=1;
}


//______________________________________________________________________________
/*void AliGenDPMjet::EvaluateCrossSections()
{
//     Glauber Calculation of geometrical x-section
//
    Float_t xTot       = 0.;          // barn
    Float_t xTotHard   = 0.;          // barn 
    Float_t xPart      = 0.;          // barn
    Float_t xPartHard  = 0.;          // barn 
    Float_t sigmaHard  = 0.1;         // mbarn
    Float_t bMin       = 0.;
    Float_t bMax       = fDPMjet->GetProjRadius()+fDPMjet->GetTargRadius();
    const Float_t kdib = 0.2;
    Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;


    printf("\n Projectile Radius (fm): %f \n",fDPMjet->GetProjRadius());
    printf("\n Target     Radius (fm): %f \n",fDPMjet->GetTargRadius());    
    Int_t i;
    Float_t oldvalue= 0.;

    Float_t* b   = new Float_t[kMax];
    Float_t* si1 = new Float_t[kMax];    
    Float_t* si2 = new Float_t[kMax];    
    
    for (i = 0; i < kMax; i++)
    {
	Float_t xb  = bMin+i*kdib;
	Float_t ov;
	ov=fDPMjet->Profile(xb);
	// ATT!->Manca la x-sec anel. nucleone-nucleone
	Float_t gb  =  2.*0.01*fDPMjet->TMath::Pi()*kdib*xb*(1.-TMath::Exp(-fDPMjet->GetXSFrac()*ov));
	Float_t gbh =  2.*0.01*fDPMjet->TMath::Pi()*kdib*xb*sigmaHard*ov;
	xTot+=gb;
	xTotHard += gbh;
	if (xb > fMinImpactParam && xb < fMaxImpactParam)
	{
	    xPart += gb;
	    xPartHard += gbh;
	}
	
	if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
	oldvalue = xTot;
	printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
	printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
	if (i>0) {
	    si1[i] = gb/kdib;
	    si2[i] = gbh/gb;
	    b[i]  = xb;
	}
    }

    printf("\n Total cross section (barn): %f \n",xTot);
    printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
    printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
    printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);

//  Store result as a graph
    b[0] = 0;
    si1[0] = 0;
    si2[0]=si2[1];
    
    fDsigmaDb  = new TGraph(i, b, si1);
    fDnDb      = new TGraph(i, b, si2);
}*/


//______________________________________________________________________________
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());

// Bookkeeping for kinematic bias
    ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
// Event Vertex
    header->SetPrimaryVertex(fVertex);
    gAlice->SetGenEventHeader(header);    
}


//______________________________________________________________________________
AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
{
// Assignment operator
    return *this;
}


//______________________________________________________________________________
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
	32	#include "AliGenDPMjet.h"
	33	#include "AliGenDPMjetEventHeader.h"
	34	#include "AliPythia.h"
	35	#include "AliRun.h"
	36	#include "AliDpmJetRndm.h"
	37
	38	ClassImp(AliGenDPMjet)
	39
	40
	41	//______________________________________________________________________________
	42	AliGenDPMjet::AliGenDPMjet()
	43	:AliGenMC()
	44	{
	45	// Constructor
	46	fParticles = 0;
	47	fDPMjet = 0;
	48	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
	49	}
	50
	51
	52	//______________________________________________________________________________
	53	AliGenDPMjet::AliGenDPMjet(Int_t npart)
	54	:AliGenMC(npart)
	55	{
	56	// Default PbPb collisions at 5. 5 TeV
	57	//
	58	fName = "DPMJET";
	59	fTitle= "Particle Generator using DPMJET";
	60
	61	SetEnergyCMS();
	62	SetTarget();
	63	SetProjectile();
	64	SetCentral();
	65	SetImpactParameterRange();
	66	SetBoostLHC();
	67
	68	fKeep = 0;
	69	fDecaysOff = 1;
	70	fEvaluate = 0;
	71	fSelectAll = 0;
	72	fFlavor = 0;
	73	fSpectators = 1;
	74	fVertex.Set(3);
	75
	76	fParticles = new TClonesArray("TParticle",10000);
	77
	78	// Set random number generator
	79	fDPMjet = 0;
	80	// Instance AliPythia
	81	AliPythia::Instance();
	82	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
	83	}
	84
	85	AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/Dpmjet/)
	86	:AliGenMC()
	87	{
	88	}
	89
	90	//______________________________________________________________________________
	91	AliGenDPMjet::~AliGenDPMjet()
	92	{
	93	// Destructor
	94	delete fParticles;
	95	}
	96
	97
	98	//______________________________________________________________________________
	99	void AliGenDPMjet::Init()
	100	{
	101	// Initialization
	102
	103	SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
	104	fBeamEn,fEnergyCMS));
	105
	106	fDPMjet=(TDPMjet*) fMCEvGen;
	107	//
	108	// **** Flag to force central production
	109	// fICentr=1. central production forced
	110	// fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
	111	// fICentr<-99 -> fraction of x-sec. = XSFRAC
	112	// fICentr=-1. -> evaporation/fzc suppressed
	113	// fICentr<-1. -> evaporation/fzc suppressed
	114	if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
	115
	116	fDPMjet->SetfFCentr(fICentr);
	117	fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
	118
	119	//
	120	// Initialize DPMjet
	121	//
	122	fDPMjet->Initialize();
	123	}
	124
	125
	126	//______________________________________________________________________________
	127	void AliGenDPMjet::Generate()
	128	{
	129	// Generate one event
	130
	131	Float_t polar[3] = {0,0,0};
	132	Float_t origin[3] = {0,0,0};
	133	Float_t origin0[3] = {0,0,0};
	134	Float_t p[3];
	135	Float_t tof;
	136
	137	// converts from mm/c to s
	138	const Float_t kconv = 0.001/2.999792458e8;
	139	Int_t nt = 0;
	140	Int_t jev = 0;
	141	Int_t kf, ks, imo;
	142	kf = 0;
	143	fTrials = 0;
	144	// Set collision vertex position
	145	if (fVertexSmear == kPerEvent) Vertex();
	146
	147	while(1)
	148	{
	149	// Generate one event
	150	// --------------------------------------------------------------------------
	151	fSpecn = 0;
	152	fSpecp = 0;
	153	// --------------------------------------------------------------------------
	154	fDPMjet->GenerateEvent();
	155	fTrials++;
	156
	157	fDPMjet->ImportParticles(fParticles,"All");
	158	if (fLHC) Boost();
	159
	160	// Temporaneo
	161	fGenImpPar = fDPMjet->GetBImpac();
	162
	163	Int_t np = fParticles->GetEntriesFast();
	164	printf("\n **************************************************%d\n",np);
	165	Int_t nc=0;
	166	if (np==0) continue;
	167	Int_t i;
	168	Int_t* newPos = new Int_t[np];
	169	Int_t* pSelected = new Int_t[np];
	170
	171	for (i = 0; i<np; i++) {
	172	newPos[i] = i;
	173	pSelected[i] = 0;
	174	}
	175
	176	// Get event vertex
	177
	178	fVertex[0] = origin0[0];
	179	fVertex[1] = origin0[1];
	180	fVertex[2] = origin0[2];
	181
	182	// First select parent particles
	183
	184	for (i = 0; i<np; i++) {
	185	TParticle iparticle = (TParticle ) fParticles->At(i);
	186
	187	// Is this a parent particle ?
	188
	189	if (Stable(iparticle)) continue;
	190
	191	Bool_t selected = kTRUE;
	192	Bool_t hasSelectedDaughters = kFALSE;
	193
	194	kf = iparticle->GetPdgCode();
	195	if (kf == 92) continue;
	196	ks = iparticle->GetStatusCode();
	197	// No initial state partons
	198	if (ks==21) continue;
	199
	200	if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
	201	SelectFlavor(kf);
	202	hasSelectedDaughters = DaughtersSelection(iparticle);
	203
	204	// Put particle on the stack if it is either selected or
	205	// it is the mother of at least one seleted particle
	206
	207	if (selected \|\| hasSelectedDaughters) {
	208	nc++;
	209	pSelected[i] = 1;
	210	} // selected
	211	} // particle loop parents
	212
	213	// Now select the final state particles
	214
	215
	216	for (i=0; i<np; i++) {
	217	TParticle iparticle = (TParticle ) fParticles->At(i);
	218
	219	// Is this a final state particle ?
	220
	221	if (!Stable(iparticle)) continue;
	222
	223	Bool_t selected = kTRUE;
	224	kf = iparticle->GetPdgCode();
	225	ks = iparticle->GetStatusCode();
	226
	227	// --------------------------------------------------------------------------
	228	// Count spectator neutrons and protons (ks == 13, 14)
	229	if(ks == 13 \|\| ks == 14){
	230	if(kf == kNeutron) fSpecn += 1;
	231	if(kf == kProton) fSpecp += 1;
	232	}
	233	// --------------------------------------------------------------------------
	234
	235	if (!fSelectAll) {
	236	selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
	237	if (!fSpectators && selected) selected = (ks == 13 \|\| ks == 14);
	238	}
	239
	240	// Put particle on the stack if selected
	241
	242	if (selected) {
	243	nc++;
	244	pSelected[i] = 1;
	245	} // selected
	246	} // particle loop final state
	247
	248	// Write particles to stack
	249
	250	for (i = 0; i<np; i++) {
	251	TParticle * iparticle = (TParticle *) fParticles->At(i);
	252	Bool_t hasMother = (iparticle->GetFirstMother()>=0);
	253	if (pSelected[i]) {
	254
	255	kf = iparticle->GetPdgCode();
	256	ks = iparticle->GetStatusCode();
	257
	258	p[0] = iparticle->Px();
	259	p[1] = iparticle->Py();
	260	p[2] = iparticle->Pz();
	261	origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
	262	origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
	263	origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
	264
	265	tof = kconv*iparticle->T();
	266
	267	imo = -1;
	268	TParticle* mother = 0;
	269	if (hasMother) {
	270	imo = iparticle->GetFirstMother();
	271	mother = (TParticle *) fParticles->At(imo);
	272	imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
	273	} // if has mother
	274
	275	Bool_t tFlag = (fTrackIt && (ks == 1));
	276
	277	PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
	278	KeepTrack(nt);
	279	newPos[i] = nt;
	280	} // if selected
	281	} // particle loop
	282	delete[] newPos;
	283	delete[] pSelected;
	284	if (nc>0) {
	285	jev += nc;
	286	if (jev >= fNpart \|\| fNpart == -1) {
	287	break;
	288	}
	289	}
	290	} // event loop
	291	MakeHeader();
	292	SetHighWaterMark(nt);
	293	}
	294
	295
	296	//______________________________________________________________________________
	297	void AliGenDPMjet::KeepFullEvent()
	298	{
	299	fKeep=1;
	300	}
	301
	302
	303	//______________________________________________________________________________
	304	/*void AliGenDPMjet::EvaluateCrossSections()
	305	{
	306	// Glauber Calculation of geometrical x-section
	307	//
	308	Float_t xTot = 0.; // barn
	309	Float_t xTotHard = 0.; // barn
	310	Float_t xPart = 0.; // barn
	311	Float_t xPartHard = 0.; // barn
	312	Float_t sigmaHard = 0.1; // mbarn
	313	Float_t bMin = 0.;
	314	Float_t bMax = fDPMjet->GetProjRadius()+fDPMjet->GetTargRadius();
	315	const Float_t kdib = 0.2;
	316	Int_t kMax = Int_t((bMax-bMin)/kdib)+1;
	317
	318
	319	printf("\n Projectile Radius (fm): %f \n",fDPMjet->GetProjRadius());
	320	printf("\n Target Radius (fm): %f \n",fDPMjet->GetTargRadius());
	321	Int_t i;
	322	Float_t oldvalue= 0.;
	323
	324	Float_t* b = new Float_t[kMax];
	325	Float_t* si1 = new Float_t[kMax];
	326	Float_t* si2 = new Float_t[kMax];
	327
	328	for (i = 0; i < kMax; i++)
	329	{
	330	Float_t xb = bMin+i*kdib;
	331	Float_t ov;
	332	ov=fDPMjet->Profile(xb);
	333	// ATT!->Manca la x-sec anel. nucleone-nucleone
	334	Float_t gb = 2.0.01fDPMjet->TMath::Pi()kdibxb(1.-TMath::Exp(-fDPMjet->GetXSFrac()ov));
	335	Float_t gbh = 2.0.01fDPMjet->TMath::Pi()kdibxbsigmaHardov;
	336	xTot+=gb;
	337	xTotHard += gbh;
	338	if (xb > fMinImpactParam && xb < fMaxImpactParam)
	339	{
	340	xPart += gb;
	341	xPartHard += gbh;
	342	}
	343
	344	if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
	345	oldvalue = xTot;
	346	printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
	347	printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
	348	if (i>0) {
	349	si1[i] = gb/kdib;
	350	si2[i] = gbh/gb;
	351	b[i] = xb;
	352	}
	353	}
	354
	355	printf("\n Total cross section (barn): %f \n",xTot);
	356	printf("\n Hard cross section (barn): %f \n \n",xTotHard);
	357	printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
	358	printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
	359
	360	// Store result as a graph
	361	b[0] = 0;
	362	si1[0] = 0;
	363	si2[0]=si2[1];
	364
	365	fDsigmaDb = new TGraph(i, b, si1);
	366	fDnDb = new TGraph(i, b, si2);
	367	}*/
	368
	369
	370	//______________________________________________________________________________
	371	Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
	372	{
	373	//
	374	// Looks recursively if one of the daughters has been selected
	375	//
	376	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	377	Int_t imin = -1;
	378	Int_t imax = -1;
	379	Int_t i;
	380	Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
	381	Bool_t selected = kFALSE;
	382	if (hasDaughters) {
	383	imin = iparticle->GetFirstDaughter();
	384	imax = iparticle->GetLastDaughter();
	385	for (i = imin; i <= imax; i++){
	386	TParticle * jparticle = (TParticle *) fParticles->At(i);
	387	Int_t ip = jparticle->GetPdgCode();
	388	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	389	selected=kTRUE; break;
	390	}
	391	if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
	392	}
	393	} else {
	394	return kFALSE;
	395	}
	396	return selected;
	397	}
	398
	399
	400
	401	//______________________________________________________________________________
	402	Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
	403	{
	404	// Select flavor of particle
	405	// 0: all
	406	// 4: charm and beauty
	407	// 5: beauty
	408	Bool_t res = 0;
	409
	410	if (fFlavor == 0) {
	411	res = kTRUE;
	412	} else {
	413	Int_t ifl = TMath::Abs(pid/100);
	414	if (ifl > 10) ifl/=10;
	415	res = (fFlavor == ifl);
	416	}
	417	//
	418	// This part if gamma writing is inhibited
	419	if (fNoGammas)
	420	res = res && (pid != kGamma && pid != kPi0);
	421	//
	422	return res;
	423	}
	424
	425	//______________________________________________________________________________
	426	Bool_t AliGenDPMjet::Stable(TParticle* particle)
	427	{
	428	// Return true for a stable particle
	429	//
	430
	431	// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
	432	if (particle->GetStatusCode() == 1) return kTRUE;
	433	else return kFALSE;
	434
	435	}
	436
	437	//______________________________________________________________________________
	438	void AliGenDPMjet::MakeHeader()
	439	{
	440	// Builds the event header, to be called after each event
	441	AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
	442	((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
	443	((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
	444	((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
	445	((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
	446	fDPMjet->GetfIt());
	447
	448	// Bookkeeping for kinematic bias
	449	((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
	450	// Event Vertex
	451	header->SetPrimaryVertex(fVertex);
	452	gAlice->SetGenEventHeader(header);
	453	}
	454
	455
	456
	457	//______________________________________________________________________________
	458	AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /rhs/)
	459	{
	460	// Assignment operator
	461	return *this;
	462	}
	463
	464
	465
	466	//______________________________________________________________________________