[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";

    SetBeamEnergy();
    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(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();
    //if (fEvaluate) EvaluateCrossSections();

    //  Initialize AliIonPDGCodes to define ion PDG codes
    /*AliIonPDGCodes *PDGcodes = new AliIonPDGCodes();
    PDGcodes->AddParticlesToPdgDataBase();
    PDGcodes->MapPDGGEant3Codes();*/
}


//______________________________________________________________________________
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], random[6];
  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 j, kf, ks, imo;
  kf = 0;
        
  fTrials = 0;
  for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
  if(fVertexSmear == kPerEvent) {
      Float_t dv[3];
      dv[2] = 1.e10;
      while(TMath::Abs(dv[2])>fCutVertexZ*fOsigma[2]) {
	  Rndm(random,6);
	  for (j=0; j < 3; j++) {
	      dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
		  TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
	  }
      }
      for (j=0; j < 3; j++) origin0[j] += dv[j];
  }
  
  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
      
      //TParticle *iparticle = (TParticle*) fParticles->At(0);
      fVertex[0] = origin0[0];
      fVertex[1] = origin0[1];	
      fVertex[2] = origin0[2];
      //for(int jj=0; jj<3; jj++) printf("	fEventVertex[%d] = %f\n",jj,fEventVertex[jj]);
      
//      First select parent particles

      for (i = 0; i<np; i++) {
	  TParticle *iparticle = (TParticle *) fParticles->At(i);
          //if(!iparticle) iparticle->Print("");

// 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);
	  //if(!iparticle) printf("	i = %d -> iparticle==0\n",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);
//	  Bool_t  hasDaughter = (iparticle->GetFirstDaughter()>=0);

	  if (pSelected[i]) {
	      
	      kf   = iparticle->GetPdgCode();	      
	      if(kf==99999) continue;
	      ks   = iparticle->GetStatusCode();	      
	      
	      p[0] = iparticle->Px();
	      p[1] = iparticle->Py();
	      p[2] = iparticle->Pz();
	      origin[0] = origin0[0]+iparticle->Vx()/10;
	      origin[1] = origin0[1]+iparticle->Vy()/10;
	      origin[2] = origin0[2]+iparticle->Vz()/10;
	      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   

	      // --- Nucleons and photons from nucleus evaporation
	      //if(ks==-1 && (kf==2212 || kf==2112 || kf==22)) imo=-1;
	      
	      // --- Offset added to nuclei, alpha particles, deuteron 
	      // --- and triton PDG codes (according to TGeant3 PDG codes)
	      //if(ks==1000 || ks==1001) kf += 10000000;
	      if(kf>10000) kf += 10000000;
	      
	      if(kf>10000 || ks==1000 || ks==1001) printf(" kf = %d, ks = %d imo = %d \n", kf,ks,imo);
	     

	      Bool_t tFlag = (fTrackIt && (ks == 1));
	      
/*	      // --- Particle NOT to be tracked:
	      // --- (1) chains from valence quark system (idhkk=99999)
	      // --- (2) quark, diquarks and gluons (idhkk=1->8,...,idhkk=21)
	      if(kf==99999 || kf==1 || kf==2 || kf==3 || kf==4 || kf==5 ||
	         kf==6 || kf==7 || kf==8 || kf==21 || 
		 kf==1103 || kf==2101 || kf==2103 || kf==2203 || kf==3101 || 
		 kf==3103 || kf==3201 || kf==3203 || kf==3303 || kf==4101 || 
		 kf==4103 || kf==4201 || kf==4203 || kf==4301 || kf==4303 || 
		 kf==4403 || kf==5101 || kf==5103 || kf==5201 || kf==5203 || 
		 kf==5301 || kf==5303 || kf==5401 || kf==5403 || kf==5503) 
		 tFlag=kFALSE;
*/	      
	      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;
      
      printf("\n I've put %i particles on the stack \n",nc);
      if (nc>0) {
	  jev += nc;
	  if (jev >= fNpart || fNpart == -1) {
	      printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
	      break;
	  }
	  printf("\n Not breaking ### Trials: %i %i %i\n",fTrials, fNpart, jev);
      }
  } // 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::Boost(TClonesArray* particles)
{
//
// Boost cms into LHC lab frame
//
    Double_t dy    = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) / 
				      (Double_t(fZTarget)    * Double_t(fAProjectile)));
    Double_t beta  = TMath::TanH(dy);
    Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
    Double_t gb    = gamma * beta;

    printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma);
    
    Int_t i;
    Int_t np = particles->GetEntriesFast();
    for (i = 0; i < np; i++) 
    {
	TParticle* iparticle = (TParticle*) fParticles->At(i);

	Double_t e   = iparticle->Energy();
	Double_t px  = iparticle->Px();
	Double_t py  = iparticle->Py();
	Double_t pz  = iparticle->Pz();

	Double_t eb  = gamma * e -      gb * pz;
	Double_t pzb =   -gb * e +   gamma * pz;

	iparticle->SetMomentum(px, py, pzb, eb);
    }
}*/


//______________________________________________________________________________
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)->SetWounded(fDPMjet->GetProjWounded(),
//    						    fDPMjet->GetTargWounded());
    ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
    							 fDPMjet->GetfIt());
    /*((AliGenDPMjetEventHeader*) header)->SetCollisions(fDPMjet->GetN0(),
						       fDPMjet->GetN01(),
						       fDPMjet->GetN10(),
						       fDPMjet->GetN11());*/
    ((AliGenDPMjetEventHeader*) header)->SetSpectators(fSpecn, fSpecp);

// 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
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>
	31
	32	#include "AliGenDPMjet.h"
	33	#include "AliGenDPMjetEventHeader.h"
	34	#include "AliPythia.h"
	35	#include "AliRun.h"
f2ba02a4	36	#include "AliDpmJetRndm.h"
1a52e0ed	37
	38	ClassImp(AliGenDPMjet)
	39
	40
	41	//______________________________________________________________________________
	42	AliGenDPMjet::AliGenDPMjet()
	43	:AliGenMC()
	44	{
	45	// Constructor
	46	fParticles = 0;
	47	fDPMjet = 0;
f2ba02a4	48	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed	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	SetBeamEnergy();
	62	SetEnergyCMS();
	63	SetTarget();
	64	SetProjectile();
	65	SetCentral();
	66	SetImpactParameterRange();
	67	SetBoostLHC();
	68
	69	fKeep = 0;
	70	fDecaysOff = 1;
	71	fEvaluate = 0;
	72	fSelectAll = 0;
	73	fFlavor = 0;
	74	fSpectators = 1;
f2ba02a4	75	fVertex.Set(3);
1a52e0ed	76
	77	fParticles = new TClonesArray("TParticle",10000);
	78
	79	// Set random number generator
1a52e0ed	80	fDPMjet = 0;
	81	// Instance AliPythia
	82	AliPythia::Instance();
f2ba02a4	83	AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed	84	}
1a52e0ed	85
ba758f5a	86	AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/Dpmjet/)
	87	:AliGenMC()
	88	{
	89	}
1a52e0ed	90
	91	//______________________________________________________________________________
	92	AliGenDPMjet::~AliGenDPMjet()
	93	{
	94	// Destructor
	95	delete fParticles;
	96	}
	97
	98
	99	//______________________________________________________________________________
	100	void AliGenDPMjet::Init()
	101	{
	102	// Initialization
	103
	104	SetMC(new TDPMjet(fAProjectile, fZProjectile, fATarget, fZTarget,
	105	fBeamEn,fEnergyCMS));
	106
35c87559	107	fDPMjet=(TDPMjet*) fMCEvGen;
1a52e0ed	108	//
	109	// **** Flag to force central production
	110	// fICentr=1. central production forced
	111	// fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
	112	// fICentr<-99 -> fraction of x-sec. = XSFRAC
	113	// fICentr=-1. -> evaporation/fzc suppressed
	114	// fICentr<-1. -> evaporation/fzc suppressed
ba758f5a	115	if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
ba758f5a	116
1a52e0ed	117	fDPMjet->SetfFCentr(fICentr);
	118	fDPMjet->SetbRange(fMinImpactParam,fMaxImpactParam);
	119
	120	//
	121	// Initialize DPMjet
	122	//
	123	fDPMjet->Initialize();
	124	//if (fEvaluate) EvaluateCrossSections();
	125
	126	// Initialize AliIonPDGCodes to define ion PDG codes
	127	/AliIonPDGCodes PDGcodes = new AliIonPDGCodes();
	128	PDGcodes->AddParticlesToPdgDataBase();
	129	PDGcodes->MapPDGGEant3Codes();*/
	130	}
	131
	132
	133	//______________________________________________________________________________
	134	void AliGenDPMjet::Generate()
	135	{
	136	// Generate one event
	137
	138	Float_t polar[3] = {0,0,0};
	139	Float_t origin[3] = {0,0,0};
	140	Float_t origin0[3] = {0,0,0};
	141	Float_t p[3], random[6];
	142	Float_t tof;
	143
	144	// converts from mm/c to s
	145	const Float_t kconv = 0.001/2.999792458e8;
	146	Int_t nt = 0;
	147	Int_t jev = 0;
	148	Int_t j, kf, ks, imo;
	149	kf = 0;
	150
	151	fTrials = 0;
	152	for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
	153	if(fVertexSmear == kPerEvent) {
	154	Float_t dv[3];
	155	dv[2] = 1.e10;
	156	while(TMath::Abs(dv[2])>fCutVertexZ*fOsigma[2]) {
	157	Rndm(random,6);
	158	for (j=0; j < 3; j++) {
	159	dv[j] = fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	160	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	161	}
	162	}
	163	for (j=0; j < 3; j++) origin0[j] += dv[j];
	164	}
	165
	166	while(1)
	167	{
	168	// Generate one event
	169	// --------------------------------------------------------------------------
	170	fSpecn = 0;
	171	fSpecp = 0;
	172	// --------------------------------------------------------------------------
	173	fDPMjet->GenerateEvent();
	174	fTrials++;
	175
	176	fDPMjet->ImportParticles(fParticles,"All");
	177	if (fLHC) Boost();
	178
	179	// Temporaneo
	180	fGenImpPar = fDPMjet->GetBImpac();
181
182	Int_t np = fParticles->GetEntriesFast();
183	printf("\n **************************************************%d\n",np);
184	Int_t nc=0;
185	if (np==0) continue;
186	Int_t i;
187	Int_t* newPos = new Int_t[np];
188	Int_t* pSelected = new Int_t[np];
189
190	for (i = 0; i<np; i++) {
191	newPos[i] = i;
192	pSelected[i] = 0;
193	}
194
195	// Get event vertex
196
197	//TParticle iparticle = (TParticle) fParticles->At(0);
f2ba02a4	198	fVertex[0] = origin0[0];
	199	fVertex[1] = origin0[1];
	200	fVertex[2] = origin0[2];
1a52e0ed	201	//for(int jj=0; jj<3; jj++) printf(" fEventVertex[%d] = %f\n",jj,fEventVertex[jj]);
	202
	203	// First select parent particles
	204
	205	for (i = 0; i<np; i++) {
	206	TParticle iparticle = (TParticle ) fParticles->At(i);
	207	//if(!iparticle) iparticle->Print("");
	208
	209	// Is this a parent particle ?
	210
	211	if (Stable(iparticle)) continue;
	212
	213	Bool_t selected = kTRUE;
	214	Bool_t hasSelectedDaughters = kFALSE;
	215
	216	kf = iparticle->GetPdgCode();
	217	if (kf == 92) continue;
	218	ks = iparticle->GetStatusCode();
	219	// No initial state partons
	220	if (ks==21) continue;
	221
	222	if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
	223	SelectFlavor(kf);
	224	hasSelectedDaughters = DaughtersSelection(iparticle);
	225
	226	// Put particle on the stack if it is either selected or
	227	// it is the mother of at least one seleted particle
	228
	229	if (selected \|\| hasSelectedDaughters) {
	230	nc++;
	231	pSelected[i] = 1;
	232	} // selected
	233	} // particle loop parents
	234
	235	// Now select the final state particles
	236
	237
	238	for (i=0; i<np; i++) {
	239	TParticle iparticle = (TParticle ) fParticles->At(i);
	240	//if(!iparticle) printf(" i = %d -> iparticle==0\n",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);
ba758f5a	276	// Bool_t hasDaughter = (iparticle->GetFirstDaughter()>=0);
1a52e0ed	277
	278	if (pSelected[i]) {
	279
	280	kf = iparticle->GetPdgCode();
	281	if(kf==99999) continue;
	282	ks = iparticle->GetStatusCode();
	283
	284	p[0] = iparticle->Px();
	285	p[1] = iparticle->Py();
	286	p[2] = iparticle->Pz();
	287	origin[0] = origin0[0]+iparticle->Vx()/10;
	288	origin[1] = origin0[1]+iparticle->Vy()/10;
	289	origin[2] = origin0[2]+iparticle->Vz()/10;
	290	tof = kconv*iparticle->T();
	291
	292	imo = -1;
	293	TParticle* mother = 0;
	294	if (hasMother) {
	295	imo = iparticle->GetFirstMother();
	296	mother = (TParticle *) fParticles->At(imo);
	297	imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
	298	} // if has mother
	299
	300	// --- Nucleons and photons from nucleus evaporation
	301	//if(ks==-1 && (kf==2212 \|\| kf==2112 \|\| kf==22)) imo=-1;
	302
	303	// --- Offset added to nuclei, alpha particles, deuteron
	304	// --- and triton PDG codes (according to TGeant3 PDG codes)
	305	//if(ks==1000 \|\| ks==1001) kf += 10000000;
	306	if(kf>10000) kf += 10000000;
	307
	308	if(kf>10000 \|\| ks==1000 \|\| ks==1001) printf(" kf = %d, ks = %d imo = %d \n", kf,ks,imo);
	309
	310
	311	Bool_t tFlag = (fTrackIt && (ks == 1));
	312
	313	/* // --- Particle NOT to be tracked:
	314	// --- (1) chains from valence quark system (idhkk=99999)
	315	// --- (2) quark, diquarks and gluons (idhkk=1->8,...,idhkk=21)
	316	if(kf==99999 \|\| kf==1 \|\| kf==2 \|\| kf==3 \|\| kf==4 \|\| kf==5 \|\|
	317	kf==6 \|\| kf==7 \|\| kf==8 \|\| kf==21 \|\|
	318	kf==1103 \|\| kf==2101 \|\| kf==2103 \|\| kf==2203 \|\| kf==3101 \|\|
	319	kf==3103 \|\| kf==3201 \|\| kf==3203 \|\| kf==3303 \|\| kf==4101 \|\|
	320	kf==4103 \|\| kf==4201 \|\| kf==4203 \|\| kf==4301 \|\| kf==4303 \|\|
	321	kf==4403 \|\| kf==5101 \|\| kf==5103 \|\| kf==5201 \|\| kf==5203 \|\|
	322	kf==5301 \|\| kf==5303 \|\| kf==5401 \|\| kf==5403 \|\| kf==5503)
	323	tFlag=kFALSE;
	324	*/
642f15cf	325	PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
1a52e0ed	326	KeepTrack(nt);
	327	newPos[i] = nt;
	328	} // if selected
	329	} // particle loop
	330	delete[] newPos;
	331	delete[] pSelected;
	332
	333	printf("\n I've put %i particles on the stack \n",nc);
	334	if (nc>0) {
	335	jev += nc;
	336	if (jev >= fNpart \|\| fNpart == -1) {
	337	printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
	338	break;
	339	}
	340	printf("\n Not breaking ### Trials: %i %i %i\n",fTrials, fNpart, jev);
	341	}
	342	} // event loop
	343	MakeHeader();
	344	SetHighWaterMark(nt);
	345	}
	346
	347
	348	//______________________________________________________________________________
	349	void AliGenDPMjet::KeepFullEvent()
	350	{
	351	fKeep=1;
	352	}
	353
	354
	355	//______________________________________________________________________________
	356	/*void AliGenDPMjet::EvaluateCrossSections()
	357	{
	358	// Glauber Calculation of geometrical x-section
	359	//
	360	Float_t xTot = 0.; // barn
	361	Float_t xTotHard = 0.; // barn
	362	Float_t xPart = 0.; // barn
	363	Float_t xPartHard = 0.; // barn
	364	Float_t sigmaHard = 0.1; // mbarn
	365	Float_t bMin = 0.;
	366	Float_t bMax = fDPMjet->GetProjRadius()+fDPMjet->GetTargRadius();
	367	const Float_t kdib = 0.2;
	368	Int_t kMax = Int_t((bMax-bMin)/kdib)+1;
	369
	370
	371	printf("\n Projectile Radius (fm): %f \n",fDPMjet->GetProjRadius());
	372	printf("\n Target Radius (fm): %f \n",fDPMjet->GetTargRadius());
	373	Int_t i;
	374	Float_t oldvalue= 0.;
	375
	376	Float_t* b = new Float_t[kMax];
	377	Float_t* si1 = new Float_t[kMax];
	378	Float_t* si2 = new Float_t[kMax];
	379
	380	for (i = 0; i < kMax; i++)
	381	{
	382	Float_t xb = bMin+i*kdib;
	383	Float_t ov;
	384	ov=fDPMjet->Profile(xb);
	385	// ATT!->Manca la x-sec anel. nucleone-nucleone
	386	Float_t gb = 2.0.01fDPMjet->TMath::Pi()kdibxb(1.-TMath::Exp(-fDPMjet->GetXSFrac()ov));
	387	Float_t gbh = 2.0.01fDPMjet->TMath::Pi()kdibxbsigmaHardov;
	388	xTot+=gb;
	389	xTotHard += gbh;
390	if (xb > fMinImpactParam && xb < fMaxImpactParam)
391	{
392	xPart += gb;
393	xPartHard += gbh;
394	}
395
396	if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
397	oldvalue = xTot;
398	printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
399	printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
400	if (i>0) {
401	si1[i] = gb/kdib;
402	si2[i] = gbh/gb;
403	b[i] = xb;
404	}
405	}
406
407	printf("\n Total cross section (barn): %f \n",xTot);
408	printf("\n Hard cross section (barn): %f \n \n",xTotHard);
409	printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
410	printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
411
412	// Store result as a graph
413	b[0] = 0;
414	si1[0] = 0;
415	si2[0]=si2[1];
416
417	fDsigmaDb = new TGraph(i, b, si1);
418	fDnDb = new TGraph(i, b, si2);
419	}*/
420
421
422	//______________________________________________________________________________
423	Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
424	{
425	//
426	// Looks recursively if one of the daughters has been selected
427	//
428	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
429	Int_t imin = -1;
430	Int_t imax = -1;
431	Int_t i;
432	Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
433	Bool_t selected = kFALSE;
434	if (hasDaughters) {
435	imin = iparticle->GetFirstDaughter();
436	imax = iparticle->GetLastDaughter();
437	for (i = imin; i <= imax; i++){
438	TParticle * jparticle = (TParticle *) fParticles->At(i);
439	Int_t ip = jparticle->GetPdgCode();
440	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
441	selected=kTRUE; break;
442	}
443	if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
444	}
445	} else {
446	return kFALSE;
447	}
448	return selected;
449	}
450
451
452
453	//______________________________________________________________________________
454	Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
455	{
456	// Select flavor of particle
457	// 0: all
458	// 4: charm and beauty
459	// 5: beauty
460	Bool_t res = 0;
461
462	if (fFlavor == 0) {
463	res = kTRUE;
464	} else {
465	Int_t ifl = TMath::Abs(pid/100);
466	if (ifl > 10) ifl/=10;
467	res = (fFlavor == ifl);
468	}
469	//
470	// This part if gamma writing is inhibited
471	if (fNoGammas)
472	res = res && (pid != kGamma && pid != kPi0);
473	//
474	return res;
475	}
476
477	//______________________________________________________________________________
478	Bool_t AliGenDPMjet::Stable(TParticle* particle)
479	{
480	// Return true for a stable particle
481	//
482
483	// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
484	if (particle->GetStatusCode() == 1) return kTRUE;
485	else return kFALSE;
486
487	}
488
489
490
491	//______________________________________________________________________________
492	/void AliGenDPMjet::Boost(TClonesArray particles)
493	{
494	//
495	// Boost cms into LHC lab frame
496	//
497	Double_t dy = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) /
498	(Double_t(fZTarget) * Double_t(fAProjectile)));
499	Double_t beta = TMath::TanH(dy);
500	Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
501	Double_t gb = gamma * beta;
502
503	printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma);
504
505	Int_t i;
506	Int_t np = particles->GetEntriesFast();
507	for (i = 0; i < np; i++)
508	{
509	TParticle* iparticle = (TParticle*) fParticles->At(i);
510
511	Double_t e = iparticle->Energy();
512	Double_t px = iparticle->Px();
513	Double_t py = iparticle->Py();
514	Double_t pz = iparticle->Pz();
515
516	Double_t eb = gamma * e - gb * pz;
517	Double_t pzb = -gb * e + gamma * pz;
518
519	iparticle->SetMomentum(px, py, pzb, eb);
520	}
521	}*/
522
523
524
525	//______________________________________________________________________________
526	void AliGenDPMjet::MakeHeader()
527	{
528	// Builds the event header, to be called after each event
529	AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
530	((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
531	((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
532	((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
533	// ((AliGenDPMjetEventHeader*) header)->SetWounded(fDPMjet->GetProjWounded(),
534	// fDPMjet->GetTargWounded());
535	((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
536	fDPMjet->GetfIt());
537	/((AliGenDPMjetEventHeader) header)->SetCollisions(fDPMjet->GetN0(),
538	fDPMjet->GetN01(),
539	fDPMjet->GetN10(),
540	fDPMjet->GetN11());*/
541	((AliGenDPMjetEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
542
543	// Bookkeeping for kinematic bias
544	((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
545	// Event Vertex
f2ba02a4	546	header->SetPrimaryVertex(fVertex);
1a52e0ed	547	gAlice->SetGenEventHeader(header);
	548	}
	549
	550
	551
	552	//______________________________________________________________________________
ba758f5a	553	AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /rhs/)
1a52e0ed	554	{
	555	// Assignment operator
	556	return *this;
	557	}
	558
	559
	560
	561	//______________________________________________________________________________