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>
29 #include "AliGenHijingEventHeader.h"
30 #define AliGenAmptEventHeader AliGenHijingEventHeader
31 #include "AliAmptRndm.h"
37 AliGenAmpt::AliGenAmpt()
63 fPhiMaxJet(TMath::TwoPi()),
73 fNoHeavyQuarks(kFALSE),
82 fHeader(new AliGenAmptEventHeader("Ampt"))
86 AliAmptRndm::SetAmptRandom(GetRandom());
89 AliGenAmpt::AliGenAmpt(Int_t npart)
115 fPhiMaxJet(2. * TMath::Pi()),
125 fNoHeavyQuarks(kFALSE),
134 fHeader(new AliGenAmptEventHeader("Ampt"))
136 // Default PbPb collisions at 5.5 TeV
140 fTitle= "Particle Generator using AMPT";
141 AliAmptRndm::SetAmptRandom(GetRandom());
144 AliGenAmpt::~AliGenAmpt()
147 if ( fDsigmaDb) delete fDsigmaDb;
148 if ( fDnDb) delete fDnDb;
149 if ( fHeader) delete fHeader;
152 void AliGenAmpt::Init()
158 fProjectile.Resize(8);
160 fAmpt = new TAmpt(fEnergyCMS, fFrame, fProjectile, fTarget,
161 fAProjectile, fZProjectile, fATarget, fZTarget,
162 fMinImpactParam, fMaxImpactParam);
165 fAmpt->SetIHPR2(2, fRadiation);
166 fAmpt->SetIHPR2(3, fTrigger);
167 fAmpt->SetIHPR2(6, fShadowing);
168 fAmpt->SetIHPR2(12, fDecaysOff);
169 fAmpt->SetIHPR2(21, fKeep);
170 fAmpt->SetHIPR1(8, fPtHardMin);
171 fAmpt->SetHIPR1(9, fPtHardMax);
172 fAmpt->SetHIPR1(10, fPtMinJet);
173 fAmpt->SetHIPR1(50, fSimpleJet);
176 // fQuench = 0: no quenching
177 // fQuench = 1: Hijing default
178 // fQuench = 2: new LHC parameters for HIPR1(11) and HIPR1(14)
179 // fQuench = 3: new RHIC parameters for HIPR1(11) and HIPR1(14)
180 // fQuench = 4: new LHC parameters with log(e) dependence
181 // fQuench = 5: new RHIC parameters with log(e) dependence
182 fAmpt->SetIHPR2(50, 0);
184 fAmpt->SetIHPR2(4, 1);
186 fAmpt->SetIHPR2(4, 0);
189 fAmpt->SetHIPR1(14, 1.1);
190 fAmpt->SetHIPR1(11, 3.7);
191 } else if (fQuench == 3) {
192 fAmpt->SetHIPR1(14, 0.20);
193 fAmpt->SetHIPR1(11, 2.5);
194 } else if (fQuench == 4) {
195 fAmpt->SetIHPR2(50, 1);
196 fAmpt->SetHIPR1(14, 4.*0.34);
197 fAmpt->SetHIPR1(11, 3.7);
198 } else if (fQuench == 5) {
199 fAmpt->SetIHPR2(50, 1);
200 fAmpt->SetHIPR1(14, 0.34);
201 fAmpt->SetHIPR1(11, 2.5);
205 if (fNoHeavyQuarks) {
206 fAmpt->SetIHPR2(49, 1);
208 fAmpt->SetIHPR2(49, 0);
212 fAmpt->SetIsoft(fIsoft);
213 fAmpt->SetNtMax(fNtMax);
214 fAmpt->SetIpop(fIpop);
216 fAmpt->SetAlpha(fAlpha);
217 fAmpt->SetStringFrag(fStringA, fStringB);
224 EvaluateCrossSections();
227 void AliGenAmpt::Generate()
229 // Generate one event
231 Float_t polar[3] = {0,0,0};
232 Float_t origin[3] = {0,0,0};
233 Float_t origin0[3] = {0,0,0};
237 // converts from mm/c to s
238 const Float_t kconv = 0.001/2.99792458e8;
242 Int_t j, kf, ks, ksp, imo;
246 for (j = 0;j < 3; j++)
247 origin0[j] = fOrigin[j];
249 if(fVertexSmear == kPerEvent) {
251 for (j=0; j < 3; j++)
252 origin0[j] = fVertex[j];
255 Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
258 // Generate one event
259 Int_t fpemask = gSystem->GetFPEMask();
260 gSystem->SetFPEMask(0);
261 fAmpt->GenerateEvent();
262 gSystem->SetFPEMask(fpemask);
265 fAmpt->ImportParticles(&fParticles,"All");
266 Int_t np = fParticles.GetEntriesFast();
270 if (fTrigger != kNoTrigger) {
278 Int_t* newPos = new Int_t[np];
279 Int_t* pSelected = new Int_t[np];
281 for (Int_t i = 0; i < np; i++) {
287 //TParticle * iparticle = (TParticle *) fParticles.At(0);
288 fVertex[0] = origin0[0];
289 fVertex[1] = origin0[1];
290 fVertex[2] = origin0[2];
292 // First select parent particles
293 for (Int_t i = 0; i < np; i++) {
294 TParticle *iparticle = (TParticle *) fParticles.At(i);
296 // Is this a parent particle ?
297 if (Stable(iparticle)) continue;
298 Bool_t selected = kTRUE;
299 Bool_t hasSelectedDaughters = kFALSE;
300 kf = iparticle->GetPdgCode();
301 ks = iparticle->GetStatusCode();
306 selected = KinematicSelection(iparticle, 0) && SelectFlavor(kf);
307 hasSelectedDaughters = DaughtersSelection(iparticle);
309 // Put particle on the stack if it is either selected or
310 // it is the mother of at least one seleted particle
311 if (selected || hasSelectedDaughters) {
315 } // particle loop parents
317 // Now select the final state particles
318 fProjectileSpecn = 0;
319 fProjectileSpecp = 0;
322 for (Int_t i = 0; i<np; i++) {
323 TParticle *iparticle = (TParticle *) fParticles.At(i);
324 // Is this a final state particle ?
325 if (!Stable(iparticle))
327 Bool_t selected = kTRUE;
328 kf = iparticle->GetPdgCode();
329 ks = iparticle->GetStatusCode();
330 ksp = iparticle->GetUniqueID();
332 // --------------------------------------------------------------------------
333 // Count spectator neutrons and protons
334 if(ksp == 0 || ksp == 1) {
335 if(kf == kNeutron) fProjectileSpecn += 1;
336 if(kf == kProton) fProjectileSpecp += 1;
337 } else if(ksp == 10 || ksp == 11) {
338 if(kf == kNeutron) fTargetSpecn += 1;
339 if(kf == kProton) fTargetSpecp += 1;
341 // --------------------------------------------------------------------------
343 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
344 if (!fSpectators && selected)
345 selected = (ksp != 0 && ksp != 1 && ksp != 10 && ksp != 11);
348 // Put particle on the stack if selected
353 } // particle loop final state
355 //Time of the interactions
357 if (fPileUpTimeWindow > 0.)
358 tInt = fPileUpTimeWindow * (2. * gRandom->Rndm() - 1.);
360 // Write particles to stack
361 for (Int_t i = 0; i<np; i++) {
362 TParticle *iparticle = (TParticle *) fParticles.At(i);
363 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
364 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
366 kf = iparticle->GetPdgCode();
367 ks = iparticle->GetStatusCode();
368 p[0] = iparticle->Px();
369 p[1] = iparticle->Py();
370 p[2] = iparticle->Pz() * sign;
371 origin[0] = origin0[0]+iparticle->Vx()/10;
372 origin[1] = origin0[1]+iparticle->Vy()/10;
373 origin[2] = origin0[2]+iparticle->Vz()/10;
376 if (TestBit(kVertexRange)) {
377 fEventTime = sign * origin0[2] / 2.99792458e10;
378 tof = kconv * iparticle->T() + fEventTime;
380 tof = kconv * iparticle->T();
382 if (fPileUpTimeWindow > 0.) tof += tInt;
385 TParticle* mother = 0;
387 imo = iparticle->GetFirstMother();
388 mother = (TParticle *) fParticles.At(imo);
389 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
391 Bool_t tFlag = (fTrackIt && !hasDaughter);
392 PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
401 AliInfo(Form("\n I've put %i particles on the stack \n",nc));
404 if (jev >= fNpart || fNpart == -1) {
405 fKineBias = Float_t(fNpart)/Float_t(fTrials);
406 AliInfo(Form("\n Trials: %i %i %i\n",fTrials, fNpart, jev));
412 SetHighWaterMark(nt);
415 void AliGenAmpt::EvaluateCrossSections()
417 // Glauber Calculation of geometrical x-section
419 Float_t xTot = 0.; // barn
420 Float_t xTotHard = 0.; // barn
421 Float_t xPart = 0.; // barn
422 Float_t xPartHard = 0.; // barn
423 Float_t sigmaHard = 0.1; // mbarn
425 Float_t bMax = fAmpt->GetHIPR1(34)+fAmpt->GetHIPR1(35);
426 const Float_t kdib = 0.2;
427 Int_t kMax = Int_t((bMax-bMin)/kdib)+1;
429 printf("\n Projectile Radius (fm): %f \n",fAmpt->GetHIPR1(34));
430 printf("\n Target Radius (fm): %f \n",fAmpt->GetHIPR1(35));
433 Float_t oldvalue= 0.;
434 Float_t* b = new Float_t[kMax]; memset(b,0,kMax*sizeof(Float_t));
435 Float_t* si1 = new Float_t[kMax]; memset(si1,0,kMax*sizeof(Float_t));
436 Float_t* si2 = new Float_t[kMax]; memset(si2,0,kMax*sizeof(Float_t));
437 for (i = 0; i < kMax; i++) {
438 Float_t xb = bMin+i*kdib;
439 Float_t ov=fAmpt->Profile(xb);
440 Float_t gb = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fAmpt->GetHINT1(12)*ov));
441 Float_t gbh = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
444 printf("profile %f %f %f\n", xb, ov, fAmpt->GetHINT1(12));
446 if (xb > fMinImpactParam && xb < fMaxImpactParam) {
451 if ((oldvalue) && ((xTot-oldvalue)/oldvalue<0.0001))
454 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
455 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
463 printf("\n Total cross section (barn): %f \n",xTot);
464 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
465 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
466 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
468 // Store result as a graph
473 fDsigmaDb = new TGraph(i, b, si1);
475 fDnDb = new TGraph(i, b, si2);
478 Bool_t AliGenAmpt::DaughtersSelection(TParticle* iparticle)
480 // Looks recursively if one of the daughters has been selected
481 //printf("\n Consider daughters %d:",iparticle->GetPdgCode());
484 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
485 Bool_t selected = kFALSE;
487 imin = iparticle->GetFirstDaughter();
488 imax = iparticle->GetLastDaughter();
489 for (Int_t i = imin; i <= imax; i++){
490 TParticle * jparticle = (TParticle *) fParticles.At(i);
491 Int_t ip = jparticle->GetPdgCode();
492 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
493 selected=kTRUE; break;
495 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
503 Bool_t AliGenAmpt::SelectFlavor(Int_t pid)
505 // Select flavor of particle
507 // 4: charm and beauty
514 Int_t ifl = TMath::Abs(pid/100);
515 if (ifl > 10) ifl/=10;
516 res = (fFlavor == ifl);
519 // This part if gamma writing is inhibited
521 res = res && (pid != kGamma && pid != kPi0);
526 Bool_t AliGenAmpt::Stable(TParticle* particle) const
528 // Return true for a stable particle
530 if (particle->GetFirstDaughter() < 0 )
538 void AliGenAmpt::MakeHeader()
540 // Fills the event header, to be called after each event
542 fHeader->SetNProduced(fNprimaries);
543 fHeader->SetImpactParameter(fAmpt->GetHINT1(19));
544 fHeader->SetTotalEnergy(fAmpt->GetEATT());
545 fHeader->SetHardScatters(fAmpt->GetJATT());
546 fHeader->SetParticipants(fAmpt->GetNP(), fAmpt->GetNT());
547 fHeader->SetCollisions(fAmpt->GetN0(),
551 fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
552 fTargetSpecn,fTargetSpecp);
553 fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
554 //printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
556 // 4-momentum vectors of the triggered jets.
557 // Before final state gluon radiation.
558 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(21),
561 fAmpt->GetHINT1(24));
563 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(31),
566 fAmpt->GetHINT1(34));
567 // After final state gluon radiation.
568 TLorentzVector* jet3 = new TLorentzVector(fAmpt->GetHINT1(26),
571 fAmpt->GetHINT1(29));
573 TLorentzVector* jet4 = new TLorentzVector(fAmpt->GetHINT1(36),
576 fAmpt->GetHINT1(39));
577 fHeader->SetJets(jet1, jet2, jet3, jet4);
578 // Bookkeeping for kinematic bias
579 fHeader->SetTrials(fTrials);
581 fHeader->SetPrimaryVertex(fVertex);
582 fHeader->SetInteractionTime(fEventTime);
584 fCollisionGeometry = fHeader;
589 Bool_t AliGenAmpt::CheckTrigger()
591 // Check the kinematic trigger condition
593 Bool_t triggered = kFALSE;
597 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(26),
600 fAmpt->GetHINT1(29));
602 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(36),
605 fAmpt->GetHINT1(39));
606 Double_t eta1 = jet1->Eta();
607 Double_t eta2 = jet2->Eta();
608 Double_t phi1 = jet1->Phi();
609 Double_t phi2 = jet2->Phi();
610 //printf("\n Trigger: %f %f %f %f", fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
611 if ( (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&
612 phi1 < fPhiMaxJet && phi1 > fPhiMinJet)
614 (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&
615 phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
618 } else if (fTrigger == 2) {
620 Int_t np = fParticles.GetEntriesFast();
621 for (Int_t i = 0; i < np; i++) {
622 TParticle* part = (TParticle*) fParticles.At(i);
623 Int_t kf = part->GetPdgCode();
624 Int_t ksp = part->GetUniqueID();
625 if (kf == 22 && ksp == 40) {
626 Float_t phi = part->Phi();
627 Float_t eta = part->Eta();
628 if (eta < fEtaMaxJet &&
634 } // check phi,eta within limits