X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TAmpt%2FAliGenAmpt.cxx;h=0424e0be49c057e59bc07c805adff84320c3a3a0;hb=5a9dc560b8bb65e65c88f8cb8278ee6b9c25e9c5;hp=20b5083a96c81c24cbf37f51a4216cd8f2837296;hpb=a004b331996aee37f6f5c2f354b835b20b19170a;p=u%2Fmrichter%2FAliRoot.git diff --git a/TAmpt/AliGenAmpt.cxx b/TAmpt/AliGenAmpt.cxx index 20b5083a96c..0424e0be49c 100644 --- a/TAmpt/AliGenAmpt.cxx +++ b/TAmpt/AliGenAmpt.cxx @@ -25,17 +25,20 @@ #include #include #include - +#include +#include #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.), @@ -71,23 +74,26 @@ AliGenAmpt::AliGenAmpt() 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.), @@ -123,7 +129,6 @@ AliGenAmpt::AliGenAmpt(Int_t npart) fLHC(kFALSE), fRandomPz(kFALSE), fNoHeavyQuarks(kFALSE), - fEventTime(0.), fIsoft(1), fNtMax(150), fIpop(1), @@ -131,11 +136,14 @@ AliGenAmpt::AliGenAmpt(Int_t npart) 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()); @@ -222,6 +230,9 @@ void AliGenAmpt::Init() fAmpt->Initialize(); if (fEvaluate) EvaluateCrossSections(); + + fAmpt->SetReactionPlaneAngle(0.0); + fRotating=kFALSE; } void AliGenAmpt::Generate() @@ -231,12 +242,10 @@ 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; @@ -245,16 +254,25 @@ void AliGenAmpt::Generate() 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); @@ -262,15 +280,121 @@ void AliGenAmpt::Generate() 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(); @@ -288,17 +412,18 @@ void AliGenAmpt::Generate() 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; @@ -322,12 +447,13 @@ void AliGenAmpt::Generate() for (Int_t i = 0; iGetPdgCode(); + 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 @@ -349,47 +475,41 @@ void AliGenAmpt::Generate() 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; iGetFirstMother() >=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; @@ -412,6 +532,33 @@ void AliGenAmpt::Generate() 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 @@ -431,9 +578,9 @@ void AliGenAmpt::EvaluateCrossSections() 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); @@ -523,16 +670,18 @@ Bool_t AliGenAmpt::SelectFlavor(Int_t pid) 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() @@ -550,7 +699,8 @@ 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. @@ -580,7 +730,7 @@ void AliGenAmpt::MakeHeader() // Event Vertex fHeader->SetPrimaryVertex(fVertex); fHeader->SetInteractionTime(fEventTime); - + fCollisionGeometry = fHeader; AddHeader(fHeader); }