]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TAmpt/AliGenAmpt.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added new triggers
[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),
79 fHeader(new AliGenAmptEventHeader("Ampt"))
80{
81 // Constructor
82 fEnergyCMS = 5500.;
83 AliAmptRndm::SetAmptRandom(GetRandom());
84}
85
86AliGenAmpt::AliGenAmpt(Int_t npart)
87 : AliGenMC(npart),
88 fFrame("CMS"),
89 fMinImpactParam(0.),
90 fMaxImpactParam(5.),
91 fKeep(0),
92 fQuench(0),
93 fShadowing(1),
94 fDecaysOff(1),
95 fTrigger(0),
96 fEvaluate(0),
97 fSelectAll(0),
98 fFlavor(0),
99 fKineBias(0.),
100 fTrials(0),
101 fXsection(0.),
102 fAmpt(0),
103 fPtHardMin(2.0),
104 fPtHardMax(-1),
105 fSpectators(1),
106 fDsigmaDb(0),
107 fDnDb(0),
108 fPtMinJet(-2.5),
109 fEtaMinJet(-20.),
110 fEtaMaxJet(+20.),
111 fPhiMinJet(0.),
112 fPhiMaxJet(2. * TMath::Pi()),
113 fRadiation(3),
114 fSimpleJet(kFALSE),
115 fNoGammas(kFALSE),
116 fProjectileSpecn(0),
117 fProjectileSpecp(0),
118 fTargetSpecn(0),
119 fTargetSpecp(0),
120 fLHC(kFALSE),
121 fRandomPz(kFALSE),
122 fNoHeavyQuarks(kFALSE),
123 fEventTime(0.),
124 fIsoft(1),
125 fNtMax(150),
126 fIpop(1),
127 fXmu(3.2264),
128 fHeader(new AliGenAmptEventHeader("Ampt"))
129{
130 // Default PbPb collisions at 5.5 TeV
131
132 fEnergyCMS = 5500.;
133 fName = "Ampt";
134 fTitle= "Particle Generator using AMPT";
135 AliAmptRndm::SetAmptRandom(GetRandom());
136}
137
138AliGenAmpt::~AliGenAmpt()
139{
140 // Destructor
141 if ( fDsigmaDb) delete fDsigmaDb;
142 if ( fDnDb) delete fDnDb;
143 if ( fHeader) delete fHeader;
144}
145
146void AliGenAmpt::Init()
147{
148 // Initialisation
149
150 fFrame.Resize(8);
151 fTarget.Resize(8);
152 fProjectile.Resize(8);
153
154 fAmpt = new TAmpt(fEnergyCMS, fFrame, fProjectile, fTarget,
155 fAProjectile, fZProjectile, fATarget, fZTarget,
156 fMinImpactParam, fMaxImpactParam);
157 SetMC(fAmpt);
158
159 fAmpt->SetIHPR2(2, fRadiation);
160 fAmpt->SetIHPR2(3, fTrigger);
161 fAmpt->SetIHPR2(6, fShadowing);
162 fAmpt->SetIHPR2(12, fDecaysOff);
163 fAmpt->SetIHPR2(21, fKeep);
164 fAmpt->SetHIPR1(8, fPtHardMin);
165 fAmpt->SetHIPR1(9, fPtHardMax);
166 fAmpt->SetHIPR1(10, fPtMinJet);
167 fAmpt->SetHIPR1(50, fSimpleJet);
168
169 // Quenching
170 // fQuench = 0: no quenching
171 // fQuench = 1: Hijing default
172 // fQuench = 2: new LHC parameters for HIPR1(11) and HIPR1(14)
173 // fQuench = 3: new RHIC parameters for HIPR1(11) and HIPR1(14)
174 // fQuench = 4: new LHC parameters with log(e) dependence
175 // fQuench = 5: new RHIC parameters with log(e) dependence
176 fAmpt->SetIHPR2(50, 0);
177 if (fQuench > 0)
178 fAmpt->SetIHPR2(4, 1);
179 else
180 fAmpt->SetIHPR2(4, 0);
181
182 if (fQuench == 2) {
183 fAmpt->SetHIPR1(14, 1.1);
184 fAmpt->SetHIPR1(11, 3.7);
185 } else if (fQuench == 3) {
186 fAmpt->SetHIPR1(14, 0.20);
187 fAmpt->SetHIPR1(11, 2.5);
188 } else if (fQuench == 4) {
189 fAmpt->SetIHPR2(50, 1);
190 fAmpt->SetHIPR1(14, 4.*0.34);
191 fAmpt->SetHIPR1(11, 3.7);
192 } else if (fQuench == 5) {
193 fAmpt->SetIHPR2(50, 1);
194 fAmpt->SetHIPR1(14, 0.34);
195 fAmpt->SetHIPR1(11, 2.5);
196 }
197
198 // Heavy quarks
199 if (fNoHeavyQuarks) {
200 fAmpt->SetIHPR2(49, 1);
201 } else {
202 fAmpt->SetIHPR2(49, 0);
203 }
204
205 // Ampt specific
206 fAmpt->SetIsoft(fIsoft);
207 fAmpt->SetNtMax(fNtMax);
208 fAmpt->SetIpop(fIpop);
209 fAmpt->SetXmu(fXmu);
210
211 AliGenMC::Init();
212
213 // Initialize Ampt
214 fAmpt->Initialize();
215 if (fEvaluate)
216 EvaluateCrossSections();
217}
218
219void AliGenAmpt::Generate()
220{
221 // Generate one event
222
223 Float_t polar[3] = {0,0,0};
224 Float_t origin[3] = {0,0,0};
225 Float_t origin0[3] = {0,0,0};
226 Float_t p[3];
227 Float_t tof;
228
229 // converts from mm/c to s
230 const Float_t kconv = 0.001/2.99792458e8;
231
232 Int_t nt = 0;
233 Int_t jev = 0;
234 Int_t j, kf, ks, ksp, imo;
235 kf = 0;
236
237 fTrials = 0;
238 for (j = 0;j < 3; j++)
239 origin0[j] = fOrigin[j];
240
241 if(fVertexSmear == kPerEvent) {
242 Vertex();
243 for (j=0; j < 3; j++)
244 origin0[j] = fVertex[j];
245 }
246
247 Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
248
249 while(1) {
250 // Generate one event
251 Int_t fpemask = gSystem->GetFPEMask();
252 gSystem->SetFPEMask(0);
253 fAmpt->GenerateEvent();
254 gSystem->SetFPEMask(fpemask);
255 fTrials++;
256 fNprimaries = 0;
257 fAmpt->ImportParticles(&fParticles,"All");
258 Int_t np = fParticles.GetEntriesFast();
259 if (np == 0 )
260 continue;
261
262 if (fTrigger != kNoTrigger) {
263 if (!CheckTrigger())
264 continue;
265 }
266 if (fLHC)
267 Boost();
268
269 Int_t nc = 0;
270 Int_t* newPos = new Int_t[np];
271 Int_t* pSelected = new Int_t[np];
272
273 for (Int_t i = 0; i < np; i++) {
274 newPos[i] = i;
275 pSelected[i] = 0;
276 }
277
278 // Get event vertex
279 //TParticle * iparticle = (TParticle *) fParticles.At(0);
280 fVertex[0] = origin0[0];
281 fVertex[1] = origin0[1];
282 fVertex[2] = origin0[2];
283
284 // First select parent particles
285 for (Int_t i = 0; i < np; i++) {
286 TParticle *iparticle = (TParticle *) fParticles.At(i);
287
288 // Is this a parent particle ?
289 if (Stable(iparticle)) continue;
290 Bool_t selected = kTRUE;
291 Bool_t hasSelectedDaughters = kFALSE;
292 kf = iparticle->GetPdgCode();
293 ks = iparticle->GetStatusCode();
294 if (kf == 92)
295 continue;
296
297 if (!fSelectAll)
298 selected = KinematicSelection(iparticle, 0) && SelectFlavor(kf);
299 hasSelectedDaughters = DaughtersSelection(iparticle);
300
301 // Put particle on the stack if it is either selected or
302 // it is the mother of at least one seleted particle
303 if (selected || hasSelectedDaughters) {
304 nc++;
305 pSelected[i] = 1;
306 } // selected
307 } // particle loop parents
308
309 // Now select the final state particles
310 fProjectileSpecn = 0;
311 fProjectileSpecp = 0;
312 fTargetSpecn = 0;
313 fTargetSpecp = 0;
314 for (Int_t i = 0; i<np; i++) {
315 TParticle *iparticle = (TParticle *) fParticles.At(i);
316 // Is this a final state particle ?
317 if (!Stable(iparticle))
318 continue;
319 Bool_t selected = kTRUE;
320 kf = iparticle->GetPdgCode();
321 ks = iparticle->GetStatusCode();
322 ksp = iparticle->GetUniqueID();
323
324 // --------------------------------------------------------------------------
325 // Count spectator neutrons and protons
326 if(ksp == 0 || ksp == 1) {
327 if(kf == kNeutron) fProjectileSpecn += 1;
328 if(kf == kProton) fProjectileSpecp += 1;
329 } else if(ksp == 10 || ksp == 11) {
330 if(kf == kNeutron) fTargetSpecn += 1;
331 if(kf == kProton) fTargetSpecp += 1;
332 }
333 // --------------------------------------------------------------------------
334 if (!fSelectAll) {
335 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
336 if (!fSpectators && selected)
337 selected = (ksp != 0 && ksp != 1 && ksp != 10 && ksp != 11);
338 }
339
340 // Put particle on the stack if selected
341 if (selected) {
342 nc++;
343 pSelected[i] = 1;
344 } // selected
345 } // particle loop final state
346
347 //Time of the interactions
348 Float_t tInt = 0.;
349 if (fPileUpTimeWindow > 0.)
350 tInt = fPileUpTimeWindow * (2. * gRandom->Rndm() - 1.);
351
352 // Write particles to stack
353 for (Int_t i = 0; i<np; i++) {
354 TParticle *iparticle = (TParticle *) fParticles.At(i);
355 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
356 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
357 if (pSelected[i]) {
358 kf = iparticle->GetPdgCode();
359 ks = iparticle->GetStatusCode();
360 p[0] = iparticle->Px();
361 p[1] = iparticle->Py();
362 p[2] = iparticle->Pz() * sign;
363 origin[0] = origin0[0]+iparticle->Vx()/10;
364 origin[1] = origin0[1]+iparticle->Vy()/10;
365 origin[2] = origin0[2]+iparticle->Vz()/10;
366 fEventTime = 0.;
367
368 if (TestBit(kVertexRange)) {
369 fEventTime = sign * origin0[2] / 2.99792458e10;
370 tof = kconv * iparticle->T() + fEventTime;
371 } else {
372 tof = kconv * iparticle->T();
373 fEventTime = tInt;
374 if (fPileUpTimeWindow > 0.) tof += tInt;
375 }
376 imo = -1;
377 TParticle* mother = 0;
378 if (hasMother) {
379 imo = iparticle->GetFirstMother();
380 mother = (TParticle *) fParticles.At(imo);
381 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
382 } // if has mother
383 Bool_t tFlag = (fTrackIt && !hasDaughter);
384 PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
385 fNprimaries++;
386 KeepTrack(nt);
387 newPos[i] = nt;
388 } // if selected
389 } // particle loop
390 delete[] newPos;
391 delete[] pSelected;
392
393 AliInfo(Form("\n I've put %i particles on the stack \n",nc));
394 if (nc > 0) {
395 jev += nc;
396 if (jev >= fNpart || fNpart == -1) {
397 fKineBias = Float_t(fNpart)/Float_t(fTrials);
398 AliInfo(Form("\n Trials: %i %i %i\n",fTrials, fNpart, jev));
399 break;
400 }
401 }
402 } // event loop
403 MakeHeader();
404 SetHighWaterMark(nt);
405}
406
407void AliGenAmpt::EvaluateCrossSections()
408{
409 // Glauber Calculation of geometrical x-section
410
411 Float_t xTot = 0.; // barn
412 Float_t xTotHard = 0.; // barn
413 Float_t xPart = 0.; // barn
414 Float_t xPartHard = 0.; // barn
415 Float_t sigmaHard = 0.1; // mbarn
416 Float_t bMin = 0.;
417 Float_t bMax = fAmpt->GetHIPR1(34)+fAmpt->GetHIPR1(35);
418 const Float_t kdib = 0.2;
419 Int_t kMax = Int_t((bMax-bMin)/kdib)+1;
420
421 printf("\n Projectile Radius (fm): %f \n",fAmpt->GetHIPR1(34));
422 printf("\n Target Radius (fm): %f \n",fAmpt->GetHIPR1(35));
423
424 Int_t i;
425 Float_t oldvalue= 0.;
426 Float_t* b = new Float_t[kMax];
427 Float_t* si1 = new Float_t[kMax];
428 Float_t* si2 = new Float_t[kMax];
429 for (i = 0; i < kMax; i++) {
430 Float_t xb = bMin+i*kdib;
431 Float_t ov=fAmpt->Profile(xb);
432 Float_t gb = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fAmpt->GetHINT1(12)*ov));
433 Float_t gbh = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
434 xTot+=gb;
435 xTotHard += gbh;
436 printf("profile %f %f %f\n", xb, ov, fAmpt->GetHINT1(12));
437
438 if (xb > fMinImpactParam && xb < fMaxImpactParam) {
439 xPart += gb;
440 xPartHard += gbh;
441 }
442
443 if ((oldvalue) && ((xTot-oldvalue)/oldvalue<0.0001))
444 break;
445 oldvalue = xTot;
446 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
447 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
448 if (i>0) {
449 si1[i] = gb/kdib;
450 si2[i] = gbh/gb;
451 b[i] = xb;
452 }
453 }
454
455 printf("\n Total cross section (barn): %f \n",xTot);
456 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
457 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
458 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
459
460 // Store result as a graph
461 b[0] = 0;
462 si1[0] = 0;
463 si2[0]=si2[1];
464 delete fDsigmaDb;
465 fDsigmaDb = new TGraph(i, b, si1);
466 delete fDnDb;
467 fDnDb = new TGraph(i, b, si2);
468}
469
470Bool_t AliGenAmpt::DaughtersSelection(TParticle* iparticle)
471{
472 // Looks recursively if one of the daughters has been selected
473 //printf("\n Consider daughters %d:",iparticle->GetPdgCode());
474 Int_t imin = -1;
475 Int_t imax = -1;
476 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
477 Bool_t selected = kFALSE;
478 if (hasDaughters) {
479 imin = iparticle->GetFirstDaughter();
480 imax = iparticle->GetLastDaughter();
481 for (Int_t i = imin; i <= imax; i++){
482 TParticle * jparticle = (TParticle *) fParticles.At(i);
483 Int_t ip = jparticle->GetPdgCode();
484 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
485 selected=kTRUE; break;
486 }
487 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
488 }
489 } else {
490 return kFALSE;
491 }
492 return selected;
493}
494
495Bool_t AliGenAmpt::SelectFlavor(Int_t pid)
496{
497 // Select flavor of particle
498 // 0: all
499 // 4: charm and beauty
500 // 5: beauty
501 Bool_t res = 0;
502
503 if (fFlavor == 0) {
504 res = kTRUE;
505 } else {
506 Int_t ifl = TMath::Abs(pid/100);
507 if (ifl > 10) ifl/=10;
508 res = (fFlavor == ifl);
509 }
510
511 // This part if gamma writing is inhibited
512 if (fNoGammas)
513 res = res && (pid != kGamma && pid != kPi0);
514
515 return res;
516}
517
518Bool_t AliGenAmpt::Stable(TParticle* particle) const
519{
520 // Return true for a stable particle
521
522 if (particle->GetFirstDaughter() < 0 )
523 {
524 return kTRUE;
525 } else {
526 return kFALSE;
527 }
528}
529
530void AliGenAmpt::MakeHeader()
531{
532 // Fills the event header, to be called after each event
533
534 fHeader->SetNProduced(fNprimaries);
535 fHeader->SetImpactParameter(fAmpt->GetHINT1(19));
536 fHeader->SetTotalEnergy(fAmpt->GetEATT());
537 fHeader->SetHardScatters(fAmpt->GetJATT());
538 fHeader->SetParticipants(fAmpt->GetNP(), fAmpt->GetNT());
539 fHeader->SetCollisions(fAmpt->GetN0(),
540 fAmpt->GetN01(),
541 fAmpt->GetN10(),
542 fAmpt->GetN11());
543 fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
544 fTargetSpecn,fTargetSpecp);
545 fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
546 //printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
547
548 // 4-momentum vectors of the triggered jets.
549 // Before final state gluon radiation.
550 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(21),
551 fAmpt->GetHINT1(22),
552 fAmpt->GetHINT1(23),
553 fAmpt->GetHINT1(24));
554
555 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(31),
556 fAmpt->GetHINT1(32),
557 fAmpt->GetHINT1(33),
558 fAmpt->GetHINT1(34));
559 // After final state gluon radiation.
560 TLorentzVector* jet3 = new TLorentzVector(fAmpt->GetHINT1(26),
561 fAmpt->GetHINT1(27),
562 fAmpt->GetHINT1(28),
563 fAmpt->GetHINT1(29));
564
565 TLorentzVector* jet4 = new TLorentzVector(fAmpt->GetHINT1(36),
566 fAmpt->GetHINT1(37),
567 fAmpt->GetHINT1(38),
568 fAmpt->GetHINT1(39));
569 fHeader->SetJets(jet1, jet2, jet3, jet4);
570 // Bookkeeping for kinematic bias
571 fHeader->SetTrials(fTrials);
572 // Event Vertex
573 fHeader->SetPrimaryVertex(fVertex);
574 fHeader->SetInteractionTime(fEventTime);
575
576 fCollisionGeometry = fHeader;
577 AddHeader(fHeader);
578}
579
580
581Bool_t AliGenAmpt::CheckTrigger()
582{
583 // Check the kinematic trigger condition
584
585 Bool_t triggered = kFALSE;
586
587 if (fTrigger == 1) {
588 // jet-jet Trigger
589 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(26),
590 fAmpt->GetHINT1(27),
591 fAmpt->GetHINT1(28),
592 fAmpt->GetHINT1(29));
593
594 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(36),
595 fAmpt->GetHINT1(37),
596 fAmpt->GetHINT1(38),
597 fAmpt->GetHINT1(39));
598 Double_t eta1 = jet1->Eta();
599 Double_t eta2 = jet2->Eta();
600 Double_t phi1 = jet1->Phi();
601 Double_t phi2 = jet2->Phi();
602 //printf("\n Trigger: %f %f %f %f", fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
603 if ( (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&
604 phi1 < fPhiMaxJet && phi1 > fPhiMinJet)
605 ||
606 (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&
607 phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
608 )
609 triggered = kTRUE;
610 } else if (fTrigger == 2) {
611 // Gamma Jet
612 Int_t np = fParticles.GetEntriesFast();
613 for (Int_t i = 0; i < np; i++) {
614 TParticle* part = (TParticle*) fParticles.At(i);
615 Int_t kf = part->GetPdgCode();
616 Int_t ksp = part->GetUniqueID();
617 if (kf == 22 && ksp == 40) {
618 Float_t phi = part->Phi();
619 Float_t eta = part->Eta();
620 if (eta < fEtaMaxJet &&
621 eta > fEtaMinJet &&
622 phi < fPhiMaxJet &&
623 phi > fPhiMinJet) {
624 triggered = 1;
625 break;
626 } // check phi,eta within limits
627 } // direct gamma ?
628 } // particle loop
629 } // fTrigger == 2
630 return triggered;
631}