1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 // Generator using AMPT as an external generator
20 #include "AliGenAmpt.h"
22 #include <TClonesArray.h>
25 #include <TLorentzVector.h>
27 #include <TParticle.h>
28 #include <TVirtualMC.h>
29 #include <TParticlePDG.h>
30 #include "AliGenHijingEventHeader.h"
31 #define AliGenAmptEventHeader AliGenHijingEventHeader
32 #include "AliAmptRndm.h"
35 #include "AliDecayer.h"
39 AliGenAmpt::AliGenAmpt()
66 fPhiMaxJet(TMath::TwoPi()),
76 fNoHeavyQuarks(kFALSE),
85 fHeader(new AliGenAmptEventHeader("Ampt")),
91 AliAmptRndm::SetAmptRandom(GetRandom());
94 AliGenAmpt::AliGenAmpt(Int_t npart)
121 fPhiMaxJet(2. * TMath::Pi()),
131 fNoHeavyQuarks(kFALSE),
140 fHeader(new AliGenAmptEventHeader("Ampt")),
144 // Default PbPb collisions at 2.76 TeV
148 fTitle= "Particle Generator using AMPT";
149 AliAmptRndm::SetAmptRandom(GetRandom());
152 AliGenAmpt::~AliGenAmpt()
155 if ( fDsigmaDb) delete fDsigmaDb;
156 if ( fDnDb) delete fDnDb;
157 if ( fHeader) delete fHeader;
160 void AliGenAmpt::Init()
166 fProjectile.Resize(8);
168 fAmpt = new TAmpt(fEnergyCMS, fFrame, fProjectile, fTarget,
169 fAProjectile, fZProjectile, fATarget, fZTarget,
170 fMinImpactParam, fMaxImpactParam);
173 fAmpt->SetIHPR2(2, fRadiation);
174 fAmpt->SetIHPR2(3, fTrigger);
175 fAmpt->SetIHPR2(6, fShadowing);
176 fAmpt->SetIHPR2(12, fDecaysOff);
177 fAmpt->SetIHPR2(21, fKeep);
178 fAmpt->SetHIPR1(8, fPtHardMin);
179 fAmpt->SetHIPR1(9, fPtHardMax);
180 fAmpt->SetHIPR1(10, fPtMinJet);
181 fAmpt->SetHIPR1(50, fSimpleJet);
184 // fQuench = 0: no quenching
185 // fQuench = 1: Hijing default
186 // fQuench = 2: new LHC parameters for HIPR1(11) and HIPR1(14)
187 // fQuench = 3: new RHIC parameters for HIPR1(11) and HIPR1(14)
188 // fQuench = 4: new LHC parameters with log(e) dependence
189 // fQuench = 5: new RHIC parameters with log(e) dependence
190 fAmpt->SetIHPR2(50, 0);
192 fAmpt->SetIHPR2(4, 1);
194 fAmpt->SetIHPR2(4, 0);
197 fAmpt->SetHIPR1(14, 1.1);
198 fAmpt->SetHIPR1(11, 3.7);
199 } else if (fQuench == 3) {
200 fAmpt->SetHIPR1(14, 0.20);
201 fAmpt->SetHIPR1(11, 2.5);
202 } else if (fQuench == 4) {
203 fAmpt->SetIHPR2(50, 1);
204 fAmpt->SetHIPR1(14, 4.*0.34);
205 fAmpt->SetHIPR1(11, 3.7);
206 } else if (fQuench == 5) {
207 fAmpt->SetIHPR2(50, 1);
208 fAmpt->SetHIPR1(14, 0.34);
209 fAmpt->SetHIPR1(11, 2.5);
213 if (fNoHeavyQuarks) {
214 fAmpt->SetIHPR2(49, 1);
216 fAmpt->SetIHPR2(49, 0);
220 fAmpt->SetIsoft(fIsoft);
221 fAmpt->SetNtMax(fNtMax);
222 fAmpt->SetIpop(fIpop);
224 fAmpt->SetAlpha(fAlpha);
225 fAmpt->SetStringFrag(fStringA, fStringB);
232 EvaluateCrossSections();
234 fAmpt->SetReactionPlaneAngle(0.0);
238 void AliGenAmpt::Generate()
240 // Generate one event
242 Float_t polar[3] = {0,0,0};
243 Float_t origin[3] = {0,0,0};
244 Float_t origin0[3] = {0,0,0};
249 // converts from mm/c to s
250 const Float_t kconv = 0.001/2.99792458e8;
254 Int_t j, kf, ks, ksp, imo;
258 for (j = 0;j < 3; j++)
259 origin0[j] = fOrigin[j];
260 //time0 = fTimeOrigin;
262 if(fVertexSmear == kPerEvent) {
264 for (j=0; j < 3; j++)
265 origin0[j] = fVertex[j];
269 Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
273 // Generate random reaction plane angle if requested
275 TRandom *r=AliAmptRndm::GetAmptRandom();
276 fAmpt->SetReactionPlaneAngle(TMath::TwoPi()*r->Rndm());
279 // Generate one event
280 Int_t fpemask = gSystem->GetFPEMask();
281 gSystem->SetFPEMask(0);
282 fAmpt->GenerateEvent();
283 gSystem->SetFPEMask(fpemask);
288 fAmpt->ImportParticles(&fParticles,"All");
289 Int_t np = fParticles.GetEntriesFast();
293 if (fTrigger != kNoTrigger) {
298 AliDecayer *decayer = 0;
300 // decayer = gMC->GetDecayer();
301 decayer = fDecayer; //AMPT does not do the strong decays per dafault
303 if (decayer&&fDecay) {
304 TClonesArray arr("TParticle",100);
305 for( Int_t nLoop=0; nLoop!=2; ++nLoop) { // In order to produce more than one generation of decays: NumberOfNestedLoops set to 2
307 for (Int_t i = 0; i < np; i++) {
308 TParticle *iparticle = (TParticle *)fParticles.At(i);
309 if (!Stable(iparticle)) // true if particle has daughters already
311 kf = TMath::Abs(iparticle->GetPdgCode());
314 if( !IsThisAKnownParticle(iparticle) ) continue; // skip undesired particles
316 if (0) { // this turned out to be too cumbersome!
317 if (kf!=331&&kf!=3114&&kf!=3114&&kf!=411&&kf!=-4122&&kf!=-3324&&kf!=-3312&&kf!=-3114&&
318 kf!=-311&&kf!=3214&&kf!=-3214&&kf!=-433&&kf!=413&&kf!=3122&&kf!=-3122&&kf!=-413&&
319 kf!=-421&&kf!=-423&&kf!=3324&&kf!=-313&&kf!=213&&kf!=-213&&kf!=3314&&kf!=3222&&
320 kf!=-3222&&kf!=3224&&kf!=-3224&&kf!=-4212&&kf!=4212&&kf!=433&&kf!=423&&kf!=-3322&&
322 continue; //decay eta',Sigma*+,Sigma*-,D+,Lambda_c-,Xi*0_bar,Xi-_bar,Sigma*-,
323 // K0_bar,Sigma*0,Sigma*0_bar,D*_s-,D*+,Lambda0,Lambda0_bar,D*-
324 // D0_bar,D*0_bar,Xi*0,K*0_bar,rho+,rho-,Xi*-,Sigma-,
325 // Sigma+,Sigma*+,Sigma*-,Sigma_c-,Sigma_c+,D*_s+,D*0,Xi0_bar
327 //} else { // really only decay particles if there are not known to Geant3
328 // if (gMC->IdFromPDG(kf)>0)
331 if (0) { // defining the particle for Geant3 leads to a floating point exception.
332 TParticlePDG *pdg = iparticle->GetPDG(1);
333 //pdg->Print(); printf("%s\n",pdg->ParticleClass());
334 TString ptype(pdg->ParticleClass());
335 TMCParticleType mctype(kPTUndefined);
336 if (ptype=="Baryon" || ptype=="Meson")
338 gMC->DefineParticle(pdg->PdgCode(), pdg->GetName(), mctype, pdg->Mass(), pdg->Charge(), pdg->Lifetime(),
339 ptype,pdg->Width(), (Int_t)pdg->Spin(), (Int_t)pdg->Parity(), 0,
340 (Int_t)pdg->Isospin(), 0, 0, 0, 0, pdg->Stable());
341 gMC->SetUserDecay(pdg->PdgCode());
345 TLorentzVector pmom(iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy());
346 decayer->Decay(kf,&pmom);
347 decayer->ImportParticles(&arr);
348 Int_t ndecayed = arr.GetEntries();
350 if (np2+ndecayed>fParticles.GetSize())
351 fParticles.Expand(2*fParticles.GetSize());
353 // iparticle->SetStatusCode(2); to be compatible with Hijing
354 iparticle->SetFirstDaughter(np2);
355 for (Int_t jj = 1; jj < ndecayed; jj++) {
356 TParticle *jp = (TParticle *)arr.At(jj);
357 if (jp->GetFirstMother()!=1)
359 TParticle *newp = new(fParticles[np2]) TParticle(jp->GetPdgCode(),
360 0, //1, //to be compatible with Hijing
365 jp->Px(),jp->Py(),jp->Pz(),jp->Energy(),
366 jp->Vx(),jp->Vy(),jp->Vz(),jp->T());
367 newp->SetUniqueID( jp->GetStatusCode() );
369 } // end of jj->nDecayedParticles
370 iparticle->SetLastDaughter(np2-1);
371 } // end of nDecayedPrticles>1
373 np = fParticles.GetEntries();
375 AliError(Form("Something is fishy: %d %d\n", np,np2));
377 } // end of nLoop->NumberOfNestedLoops
380 AliError("No decayer found, but fDecay==kTRUE!");
387 Int_t* newPos = new Int_t[np];
388 Int_t* pSelected = new Int_t[np];
390 for (Int_t i = 0; i < np; i++) {
396 //TParticle * iparticle = (TParticle *) fParticles.At(0);
397 fVertex[0] = origin0[0];
398 fVertex[1] = origin0[1];
399 fVertex[2] = origin0[2];
402 // First select parent particles
403 for (Int_t i = 0; i < np; i++) {
404 TParticle *iparticle = (TParticle *) fParticles.At(i);
406 // Is this a parent particle ?
407 if (Stable(iparticle)) continue; // quit if particle has no daughters
408 Bool_t selected = kTRUE;
409 Bool_t hasSelectedDaughters = kFALSE;
410 kf = iparticle->GetPdgCode();
411 ks = iparticle->GetStatusCode();
416 selected = KinematicSelection(iparticle, 0) && SelectFlavor(kf);
417 hasSelectedDaughters = DaughtersSelection(iparticle);
419 // Put particle on the stack if it is either selected or
420 // it is the mother of at least one seleted particle
421 if (selected || hasSelectedDaughters) {
425 } // particle loop parents
427 // Now select the final state particles
428 fProjectileSpecn = 0;
429 fProjectileSpecp = 0;
432 for (Int_t i = 0; i<np; i++) {
433 TParticle *iparticle = (TParticle *) fParticles.At(i);
434 // Is this a final state particle ?
435 if (!Stable(iparticle)) continue; // quit if particle has daughters
436 Bool_t selected = kTRUE;
437 kf = iparticle->GetPdgCode();
440 ks = iparticle->GetStatusCode();
441 ksp = iparticle->GetUniqueID();
443 // --------------------------------------------------------------------------
444 // Count spectator neutrons and protons
445 if(ksp == 0 || ksp == 1) {
446 if(kf == kNeutron) fProjectileSpecn += 1;
447 if(kf == kProton) fProjectileSpecp += 1;
448 } else if(ksp == 10 || ksp == 11) {
449 if(kf == kNeutron) fTargetSpecn += 1;
450 if(kf == kProton) fTargetSpecp += 1;
452 // --------------------------------------------------------------------------
454 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
455 if (!fSpectators && selected)
456 selected = (ksp != 0 && ksp != 1 && ksp != 10 && ksp != 11);
459 // Put particle on the stack if selected
463 if (0) printf("---> %d %d %d %s\n",i,nc,kf,iparticle->GetName());
465 } // particle loop final state
467 // Write particles to stack
468 for (Int_t i = 0; i<np; i++) {
470 TParticle *iparticle = (TParticle *) fParticles.At(i);
471 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
472 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
473 kf = iparticle->GetPdgCode();
474 ks = iparticle->GetStatusCode();
475 p[0] = iparticle->Px();
476 p[1] = iparticle->Py();
477 p[2] = iparticle->Pz() * sign;
478 origin[0] = origin0[0]+iparticle->Vx()/10;
479 origin[1] = origin0[1]+iparticle->Vy()/10;
480 origin[2] = origin0[2]+iparticle->Vz()/10;
481 tof = time0+kconv * iparticle->T();
484 TParticle* mother = 0;
485 TMCProcess procID = (TMCProcess) iparticle->GetUniqueID();
487 imo = iparticle->GetFirstMother();
488 mother = (TParticle *) fParticles.At(imo);
489 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
490 } else { // if has no mothers then it was created by AMPT
492 procID = kPPrimary; // reseting to ALIROOT convention
494 procID = kPNoProcess; // for expectators
496 Bool_t tFlag = (fTrackIt && !hasDaughter);
497 PushTrack(tFlag,imo,kf,p,origin,polar,tof,procID,nt, 1., ks);
506 AliInfo(Form("\n I've put %i particles on the stack \n",nc));
509 if (jev >= fNpart || fNpart == -1) {
510 fKineBias = Float_t(fNpart)/Float_t(fTrials);
511 AliInfo(Form("\n Trials: %i %i %i\n",fTrials, fNpart, jev));
517 SetHighWaterMark(nt);
520 Bool_t AliGenAmpt::IsThisAKnownParticle(TParticle *thisGuy)
522 // In order to prevent AMPT to introduce weird particles into the decayer and transporter
523 // blame cperez@cern.ch for this method
525 Int_t pdgcode = TMath::Abs( thisGuy->GetPdgCode() );
527 Int_t myFavoriteParticles[ 38] = { 3322, 3314, 3312, 3224, 3222, // Xi0 Xi*+- Xi+- Sigma*-+ Sigma-+
528 3214, 3212, 3122, 3114, 3112, // Sigma*0 Sigma0 Lambda0 Sigma*+- Sigma+-
529 2224, 2214, 2212, 2114, 2112, // Delta--++ Delta-+ proton Delta0 neutron
530 1114, 323, 321, 313, 311, // Delta+- K*-+ K-+ K*0 K0
531 213, 211, 11, 22, 111, // rho-+ pi-+ e+- gamma pi0
532 113, 130, 221, 223, 310, // rho0 K_L0 eta omega K_S0
533 331, 333, 3324, 431, 421, // eta' phi Xi*0 Ds-+ D0
534 411, 413, 13 // D-+ D*-+ mu+-
537 Bool_t found = kFALSE;
538 for(Int_t i=0; i!=38; ++i)
539 if( myFavoriteParticles[i] == pdgcode ) {
547 void AliGenAmpt::EvaluateCrossSections()
549 // Glauber Calculation of geometrical x-section
551 Float_t xTot = 0.; // barn
552 Float_t xTotHard = 0.; // barn
553 Float_t xPart = 0.; // barn
554 Float_t xPartHard = 0.; // barn
555 Float_t sigmaHard = 0.1; // mbarn
557 Float_t bMax = fAmpt->GetHIPR1(34)+fAmpt->GetHIPR1(35);
558 const Float_t kdib = 0.2;
559 Int_t kMax = Int_t((bMax-bMin)/kdib)+1;
561 printf("\n Projectile Radius (fm): %f \n",fAmpt->GetHIPR1(34));
562 printf("\n Target Radius (fm): %f \n",fAmpt->GetHIPR1(35));
565 Float_t oldvalue= 0.;
566 Float_t* b = new Float_t[kMax]; memset(b,0,kMax*sizeof(Float_t));
567 Float_t* si1 = new Float_t[kMax]; memset(si1,0,kMax*sizeof(Float_t));
568 Float_t* si2 = new Float_t[kMax]; memset(si2,0,kMax*sizeof(Float_t));
569 for (i = 0; i < kMax; i++) {
570 Float_t xb = bMin+i*kdib;
571 Float_t ov=fAmpt->Profile(xb);
572 Float_t gb = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fAmpt->GetHINT1(12)*ov));
573 Float_t gbh = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
576 printf("profile %f %f %f\n", xb, ov, fAmpt->GetHINT1(12));
578 if (xb > fMinImpactParam && xb < fMaxImpactParam) {
583 if ((oldvalue) && ((xTot-oldvalue)/oldvalue<0.0001))
586 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
587 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
595 printf("\n Total cross section (barn): %f \n",xTot);
596 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
597 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
598 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
600 // Store result as a graph
605 fDsigmaDb = new TGraph(i, b, si1);
607 fDnDb = new TGraph(i, b, si2);
610 Bool_t AliGenAmpt::DaughtersSelection(TParticle* iparticle)
612 // Looks recursively if one of the daughters has been selected
613 //printf("\n Consider daughters %d:",iparticle->GetPdgCode());
616 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
617 Bool_t selected = kFALSE;
619 imin = iparticle->GetFirstDaughter();
620 imax = iparticle->GetLastDaughter();
621 for (Int_t i = imin; i <= imax; i++){
622 TParticle * jparticle = (TParticle *) fParticles.At(i);
623 Int_t ip = jparticle->GetPdgCode();
624 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
625 selected=kTRUE; break;
627 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
635 Bool_t AliGenAmpt::SelectFlavor(Int_t pid)
637 // Select flavor of particle
639 // 4: charm and beauty
646 Int_t ifl = TMath::Abs(pid/100);
647 if (ifl > 10) ifl/=10;
648 res = (fFlavor == ifl);
651 // This part if gamma writing is inhibited
653 res = res && (pid != kGamma && pid != kPi0);
658 Bool_t AliGenAmpt::Stable(TParticle* particle) const
660 // Return true for a stable particle
664 if (particle->GetFirstDaughter() < 0 )
672 void AliGenAmpt::MakeHeader()
674 // Fills the event header, to be called after each event
676 fHeader->SetNProduced(fNprimaries);
677 fHeader->SetImpactParameter(fAmpt->GetHINT1(19));
678 fHeader->SetTotalEnergy(fAmpt->GetEATT());
679 fHeader->SetHardScatters(fAmpt->GetJATT());
680 fHeader->SetParticipants(fAmpt->GetNP(), fAmpt->GetNT());
681 fHeader->SetCollisions(fAmpt->GetN0(),
685 fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
686 fTargetSpecn,fTargetSpecp);
687 //fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
688 fHeader->SetReactionPlaneAngle(fAmpt->GetReactionPlaneAngle());
689 //printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
691 // 4-momentum vectors of the triggered jets.
692 // Before final state gluon radiation.
693 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(21),
696 fAmpt->GetHINT1(24));
698 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(31),
701 fAmpt->GetHINT1(34));
702 // After final state gluon radiation.
703 TLorentzVector* jet3 = new TLorentzVector(fAmpt->GetHINT1(26),
706 fAmpt->GetHINT1(29));
708 TLorentzVector* jet4 = new TLorentzVector(fAmpt->GetHINT1(36),
711 fAmpt->GetHINT1(39));
712 fHeader->SetJets(jet1, jet2, jet3, jet4);
713 // Bookkeeping for kinematic bias
714 fHeader->SetTrials(fTrials);
716 fHeader->SetPrimaryVertex(fVertex);
717 fHeader->SetInteractionTime(fEventTime);
719 fCollisionGeometry = fHeader;
724 Bool_t AliGenAmpt::CheckTrigger()
726 // Check the kinematic trigger condition
728 Bool_t triggered = kFALSE;
732 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(26),
735 fAmpt->GetHINT1(29));
737 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(36),
740 fAmpt->GetHINT1(39));
741 Double_t eta1 = jet1->Eta();
742 Double_t eta2 = jet2->Eta();
743 Double_t phi1 = jet1->Phi();
744 Double_t phi2 = jet2->Phi();
745 //printf("\n Trigger: %f %f %f %f", fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
746 if ( (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&
747 phi1 < fPhiMaxJet && phi1 > fPhiMinJet)
749 (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&
750 phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
753 } else if (fTrigger == 2) {
755 Int_t np = fParticles.GetEntriesFast();
756 for (Int_t i = 0; i < np; i++) {
757 TParticle* part = (TParticle*) fParticles.At(i);
758 Int_t kf = part->GetPdgCode();
759 Int_t ksp = part->GetUniqueID();
760 if (kf == 22 && ksp == 40) {
761 Float_t phi = part->Phi();
762 Float_t eta = part->Eta();
763 if (eta < fEtaMaxJet &&
769 } // check phi,eta within limits