#include <TLorentzVector.h>
#include <TPDGCode.h>
#include <TParticle.h>
-
+#include <TVirtualMC.h>
+#include <TParticlePDG.h>
#include "AliGenHijingEventHeader.h"
#define AliGenAmptEventHeader AliGenHijingEventHeader
#include "AliAmptRndm.h"
#include "AliLog.h"
#include "AliRun.h"
+#include "AliDecayer.h"
ClassImp(AliGenAmpt)
AliGenAmpt::AliGenAmpt()
: AliGenMC(),
+ fDecayer(NULL),
fFrame("CMS"),
fMinImpactParam(0.),
fMaxImpactParam(5.),
fLHC(kFALSE),
fRandomPz(kFALSE),
fNoHeavyQuarks(kFALSE),
- fEventTime(0.),
- fIsoft(1),
+ fIsoft(4),
fNtMax(150),
fIpop(1),
fXmu(3.2264),
fAlpha(1./3),
fStringA(0.5),
fStringB(0.9),
- fHeader(new AliGenAmptEventHeader("Ampt"))
+ fEventTime(0.),
+ fHeader(new AliGenAmptEventHeader("Ampt")),
+ fDecay(kTRUE),
+ fRotating(kFALSE)
{
// Constructor
- fEnergyCMS = 5500.;
+ fEnergyCMS = 2760.;
AliAmptRndm::SetAmptRandom(GetRandom());
}
AliGenAmpt::AliGenAmpt(Int_t npart)
: AliGenMC(npart),
+ fDecayer(NULL),
fFrame("CMS"),
fMinImpactParam(0.),
fMaxImpactParam(5.),
fLHC(kFALSE),
fRandomPz(kFALSE),
fNoHeavyQuarks(kFALSE),
- fEventTime(0.),
fIsoft(1),
fNtMax(150),
fIpop(1),
fAlpha(1./3),
fStringA(0.5),
fStringB(0.9),
- fHeader(new AliGenAmptEventHeader("Ampt"))
+ fEventTime(0.),
+ fHeader(new AliGenAmptEventHeader("Ampt")),
+ fDecay(kTRUE),
+ fRotating(kFALSE)
{
- // Default PbPb collisions at 5.5 TeV
+ // Default PbPb collisions at 2.76 TeV
- fEnergyCMS = 5500.;
+ fEnergyCMS = 2760.;
fName = "Ampt";
fTitle= "Particle Generator using AMPT";
AliAmptRndm::SetAmptRandom(GetRandom());
fAmpt->Initialize();
if (fEvaluate)
EvaluateCrossSections();
+
+ fAmpt->SetReactionPlaneAngle(0.0);
+}
+
+void AliGenAmpt::SetSeed(UInt_t seed)
+{
+ AliAmptRndm::GetAmptRandom()->SetSeed(seed);
}
void AliGenAmpt::Generate()
Float_t polar[3] = {0,0,0};
Float_t origin[3] = {0,0,0};
Float_t origin0[3] = {0,0,0};
+ Float_t time0 = 0.;
Float_t p[3];
Float_t tof;
- // converts from mm/c to s
- const Float_t kconv = 0.001/2.99792458e8;
-
Int_t nt = 0;
Int_t jev = 0;
Int_t j, kf, ks, ksp, imo;
fTrials = 0;
for (j = 0;j < 3; j++)
origin0[j] = fOrigin[j];
+ //time0 = fTimeOrigin;
if(fVertexSmear == kPerEvent) {
Vertex();
for (j=0; j < 3; j++)
origin0[j] = fVertex[j];
+ //time0 = fTime;
}
Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
while(1) {
+
+ // Generate random reaction plane angle if requested
+ if( fRotating ) {
+ TRandom *r=AliAmptRndm::GetAmptRandom();
+ fAmpt->SetReactionPlaneAngle(TMath::TwoPi()*r->Rndm());
+ }
+
// Generate one event
Int_t fpemask = gSystem->GetFPEMask();
gSystem->SetFPEMask(0);
gSystem->SetFPEMask(fpemask);
fTrials++;
fNprimaries = 0;
+
+
fAmpt->ImportParticles(&fParticles,"All");
Int_t np = fParticles.GetEntriesFast();
if (np == 0 )
continue;
-
+ //
+ //RS>>: Decayers now returns cm and sec. Since TAmpt returns mm and mm/c, convert its result to cm and sec here
+ const Float_t kconvT=0.001/2.999792458e8; // mm/c to seconds conversion
+ const Float_t kconvL=1./10; // mm to cm conversion
+ for (int ip=np;ip--;) {
+ TParticle* part = (TParticle*)fParticles[ip];
+ if (!part) continue;
+ part->SetProductionVertex(part->Vx()*kconvL,part->Vy()*kconvL,part->Vz()*kconvL,kconvT*part->T());
+ }
+ // RS<<
+ //
if (fTrigger != kNoTrigger) {
if (!CheckTrigger())
continue;
}
+
+ AliDecayer *decayer = 0;
+ //if (gMC)
+ // decayer = gMC->GetDecayer();
+ decayer = fDecayer; //AMPT does not do the strong decays per dafault
+
+ if (decayer&&fDecay) {
+ TClonesArray arr("TParticle",100);
+ for( Int_t nLoop=0; nLoop!=2; ++nLoop) { // In order to produce more than one generation of decays: NumberOfNestedLoops set to 2
+ Int_t np2 = np;
+ for (Int_t i = 0; i < np; i++) {
+ TParticle *iparticle = (TParticle *)fParticles.At(i);
+ if (!Stable(iparticle)) // true if particle has daughters already
+ continue;
+ kf = TMath::Abs(iparticle->GetPdgCode());
+ if (kf==92)
+ continue;
+ if( !IsThisAKnownParticle(iparticle) ) continue; // skip undesired particles
+ /*
+ if (0) { // this turned out to be too cumbersome!
+ if (kf!=331&&kf!=3114&&kf!=3114&&kf!=411&&kf!=-4122&&kf!=-3324&&kf!=-3312&&kf!=-3114&&
+ kf!=-311&&kf!=3214&&kf!=-3214&&kf!=-433&&kf!=413&&kf!=3122&&kf!=-3122&&kf!=-413&&
+ kf!=-421&&kf!=-423&&kf!=3324&&kf!=-313&&kf!=213&&kf!=-213&&kf!=3314&&kf!=3222&&
+ kf!=-3222&&kf!=3224&&kf!=-3224&&kf!=-4212&&kf!=4212&&kf!=433&&kf!=423&&kf!=-3322&&
+ kf!=3322&&kf!=-3314)
+ continue; //decay eta',Sigma*+,Sigma*-,D+,Lambda_c-,Xi*0_bar,Xi-_bar,Sigma*-,
+ // K0_bar,Sigma*0,Sigma*0_bar,D*_s-,D*+,Lambda0,Lambda0_bar,D*-
+ // D0_bar,D*0_bar,Xi*0,K*0_bar,rho+,rho-,Xi*-,Sigma-,
+ // Sigma+,Sigma*+,Sigma*-,Sigma_c-,Sigma_c+,D*_s+,D*0,Xi0_bar
+ // Xi0,Xi*+
+ //} else { // really only decay particles if there are not known to Geant3
+ // if (gMC->IdFromPDG(kf)>0)
+ // continue;
+ }
+ if (0) { // defining the particle for Geant3 leads to a floating point exception.
+ TParticlePDG *pdg = iparticle->GetPDG(1);
+ //pdg->Print(); printf("%s\n",pdg->ParticleClass());
+ TString ptype(pdg->ParticleClass());
+ TMCParticleType mctype(kPTUndefined);
+ if (ptype=="Baryon" || ptype=="Meson")
+ mctype = kPTHadron;
+ gMC->DefineParticle(pdg->PdgCode(), pdg->GetName(), mctype, pdg->Mass(), pdg->Charge(), pdg->Lifetime(),
+ ptype,pdg->Width(), (Int_t)pdg->Spin(), (Int_t)pdg->Parity(), 0,
+ (Int_t)pdg->Isospin(), 0, 0, 0, 0, pdg->Stable());
+ gMC->SetUserDecay(pdg->PdgCode());
+ continue;
+ }
+ */
+ TLorentzVector pmom(iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy());
+ decayer->Decay(kf,&pmom);
+ decayer->ImportParticles(&arr);
+ Int_t ndecayed = arr.GetEntries();
+ if (ndecayed>1) {
+ if (np2+ndecayed>fParticles.GetSize())
+ fParticles.Expand(2*fParticles.GetSize());
+ //arr.Print();
+ // iparticle->SetStatusCode(2); to be compatible with Hijing
+ iparticle->SetFirstDaughter(np2);
+ for (Int_t jj = 1; jj < ndecayed; jj++) {
+ TParticle *jp = (TParticle *)arr.At(jj);
+ if (jp->GetFirstMother()!=1)
+ continue;
+
+ TParticle *newp = new(fParticles[np2]) TParticle(jp->GetPdgCode(),
+ 0, //1, //to be compatible with Hijing
+ i,
+ -1,
+ -1,
+ -1,
+ jp->Px(),jp->Py(),jp->Pz(),jp->Energy(),
+ jp->Vx(),jp->Vy(),jp->Vz(),jp->T());
+ //take care of the phi
+ //if((kf == 333)||(kf == 313)) {
+ if(IsThisAKnownParticle(iparticle)) {
+ //Printf("=============PANOS===================");
+ //Printf("Phi detected - daughet is: %d",jp->GetPdgCode());
+ newp->SetUniqueID(4);
+ }
+ else newp->SetUniqueID( jp->GetStatusCode() );
+ np2++;
+ } // end of jj->nDecayedParticles
+ iparticle->SetLastDaughter(np2-1);
+ } // end of nDecayedPrticles>1
+ } // end of i->np
+ np = fParticles.GetEntries();
+ if (np!=np2) {
+ AliError(Form("Something is fishy: %d %d\n", np,np2));
+ }
+ } // end of nLoop->NumberOfNestedLoops
+ } else {
+ if (fDecay)
+ AliError("No decayer found, but fDecay==kTRUE!");
+ }
+
if (fLHC)
Boost();
fVertex[0] = origin0[0];
fVertex[1] = origin0[1];
fVertex[2] = origin0[2];
+ //fTime = time0;
// First select parent particles
for (Int_t i = 0; i < np; i++) {
TParticle *iparticle = (TParticle *) fParticles.At(i);
// Is this a parent particle ?
- if (Stable(iparticle)) continue;
+ if (Stable(iparticle)) continue; // quit if particle has no daughters
Bool_t selected = kTRUE;
Bool_t hasSelectedDaughters = kFALSE;
- kf = iparticle->GetPdgCode();
- ks = iparticle->GetStatusCode();
+ kf = iparticle->GetPdgCode();
+ ks = iparticle->GetStatusCode();
if (kf == 92)
continue;
for (Int_t i = 0; i<np; i++) {
TParticle *iparticle = (TParticle *) fParticles.At(i);
// Is this a final state particle ?
- if (!Stable(iparticle))
+ if (!Stable(iparticle)) continue; // quit if particle has daughters
+ Bool_t selected = kTRUE;
+ kf = iparticle->GetPdgCode();
+ if (kf == 92)
continue;
- Bool_t selected = kTRUE;
- kf = iparticle->GetPdgCode();
- ks = iparticle->GetStatusCode();
- ksp = iparticle->GetUniqueID();
+ ks = iparticle->GetStatusCode();
+ ksp = iparticle->GetUniqueID();
// --------------------------------------------------------------------------
// Count spectator neutrons and protons
if (selected) {
nc++;
pSelected[i] = 1;
+ if (0) printf("---> %d %d %d %s\n",i,nc,kf,iparticle->GetName());
} // selected
} // particle loop final state
- //Time of the interactions
- Float_t tInt = 0.;
- if (fPileUpTimeWindow > 0.)
- tInt = fPileUpTimeWindow * (2. * gRandom->Rndm() - 1.);
-
// Write particles to stack
for (Int_t 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]) {
+ TParticle *iparticle = (TParticle *) fParticles.At(i);
+ Bool_t hasMother = (iparticle->GetFirstMother() >=0);
+ Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
kf = iparticle->GetPdgCode();
ks = iparticle->GetStatusCode();
p[0] = iparticle->Px();
p[1] = iparticle->Py();
p[2] = iparticle->Pz() * sign;
- origin[0] = origin0[0]+iparticle->Vx()/10;
- origin[1] = origin0[1]+iparticle->Vy()/10;
- origin[2] = origin0[2]+iparticle->Vz()/10;
- fEventTime = 0.;
-
- if (TestBit(kVertexRange)) {
- fEventTime = sign * origin0[2] / 2.99792458e10;
- tof = kconv * iparticle->T() + fEventTime;
- } else {
- tof = kconv * iparticle->T();
- fEventTime = tInt;
- if (fPileUpTimeWindow > 0.) tof += tInt;
- }
+ origin[0] = origin0[0]+iparticle->Vx();
+ origin[1] = origin0[1]+iparticle->Vy();
+ origin[2] = origin0[2]+iparticle->Vz();
+ tof = time0 + iparticle->T();
+
imo = -1;
TParticle* mother = 0;
+ TMCProcess procID = (TMCProcess) iparticle->GetUniqueID();
if (hasMother) {
imo = iparticle->GetFirstMother();
mother = (TParticle *) fParticles.At(imo);
imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
+ } else { // if has no mothers then it was created by AMPT
+ if(procID==999)
+ procID = kPPrimary; // reseting to ALIROOT convention
+ else
+ procID = kPNoProcess; // for expectators
} // if has mother
Bool_t tFlag = (fTrackIt && !hasDaughter);
- PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
+ PushTrack(tFlag,imo,kf,p,origin,polar,tof,procID,nt, 1., ks);
fNprimaries++;
KeepTrack(nt);
newPos[i] = nt;
SetHighWaterMark(nt);
}
+Bool_t AliGenAmpt::IsThisAKnownParticle(TParticle *thisGuy)
+{
+ // In order to prevent AMPT to introduce weird particles into the decayer and transporter
+ // blame cperez@cern.ch for this method
+
+ Int_t pdgcode = TMath::Abs( thisGuy->GetPdgCode() );
+
+ Int_t myFavoriteParticles[ 38] = { 3322, 3314, 3312, 3224, 3222, // Xi0 Xi*+- Xi+- Sigma*-+ Sigma-+
+ 3214, 3212, 3122, 3114, 3112, // Sigma*0 Sigma0 Lambda0 Sigma*+- Sigma+-
+ 2224, 2214, 2212, 2114, 2112, // Delta--++ Delta-+ proton Delta0 neutron
+ 1114, 323, 321, 313, 311, // Delta+- K*-+ K-+ K*0 K0
+ 213, 211, 11, 22, 111, // rho-+ pi-+ e+- gamma pi0
+ 113, 130, 221, 223, 310, // rho0 K_L0 eta omega K_S0
+ 331, 333, 3324, 431, 421, // eta' phi Xi*0 Ds-+ D0
+ 411, 413, 13 // D-+ D*-+ mu+-
+ };
+
+ Bool_t found = kFALSE;
+ for(Int_t i=0; i!=38; ++i)
+ if( myFavoriteParticles[i] == pdgcode ) {
+ found = kTRUE;
+ break;
+ }
+
+ return found;
+}
+
void AliGenAmpt::EvaluateCrossSections()
{
// Glauber Calculation of geometrical x-section
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];
+ Float_t* b = new Float_t[kMax]; memset(b,0,kMax*sizeof(Float_t));
+ Float_t* si1 = new Float_t[kMax]; memset(si1,0,kMax*sizeof(Float_t));
+ Float_t* si2 = new Float_t[kMax]; memset(si2,0,kMax*sizeof(Float_t));
for (i = 0; i < kMax; i++) {
Float_t xb = bMin+i*kdib;
Float_t ov=fAmpt->Profile(xb);
return res;
}
-Bool_t AliGenAmpt::Stable(TParticle* particle) const
+Bool_t AliGenAmpt::Stable(TParticle* particle) const
{
// Return true for a stable particle
+ if (!particle)
+ return kFALSE;
if (particle->GetFirstDaughter() < 0 )
- {
return kTRUE;
- } else {
- return kFALSE;
- }
+ return kFALSE;
+
+ /// ADD LIST
+
}
void AliGenAmpt::MakeHeader()
fAmpt->GetN11());
fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
fTargetSpecn,fTargetSpecp);
- fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
+ //fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
+ fHeader->SetReactionPlaneAngle(fAmpt->GetReactionPlaneAngle());
//printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
// 4-momentum vectors of the triggered jets.
// Event Vertex
fHeader->SetPrimaryVertex(fVertex);
fHeader->SetInteractionTime(fEventTime);
-
+
fCollisionGeometry = fHeader;
AddHeader(fHeader);
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff