const
[u/mrichter/AliRoot.git] / TAmpt / AliGenAmpt.cxx
CommitLineData
0119ef9a 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/* $Id$ */
17
18// Generator using AMPT as an external generator
19
20#include "AliGenAmpt.h"
21
22#include <TClonesArray.h>
23#include <TGraph.h>
24#include <TAmpt.h>
25#include <TLorentzVector.h>
26#include <TPDGCode.h>
27#include <TParticle.h>
28
29#include "AliGenHijingEventHeader.h"
30#define AliGenAmptEventHeader AliGenHijingEventHeader
31#include "AliAmptRndm.h"
32#include "AliLog.h"
33#include "AliRun.h"
34
35ClassImp(AliGenAmpt)
36
37AliGenAmpt::AliGenAmpt()
38 : AliGenMC(),
39 fFrame("CMS"),
40 fMinImpactParam(0.),
41 fMaxImpactParam(5.),
42 fKeep(0),
43 fQuench(0),
44 fShadowing(1),
45 fDecaysOff(1),
46 fTrigger(0),
47 fEvaluate(0),
48 fSelectAll(0),
49 fFlavor(0),
50 fKineBias(0.),
51 fTrials(0),
52 fXsection(0.),
53 fAmpt(0),
54 fPtHardMin(2.0),
55 fPtHardMax(-1),
56 fSpectators(1),
57 fDsigmaDb(0),
58 fDnDb(0),
59 fPtMinJet(-2.5),
60 fEtaMinJet(-20.),
61 fEtaMaxJet(+20.),
62 fPhiMinJet(0.),
63 fPhiMaxJet(TMath::TwoPi()),
64 fRadiation(3),
65 fSimpleJet(kFALSE),
66 fNoGammas(kFALSE),
67 fProjectileSpecn(0),
68 fProjectileSpecp(0),
69 fTargetSpecn(0),
70 fTargetSpecp(0),
71 fLHC(kFALSE),
72 fRandomPz(kFALSE),
73 fNoHeavyQuarks(kFALSE),
74 fEventTime(0.),
75 fIsoft(1),
76 fNtMax(150),
77 fIpop(1),
78 fXmu(3.2264),
a004b331 79 fAlpha(1./3),
80 fStringA(0.5),
81 fStringB(0.9),
0119ef9a 82 fHeader(new AliGenAmptEventHeader("Ampt"))
83{
84 // Constructor
85 fEnergyCMS = 5500.;
86 AliAmptRndm::SetAmptRandom(GetRandom());
87}
88
89AliGenAmpt::AliGenAmpt(Int_t npart)
90 : AliGenMC(npart),
91 fFrame("CMS"),
92 fMinImpactParam(0.),
93 fMaxImpactParam(5.),
94 fKeep(0),
95 fQuench(0),
96 fShadowing(1),
97 fDecaysOff(1),
98 fTrigger(0),
99 fEvaluate(0),
100 fSelectAll(0),
101 fFlavor(0),
102 fKineBias(0.),
103 fTrials(0),
104 fXsection(0.),
105 fAmpt(0),
106 fPtHardMin(2.0),
107 fPtHardMax(-1),
108 fSpectators(1),
109 fDsigmaDb(0),
110 fDnDb(0),
111 fPtMinJet(-2.5),
112 fEtaMinJet(-20.),
113 fEtaMaxJet(+20.),
114 fPhiMinJet(0.),
115 fPhiMaxJet(2. * TMath::Pi()),
116 fRadiation(3),
117 fSimpleJet(kFALSE),
118 fNoGammas(kFALSE),
119 fProjectileSpecn(0),
120 fProjectileSpecp(0),
121 fTargetSpecn(0),
122 fTargetSpecp(0),
123 fLHC(kFALSE),
124 fRandomPz(kFALSE),
125 fNoHeavyQuarks(kFALSE),
126 fEventTime(0.),
127 fIsoft(1),
128 fNtMax(150),
129 fIpop(1),
130 fXmu(3.2264),
a004b331 131 fAlpha(1./3),
132 fStringA(0.5),
133 fStringB(0.9),
0119ef9a 134 fHeader(new AliGenAmptEventHeader("Ampt"))
135{
136 // Default PbPb collisions at 5.5 TeV
137
138 fEnergyCMS = 5500.;
139 fName = "Ampt";
140 fTitle= "Particle Generator using AMPT";
141 AliAmptRndm::SetAmptRandom(GetRandom());
142}
143
144AliGenAmpt::~AliGenAmpt()
145{
146 // Destructor
147 if ( fDsigmaDb) delete fDsigmaDb;
148 if ( fDnDb) delete fDnDb;
149 if ( fHeader) delete fHeader;
150}
151
152void AliGenAmpt::Init()
153{
154 // Initialisation
155
156 fFrame.Resize(8);
157 fTarget.Resize(8);
158 fProjectile.Resize(8);
159
160 fAmpt = new TAmpt(fEnergyCMS, fFrame, fProjectile, fTarget,
161 fAProjectile, fZProjectile, fATarget, fZTarget,
162 fMinImpactParam, fMaxImpactParam);
163 SetMC(fAmpt);
164
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);
174
175 // Quenching
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);
183 if (fQuench > 0)
184 fAmpt->SetIHPR2(4, 1);
185 else
186 fAmpt->SetIHPR2(4, 0);
187
188 if (fQuench == 2) {
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);
202 }
203
204 // Heavy quarks
205 if (fNoHeavyQuarks) {
206 fAmpt->SetIHPR2(49, 1);
207 } else {
208 fAmpt->SetIHPR2(49, 0);
209 }
210
211 // Ampt specific
212 fAmpt->SetIsoft(fIsoft);
213 fAmpt->SetNtMax(fNtMax);
214 fAmpt->SetIpop(fIpop);
215 fAmpt->SetXmu(fXmu);
a004b331 216 fAmpt->SetAlpha(fAlpha);
217 fAmpt->SetStringFrag(fStringA, fStringB);
0119ef9a 218
219 AliGenMC::Init();
220
221 // Initialize Ampt
222 fAmpt->Initialize();
223 if (fEvaluate)
224 EvaluateCrossSections();
225}
226
227void AliGenAmpt::Generate()
228{
229 // Generate one event
230
231 Float_t polar[3] = {0,0,0};
232 Float_t origin[3] = {0,0,0};
233 Float_t origin0[3] = {0,0,0};
234 Float_t p[3];
235 Float_t tof;
236
237 // converts from mm/c to s
238 const Float_t kconv = 0.001/2.99792458e8;
239
240 Int_t nt = 0;
241 Int_t jev = 0;
242 Int_t j, kf, ks, ksp, imo;
243 kf = 0;
244
245 fTrials = 0;
246 for (j = 0;j < 3; j++)
247 origin0[j] = fOrigin[j];
248
249 if(fVertexSmear == kPerEvent) {
250 Vertex();
251 for (j=0; j < 3; j++)
252 origin0[j] = fVertex[j];
253 }
254
255 Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
256
257 while(1) {
258 // Generate one event
259 Int_t fpemask = gSystem->GetFPEMask();
260 gSystem->SetFPEMask(0);
261 fAmpt->GenerateEvent();
262 gSystem->SetFPEMask(fpemask);
263 fTrials++;
264 fNprimaries = 0;
265 fAmpt->ImportParticles(&fParticles,"All");
266 Int_t np = fParticles.GetEntriesFast();
267 if (np == 0 )
268 continue;
269
270 if (fTrigger != kNoTrigger) {
271 if (!CheckTrigger())
272 continue;
273 }
274 if (fLHC)
275 Boost();
276
277 Int_t nc = 0;
278 Int_t* newPos = new Int_t[np];
279 Int_t* pSelected = new Int_t[np];
280
281 for (Int_t i = 0; i < np; i++) {
282 newPos[i] = i;
283 pSelected[i] = 0;
284 }
285
286 // Get event vertex
287 //TParticle * iparticle = (TParticle *) fParticles.At(0);
288 fVertex[0] = origin0[0];
289 fVertex[1] = origin0[1];
290 fVertex[2] = origin0[2];
291
292 // First select parent particles
293 for (Int_t i = 0; i < np; i++) {
294 TParticle *iparticle = (TParticle *) fParticles.At(i);
295
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();
302 if (kf == 92)
303 continue;
304
305 if (!fSelectAll)
306 selected = KinematicSelection(iparticle, 0) && SelectFlavor(kf);
307 hasSelectedDaughters = DaughtersSelection(iparticle);
308
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) {
312 nc++;
313 pSelected[i] = 1;
314 } // selected
315 } // particle loop parents
316
317 // Now select the final state particles
318 fProjectileSpecn = 0;
319 fProjectileSpecp = 0;
320 fTargetSpecn = 0;
321 fTargetSpecp = 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))
326 continue;
327 Bool_t selected = kTRUE;
328 kf = iparticle->GetPdgCode();
329 ks = iparticle->GetStatusCode();
330 ksp = iparticle->GetUniqueID();
331
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;
340 }
341 // --------------------------------------------------------------------------
342 if (!fSelectAll) {
343 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
344 if (!fSpectators && selected)
345 selected = (ksp != 0 && ksp != 1 && ksp != 10 && ksp != 11);
346 }
347
348 // Put particle on the stack if selected
349 if (selected) {
350 nc++;
351 pSelected[i] = 1;
352 } // selected
353 } // particle loop final state
354
355 //Time of the interactions
356 Float_t tInt = 0.;
357 if (fPileUpTimeWindow > 0.)
358 tInt = fPileUpTimeWindow * (2. * gRandom->Rndm() - 1.);
359
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);
365 if (pSelected[i]) {
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;
374 fEventTime = 0.;
375
376 if (TestBit(kVertexRange)) {
377 fEventTime = sign * origin0[2] / 2.99792458e10;
378 tof = kconv * iparticle->T() + fEventTime;
379 } else {
380 tof = kconv * iparticle->T();
381 fEventTime = tInt;
382 if (fPileUpTimeWindow > 0.) tof += tInt;
383 }
384 imo = -1;
385 TParticle* mother = 0;
386 if (hasMother) {
387 imo = iparticle->GetFirstMother();
388 mother = (TParticle *) fParticles.At(imo);
389 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
390 } // if has mother
391 Bool_t tFlag = (fTrackIt && !hasDaughter);
392 PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
393 fNprimaries++;
394 KeepTrack(nt);
395 newPos[i] = nt;
396 } // if selected
397 } // particle loop
398 delete[] newPos;
399 delete[] pSelected;
400
401 AliInfo(Form("\n I've put %i particles on the stack \n",nc));
402 if (nc > 0) {
403 jev += 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));
407 break;
408 }
409 }
410 } // event loop
411 MakeHeader();
412 SetHighWaterMark(nt);
413}
414
415void AliGenAmpt::EvaluateCrossSections()
416{
417 // Glauber Calculation of geometrical x-section
418
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
424 Float_t bMin = 0.;
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;
428
429 printf("\n Projectile Radius (fm): %f \n",fAmpt->GetHIPR1(34));
430 printf("\n Target Radius (fm): %f \n",fAmpt->GetHIPR1(35));
431
432 Int_t i;
433 Float_t oldvalue= 0.;
434 Float_t* b = new Float_t[kMax];
435 Float_t* si1 = new Float_t[kMax];
436 Float_t* si2 = new Float_t[kMax];
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;
442 xTot+=gb;
443 xTotHard += gbh;
444 printf("profile %f %f %f\n", xb, ov, fAmpt->GetHINT1(12));
445
446 if (xb > fMinImpactParam && xb < fMaxImpactParam) {
447 xPart += gb;
448 xPartHard += gbh;
449 }
450
451 if ((oldvalue) && ((xTot-oldvalue)/oldvalue<0.0001))
452 break;
453 oldvalue = xTot;
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);
456 if (i>0) {
457 si1[i] = gb/kdib;
458 si2[i] = gbh/gb;
459 b[i] = xb;
460 }
461 }
462
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.);
467
468 // Store result as a graph
469 b[0] = 0;
470 si1[0] = 0;
471 si2[0]=si2[1];
472 delete fDsigmaDb;
473 fDsigmaDb = new TGraph(i, b, si1);
474 delete fDnDb;
475 fDnDb = new TGraph(i, b, si2);
476}
477
478Bool_t AliGenAmpt::DaughtersSelection(TParticle* iparticle)
479{
480 // Looks recursively if one of the daughters has been selected
481 //printf("\n Consider daughters %d:",iparticle->GetPdgCode());
482 Int_t imin = -1;
483 Int_t imax = -1;
484 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
485 Bool_t selected = kFALSE;
486 if (hasDaughters) {
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;
494 }
495 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
496 }
497 } else {
498 return kFALSE;
499 }
500 return selected;
501}
502
503Bool_t AliGenAmpt::SelectFlavor(Int_t pid)
504{
505 // Select flavor of particle
506 // 0: all
507 // 4: charm and beauty
508 // 5: beauty
509 Bool_t res = 0;
510
511 if (fFlavor == 0) {
512 res = kTRUE;
513 } else {
514 Int_t ifl = TMath::Abs(pid/100);
515 if (ifl > 10) ifl/=10;
516 res = (fFlavor == ifl);
517 }
518
519 // This part if gamma writing is inhibited
520 if (fNoGammas)
521 res = res && (pid != kGamma && pid != kPi0);
522
523 return res;
524}
525
526Bool_t AliGenAmpt::Stable(TParticle* particle) const
527{
528 // Return true for a stable particle
529
530 if (particle->GetFirstDaughter() < 0 )
531 {
532 return kTRUE;
533 } else {
534 return kFALSE;
535 }
536}
537
538void AliGenAmpt::MakeHeader()
539{
540 // Fills the event header, to be called after each event
541
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(),
548 fAmpt->GetN01(),
549 fAmpt->GetN10(),
550 fAmpt->GetN11());
551 fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
552 fTargetSpecn,fTargetSpecp);
553 fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
554 //printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
555
556 // 4-momentum vectors of the triggered jets.
557 // Before final state gluon radiation.
558 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(21),
559 fAmpt->GetHINT1(22),
560 fAmpt->GetHINT1(23),
561 fAmpt->GetHINT1(24));
562
563 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(31),
564 fAmpt->GetHINT1(32),
565 fAmpt->GetHINT1(33),
566 fAmpt->GetHINT1(34));
567 // After final state gluon radiation.
568 TLorentzVector* jet3 = new TLorentzVector(fAmpt->GetHINT1(26),
569 fAmpt->GetHINT1(27),
570 fAmpt->GetHINT1(28),
571 fAmpt->GetHINT1(29));
572
573 TLorentzVector* jet4 = new TLorentzVector(fAmpt->GetHINT1(36),
574 fAmpt->GetHINT1(37),
575 fAmpt->GetHINT1(38),
576 fAmpt->GetHINT1(39));
577 fHeader->SetJets(jet1, jet2, jet3, jet4);
578 // Bookkeeping for kinematic bias
579 fHeader->SetTrials(fTrials);
580 // Event Vertex
581 fHeader->SetPrimaryVertex(fVertex);
582 fHeader->SetInteractionTime(fEventTime);
583
584 fCollisionGeometry = fHeader;
585 AddHeader(fHeader);
586}
587
588
589Bool_t AliGenAmpt::CheckTrigger()
590{
591 // Check the kinematic trigger condition
592
593 Bool_t triggered = kFALSE;
594
595 if (fTrigger == 1) {
596 // jet-jet Trigger
597 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(26),
598 fAmpt->GetHINT1(27),
599 fAmpt->GetHINT1(28),
600 fAmpt->GetHINT1(29));
601
602 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(36),
603 fAmpt->GetHINT1(37),
604 fAmpt->GetHINT1(38),
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)
613 ||
614 (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&
615 phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
616 )
617 triggered = kTRUE;
618 } else if (fTrigger == 2) {
619 // Gamma Jet
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 &&
629 eta > fEtaMinJet &&
630 phi < fPhiMaxJet &&
631 phi > fPhiMinJet) {
632 triggered = 1;
633 break;
634 } // check phi,eta within limits
635 } // direct gamma ?
636 } // particle loop
637 } // fTrigger == 2
638 return triggered;
639}