]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TAmpt/AliGenAmpt.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
do not modify esdTrack outerParam in memory
[u/mrichter/AliRoot.git] / TAmpt / AliGenAmpt.cxx
... / ...
CommitLineData
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#include <TVirtualMC.h>
29#include <TParticlePDG.h>
30#include "AliGenHijingEventHeader.h"
31#define AliGenAmptEventHeader AliGenHijingEventHeader
32#include "AliAmptRndm.h"
33#include "AliLog.h"
34#include "AliRun.h"
35#include "AliDecayer.h"
36
37ClassImp(AliGenAmpt)
38
39AliGenAmpt::AliGenAmpt()
40 : AliGenMC(),
41 fDecayer(NULL),
42 fFrame("CMS"),
43 fMinImpactParam(0.),
44 fMaxImpactParam(5.),
45 fKeep(0),
46 fQuench(0),
47 fShadowing(1),
48 fDecaysOff(1),
49 fTrigger(0),
50 fEvaluate(0),
51 fSelectAll(0),
52 fFlavor(0),
53 fKineBias(0.),
54 fTrials(0),
55 fXsection(0.),
56 fAmpt(0),
57 fPtHardMin(2.0),
58 fPtHardMax(-1),
59 fSpectators(1),
60 fDsigmaDb(0),
61 fDnDb(0),
62 fPtMinJet(-2.5),
63 fEtaMinJet(-20.),
64 fEtaMaxJet(+20.),
65 fPhiMinJet(0.),
66 fPhiMaxJet(TMath::TwoPi()),
67 fRadiation(3),
68 fSimpleJet(kFALSE),
69 fNoGammas(kFALSE),
70 fProjectileSpecn(0),
71 fProjectileSpecp(0),
72 fTargetSpecn(0),
73 fTargetSpecp(0),
74 fLHC(kFALSE),
75 fRandomPz(kFALSE),
76 fNoHeavyQuarks(kFALSE),
77 fIsoft(4),
78 fNtMax(150),
79 fIpop(1),
80 fXmu(3.2264),
81 fAlpha(1./3),
82 fStringA(0.5),
83 fStringB(0.9),
84 fEventTime(0.),
85 fHeader(new AliGenAmptEventHeader("Ampt")),
86 fDecay(kTRUE)
87{
88 // Constructor
89 fEnergyCMS = 2760.;
90 AliAmptRndm::SetAmptRandom(GetRandom());
91}
92
93AliGenAmpt::AliGenAmpt(Int_t npart)
94 : AliGenMC(npart),
95 fDecayer(NULL),
96 fFrame("CMS"),
97 fMinImpactParam(0.),
98 fMaxImpactParam(5.),
99 fKeep(0),
100 fQuench(0),
101 fShadowing(1),
102 fDecaysOff(1),
103 fTrigger(0),
104 fEvaluate(0),
105 fSelectAll(0),
106 fFlavor(0),
107 fKineBias(0.),
108 fTrials(0),
109 fXsection(0.),
110 fAmpt(0),
111 fPtHardMin(2.0),
112 fPtHardMax(-1),
113 fSpectators(1),
114 fDsigmaDb(0),
115 fDnDb(0),
116 fPtMinJet(-2.5),
117 fEtaMinJet(-20.),
118 fEtaMaxJet(+20.),
119 fPhiMinJet(0.),
120 fPhiMaxJet(2. * TMath::Pi()),
121 fRadiation(3),
122 fSimpleJet(kFALSE),
123 fNoGammas(kFALSE),
124 fProjectileSpecn(0),
125 fProjectileSpecp(0),
126 fTargetSpecn(0),
127 fTargetSpecp(0),
128 fLHC(kFALSE),
129 fRandomPz(kFALSE),
130 fNoHeavyQuarks(kFALSE),
131 fIsoft(1),
132 fNtMax(150),
133 fIpop(1),
134 fXmu(3.2264),
135 fAlpha(1./3),
136 fStringA(0.5),
137 fStringB(0.9),
138 fEventTime(0.),
139 fHeader(new AliGenAmptEventHeader("Ampt")),
140 fDecay(kTRUE)
141{
142 // Default PbPb collisions at 2.76 TeV
143
144 fEnergyCMS = 2760.;
145 fName = "Ampt";
146 fTitle= "Particle Generator using AMPT";
147 AliAmptRndm::SetAmptRandom(GetRandom());
148}
149
150AliGenAmpt::~AliGenAmpt()
151{
152 // Destructor
153 if ( fDsigmaDb) delete fDsigmaDb;
154 if ( fDnDb) delete fDnDb;
155 if ( fHeader) delete fHeader;
156}
157
158void AliGenAmpt::Init()
159{
160 // Initialisation
161
162 fFrame.Resize(8);
163 fTarget.Resize(8);
164 fProjectile.Resize(8);
165
166 fAmpt = new TAmpt(fEnergyCMS, fFrame, fProjectile, fTarget,
167 fAProjectile, fZProjectile, fATarget, fZTarget,
168 fMinImpactParam, fMaxImpactParam);
169 SetMC(fAmpt);
170
171 fAmpt->SetIHPR2(2, fRadiation);
172 fAmpt->SetIHPR2(3, fTrigger);
173 fAmpt->SetIHPR2(6, fShadowing);
174 fAmpt->SetIHPR2(12, fDecaysOff);
175 fAmpt->SetIHPR2(21, fKeep);
176 fAmpt->SetHIPR1(8, fPtHardMin);
177 fAmpt->SetHIPR1(9, fPtHardMax);
178 fAmpt->SetHIPR1(10, fPtMinJet);
179 fAmpt->SetHIPR1(50, fSimpleJet);
180
181 // Quenching
182 // fQuench = 0: no quenching
183 // fQuench = 1: Hijing default
184 // fQuench = 2: new LHC parameters for HIPR1(11) and HIPR1(14)
185 // fQuench = 3: new RHIC parameters for HIPR1(11) and HIPR1(14)
186 // fQuench = 4: new LHC parameters with log(e) dependence
187 // fQuench = 5: new RHIC parameters with log(e) dependence
188 fAmpt->SetIHPR2(50, 0);
189 if (fQuench > 0)
190 fAmpt->SetIHPR2(4, 1);
191 else
192 fAmpt->SetIHPR2(4, 0);
193
194 if (fQuench == 2) {
195 fAmpt->SetHIPR1(14, 1.1);
196 fAmpt->SetHIPR1(11, 3.7);
197 } else if (fQuench == 3) {
198 fAmpt->SetHIPR1(14, 0.20);
199 fAmpt->SetHIPR1(11, 2.5);
200 } else if (fQuench == 4) {
201 fAmpt->SetIHPR2(50, 1);
202 fAmpt->SetHIPR1(14, 4.*0.34);
203 fAmpt->SetHIPR1(11, 3.7);
204 } else if (fQuench == 5) {
205 fAmpt->SetIHPR2(50, 1);
206 fAmpt->SetHIPR1(14, 0.34);
207 fAmpt->SetHIPR1(11, 2.5);
208 }
209
210 // Heavy quarks
211 if (fNoHeavyQuarks) {
212 fAmpt->SetIHPR2(49, 1);
213 } else {
214 fAmpt->SetIHPR2(49, 0);
215 }
216
217 // Ampt specific
218 fAmpt->SetIsoft(fIsoft);
219 fAmpt->SetNtMax(fNtMax);
220 fAmpt->SetIpop(fIpop);
221 fAmpt->SetXmu(fXmu);
222 fAmpt->SetAlpha(fAlpha);
223 fAmpt->SetStringFrag(fStringA, fStringB);
224
225 AliGenMC::Init();
226
227 // Initialize Ampt
228 fAmpt->Initialize();
229 if (fEvaluate)
230 EvaluateCrossSections();
231}
232
233void AliGenAmpt::Generate()
234{
235 // Generate one event
236
237 Float_t polar[3] = {0,0,0};
238 Float_t origin[3] = {0,0,0};
239 Float_t origin0[3] = {0,0,0};
240 Float_t time0 = 0.;
241 Float_t p[3];
242 Float_t tof;
243
244 // converts from mm/c to s
245 const Float_t kconv = 0.001/2.99792458e8;
246
247 Int_t nt = 0;
248 Int_t jev = 0;
249 Int_t j, kf, ks, ksp, imo;
250 kf = 0;
251
252 fTrials = 0;
253 for (j = 0;j < 3; j++)
254 origin0[j] = fOrigin[j];
255 //time0 = fTimeOrigin;
256
257 if(fVertexSmear == kPerEvent) {
258 Vertex();
259 for (j=0; j < 3; j++)
260 origin0[j] = fVertex[j];
261 //time0 = fTime;
262 }
263
264 Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
265
266 while(1) {
267 // Generate one event
268 Int_t fpemask = gSystem->GetFPEMask();
269 gSystem->SetFPEMask(0);
270 fAmpt->GenerateEvent();
271 gSystem->SetFPEMask(fpemask);
272 fTrials++;
273 fNprimaries = 0;
274 fAmpt->ImportParticles(&fParticles,"All");
275 Int_t np = fParticles.GetEntriesFast();
276 if (np == 0 )
277 continue;
278
279 if (fTrigger != kNoTrigger) {
280 if (!CheckTrigger())
281 continue;
282 }
283
284 AliDecayer *decayer = 0;
285 //if (gMC)
286 // decayer = gMC->GetDecayer();
287 decayer = fDecayer; //AMPT does not do the strong decays per dafault
288
289 if (decayer&&fDecay) {
290 TClonesArray arr("TParticle",100);
291 for( Int_t nLoop=0; nLoop!=2; ++nLoop) { // In order to produce more than one generation of decays: NumberOfNestedLoops set to 2
292 Int_t np2 = np;
293 for (Int_t i = 0; i < np; i++) {
294 TParticle *iparticle = (TParticle *)fParticles.At(i);
295 if (!Stable(iparticle)) // true if particle has daughters already
296 continue;
297 kf = TMath::Abs(iparticle->GetPdgCode());
298 if (kf==92)
299 continue;
300 if( !IsThisAKnownParticle(iparticle) ) continue; // skip undesired particles
301 /*
302 if (0) { // this turned out to be too cumbersome!
303 if (kf!=331&&kf!=3114&&kf!=3114&&kf!=411&&kf!=-4122&&kf!=-3324&&kf!=-3312&&kf!=-3114&&
304 kf!=-311&&kf!=3214&&kf!=-3214&&kf!=-433&&kf!=413&&kf!=3122&&kf!=-3122&&kf!=-413&&
305 kf!=-421&&kf!=-423&&kf!=3324&&kf!=-313&&kf!=213&&kf!=-213&&kf!=3314&&kf!=3222&&
306 kf!=-3222&&kf!=3224&&kf!=-3224&&kf!=-4212&&kf!=4212&&kf!=433&&kf!=423&&kf!=-3322&&
307 kf!=3322&&kf!=-3314)
308 continue; //decay eta',Sigma*+,Sigma*-,D+,Lambda_c-,Xi*0_bar,Xi-_bar,Sigma*-,
309 // K0_bar,Sigma*0,Sigma*0_bar,D*_s-,D*+,Lambda0,Lambda0_bar,D*-
310 // D0_bar,D*0_bar,Xi*0,K*0_bar,rho+,rho-,Xi*-,Sigma-,
311 // Sigma+,Sigma*+,Sigma*-,Sigma_c-,Sigma_c+,D*_s+,D*0,Xi0_bar
312 // Xi0,Xi*+
313 //} else { // really only decay particles if there are not known to Geant3
314 // if (gMC->IdFromPDG(kf)>0)
315 // continue;
316 }
317 if (0) { // defining the particle for Geant3 leads to a floating point exception.
318 TParticlePDG *pdg = iparticle->GetPDG(1);
319 //pdg->Print(); printf("%s\n",pdg->ParticleClass());
320 TString ptype(pdg->ParticleClass());
321 TMCParticleType mctype(kPTUndefined);
322 if (ptype=="Baryon" || ptype=="Meson")
323 mctype = kPTHadron;
324 gMC->DefineParticle(pdg->PdgCode(), pdg->GetName(), mctype, pdg->Mass(), pdg->Charge(), pdg->Lifetime(),
325 ptype,pdg->Width(), (Int_t)pdg->Spin(), (Int_t)pdg->Parity(), 0,
326 (Int_t)pdg->Isospin(), 0, 0, 0, 0, pdg->Stable());
327 gMC->SetUserDecay(pdg->PdgCode());
328 continue;
329 }
330 */
331 TLorentzVector pmom(iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy());
332 decayer->Decay(kf,&pmom);
333 decayer->ImportParticles(&arr);
334 Int_t ndecayed = arr.GetEntries();
335 if (ndecayed>1) {
336 if (np2+ndecayed>fParticles.GetSize())
337 fParticles.Expand(2*fParticles.GetSize());
338 //arr.Print();
339 // iparticle->SetStatusCode(2); to be compatible with Hijing
340 iparticle->SetFirstDaughter(np2);
341 for (Int_t jj = 1; jj < ndecayed; jj++) {
342 TParticle *jp = (TParticle *)arr.At(jj);
343 if (jp->GetFirstMother()!=1)
344 continue;
345 TParticle *newp = new(fParticles[np2]) TParticle(jp->GetPdgCode(),
346 0, //1, //to be compatible with Hijing
347 i,
348 -1,
349 -1,
350 -1,
351 jp->Px(),jp->Py(),jp->Pz(),jp->Energy(),
352 jp->Vx(),jp->Vy(),jp->Vz(),jp->T());
353 newp->SetUniqueID( jp->GetStatusCode() );
354 np2++;
355 } // end of jj->nDecayedParticles
356 iparticle->SetLastDaughter(np2-1);
357 } // end of nDecayedPrticles>1
358 } // end of i->np
359 np = fParticles.GetEntries();
360 if (np!=np2) {
361 AliError(Form("Something is fishy: %d %d\n", np,np2));
362 }
363 } // end of nLoop->NumberOfNestedLoops
364 } else {
365 if (fDecay)
366 AliError("No decayer found, but fDecay==kTRUE!");
367 }
368
369 if (fLHC)
370 Boost();
371
372 Int_t nc = 0;
373 Int_t* newPos = new Int_t[np];
374 Int_t* pSelected = new Int_t[np];
375
376 for (Int_t i = 0; i < np; i++) {
377 newPos[i] = i;
378 pSelected[i] = 0;
379 }
380
381 // Get event vertex
382 //TParticle * iparticle = (TParticle *) fParticles.At(0);
383 fVertex[0] = origin0[0];
384 fVertex[1] = origin0[1];
385 fVertex[2] = origin0[2];
386 //fTime = time0;
387
388 // First select parent particles
389 for (Int_t i = 0; i < np; i++) {
390 TParticle *iparticle = (TParticle *) fParticles.At(i);
391
392 // Is this a parent particle ?
393 if (Stable(iparticle)) continue; // quit if particle has no daughters
394 Bool_t selected = kTRUE;
395 Bool_t hasSelectedDaughters = kFALSE;
396 kf = iparticle->GetPdgCode();
397 ks = iparticle->GetStatusCode();
398 if (kf == 92)
399 continue;
400
401 if (!fSelectAll)
402 selected = KinematicSelection(iparticle, 0) && SelectFlavor(kf);
403 hasSelectedDaughters = DaughtersSelection(iparticle);
404
405 // Put particle on the stack if it is either selected or
406 // it is the mother of at least one seleted particle
407 if (selected || hasSelectedDaughters) {
408 nc++;
409 pSelected[i] = 1;
410 } // selected
411 } // particle loop parents
412
413 // Now select the final state particles
414 fProjectileSpecn = 0;
415 fProjectileSpecp = 0;
416 fTargetSpecn = 0;
417 fTargetSpecp = 0;
418 for (Int_t i = 0; i<np; i++) {
419 TParticle *iparticle = (TParticle *) fParticles.At(i);
420 // Is this a final state particle ?
421 if (!Stable(iparticle)) continue; // quit if particle has daughters
422 Bool_t selected = kTRUE;
423 kf = iparticle->GetPdgCode();
424 if (kf == 92)
425 continue;
426 ks = iparticle->GetStatusCode();
427 ksp = iparticle->GetUniqueID();
428
429 // --------------------------------------------------------------------------
430 // Count spectator neutrons and protons
431 if(ksp == 0 || ksp == 1) {
432 if(kf == kNeutron) fProjectileSpecn += 1;
433 if(kf == kProton) fProjectileSpecp += 1;
434 } else if(ksp == 10 || ksp == 11) {
435 if(kf == kNeutron) fTargetSpecn += 1;
436 if(kf == kProton) fTargetSpecp += 1;
437 }
438 // --------------------------------------------------------------------------
439 if (!fSelectAll) {
440 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
441 if (!fSpectators && selected)
442 selected = (ksp != 0 && ksp != 1 && ksp != 10 && ksp != 11);
443 }
444
445 // Put particle on the stack if selected
446 if (selected) {
447 nc++;
448 pSelected[i] = 1;
449 if (0) printf("---> %d %d %d %s\n",i,nc,kf,iparticle->GetName());
450 } // selected
451 } // particle loop final state
452
453 // Write particles to stack
454 for (Int_t i = 0; i<np; i++) {
455 if (pSelected[i]) {
456 TParticle *iparticle = (TParticle *) fParticles.At(i);
457 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
458 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
459 kf = iparticle->GetPdgCode();
460 ks = iparticle->GetStatusCode();
461 p[0] = iparticle->Px();
462 p[1] = iparticle->Py();
463 p[2] = iparticle->Pz() * sign;
464 origin[0] = origin0[0]+iparticle->Vx()/10;
465 origin[1] = origin0[1]+iparticle->Vy()/10;
466 origin[2] = origin0[2]+iparticle->Vz()/10;
467 tof = time0+kconv * iparticle->T();
468
469 imo = -1;
470 TParticle* mother = 0;
471 TMCProcess procID = (TMCProcess) iparticle->GetUniqueID();
472 if (hasMother) {
473 imo = iparticle->GetFirstMother();
474 mother = (TParticle *) fParticles.At(imo);
475 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
476 } else { // if has no mothers then it was created by AMPT
477 if(procID==999)
478 procID = kPPrimary; // reseting to ALIROOT convention
479 else
480 procID = kPNoProcess; // for expectators
481 } // if has mother
482 Bool_t tFlag = (fTrackIt && !hasDaughter);
483 PushTrack(tFlag,imo,kf,p,origin,polar,tof,procID,nt, 1., ks);
484 fNprimaries++;
485 KeepTrack(nt);
486 newPos[i] = nt;
487 } // if selected
488 } // particle loop
489 delete[] newPos;
490 delete[] pSelected;
491
492 AliInfo(Form("\n I've put %i particles on the stack \n",nc));
493 if (nc > 0) {
494 jev += nc;
495 if (jev >= fNpart || fNpart == -1) {
496 fKineBias = Float_t(fNpart)/Float_t(fTrials);
497 AliInfo(Form("\n Trials: %i %i %i\n",fTrials, fNpart, jev));
498 break;
499 }
500 }
501 } // event loop
502 MakeHeader();
503 SetHighWaterMark(nt);
504}
505
506Bool_t AliGenAmpt::IsThisAKnownParticle(TParticle *thisGuy)
507{
508 // In order to prevent AMPT to introduce weird particles into the decayer and transporter
509 // blame cperez@cern.ch for this method
510
511 Int_t pdgcode = TMath::Abs( thisGuy->GetPdgCode() );
512
513 Int_t myFavoriteParticles[ 38] = { 3322, 3314, 3312, 3224, 3222, // Xi0 Xi*+- Xi+- Sigma*-+ Sigma-+
514 3214, 3212, 3122, 3114, 3112, // Sigma*0 Sigma0 Lambda0 Sigma*+- Sigma+-
515 2224, 2214, 2212, 2114, 2112, // Delta--++ Delta-+ proton Delta0 neutron
516 1114, 323, 321, 313, 311, // Delta+- K*-+ K-+ K*0 K0
517 213, 211, 11, 22, 111, // rho-+ pi-+ e+- gamma pi0
518 113, 130, 221, 223, 310, // rho0 K_L0 eta omega K_S0
519 331, 333, 3324, 431, 421, // eta' phi Xi*0 Ds-+ D0
520 411, 413, 13 // D-+ D*-+ mu+-
521 };
522
523 Bool_t found = kFALSE;
524 for(Int_t i=0; i!=38; ++i)
525 if( myFavoriteParticles[i] == pdgcode ) {
526 found = kTRUE;
527 break;
528 }
529
530 return found;
531}
532
533void AliGenAmpt::EvaluateCrossSections()
534{
535 // Glauber Calculation of geometrical x-section
536
537 Float_t xTot = 0.; // barn
538 Float_t xTotHard = 0.; // barn
539 Float_t xPart = 0.; // barn
540 Float_t xPartHard = 0.; // barn
541 Float_t sigmaHard = 0.1; // mbarn
542 Float_t bMin = 0.;
543 Float_t bMax = fAmpt->GetHIPR1(34)+fAmpt->GetHIPR1(35);
544 const Float_t kdib = 0.2;
545 Int_t kMax = Int_t((bMax-bMin)/kdib)+1;
546
547 printf("\n Projectile Radius (fm): %f \n",fAmpt->GetHIPR1(34));
548 printf("\n Target Radius (fm): %f \n",fAmpt->GetHIPR1(35));
549
550 Int_t i;
551 Float_t oldvalue= 0.;
552 Float_t* b = new Float_t[kMax]; memset(b,0,kMax*sizeof(Float_t));
553 Float_t* si1 = new Float_t[kMax]; memset(si1,0,kMax*sizeof(Float_t));
554 Float_t* si2 = new Float_t[kMax]; memset(si2,0,kMax*sizeof(Float_t));
555 for (i = 0; i < kMax; i++) {
556 Float_t xb = bMin+i*kdib;
557 Float_t ov=fAmpt->Profile(xb);
558 Float_t gb = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fAmpt->GetHINT1(12)*ov));
559 Float_t gbh = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
560 xTot+=gb;
561 xTotHard += gbh;
562 printf("profile %f %f %f\n", xb, ov, fAmpt->GetHINT1(12));
563
564 if (xb > fMinImpactParam && xb < fMaxImpactParam) {
565 xPart += gb;
566 xPartHard += gbh;
567 }
568
569 if ((oldvalue) && ((xTot-oldvalue)/oldvalue<0.0001))
570 break;
571 oldvalue = xTot;
572 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
573 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
574 if (i>0) {
575 si1[i] = gb/kdib;
576 si2[i] = gbh/gb;
577 b[i] = xb;
578 }
579 }
580
581 printf("\n Total cross section (barn): %f \n",xTot);
582 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
583 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
584 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
585
586 // Store result as a graph
587 b[0] = 0;
588 si1[0] = 0;
589 si2[0]=si2[1];
590 delete fDsigmaDb;
591 fDsigmaDb = new TGraph(i, b, si1);
592 delete fDnDb;
593 fDnDb = new TGraph(i, b, si2);
594}
595
596Bool_t AliGenAmpt::DaughtersSelection(TParticle* iparticle)
597{
598 // Looks recursively if one of the daughters has been selected
599 //printf("\n Consider daughters %d:",iparticle->GetPdgCode());
600 Int_t imin = -1;
601 Int_t imax = -1;
602 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
603 Bool_t selected = kFALSE;
604 if (hasDaughters) {
605 imin = iparticle->GetFirstDaughter();
606 imax = iparticle->GetLastDaughter();
607 for (Int_t i = imin; i <= imax; i++){
608 TParticle * jparticle = (TParticle *) fParticles.At(i);
609 Int_t ip = jparticle->GetPdgCode();
610 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
611 selected=kTRUE; break;
612 }
613 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
614 }
615 } else {
616 return kFALSE;
617 }
618 return selected;
619}
620
621Bool_t AliGenAmpt::SelectFlavor(Int_t pid)
622{
623 // Select flavor of particle
624 // 0: all
625 // 4: charm and beauty
626 // 5: beauty
627 Bool_t res = 0;
628
629 if (fFlavor == 0) {
630 res = kTRUE;
631 } else {
632 Int_t ifl = TMath::Abs(pid/100);
633 if (ifl > 10) ifl/=10;
634 res = (fFlavor == ifl);
635 }
636
637 // This part if gamma writing is inhibited
638 if (fNoGammas)
639 res = res && (pid != kGamma && pid != kPi0);
640
641 return res;
642}
643
644Bool_t AliGenAmpt::Stable(TParticle* particle) const
645{
646 // Return true for a stable particle
647
648 if (!particle)
649 return kFALSE;
650 if (particle->GetFirstDaughter() < 0 )
651 return kTRUE;
652 return kFALSE;
653
654 /// ADD LIST
655
656}
657
658void AliGenAmpt::MakeHeader()
659{
660 // Fills the event header, to be called after each event
661
662 fHeader->SetNProduced(fNprimaries);
663 fHeader->SetImpactParameter(fAmpt->GetHINT1(19));
664 fHeader->SetTotalEnergy(fAmpt->GetEATT());
665 fHeader->SetHardScatters(fAmpt->GetJATT());
666 fHeader->SetParticipants(fAmpt->GetNP(), fAmpt->GetNT());
667 fHeader->SetCollisions(fAmpt->GetN0(),
668 fAmpt->GetN01(),
669 fAmpt->GetN10(),
670 fAmpt->GetN11());
671 fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
672 fTargetSpecn,fTargetSpecp);
673 fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
674 //printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
675
676 // 4-momentum vectors of the triggered jets.
677 // Before final state gluon radiation.
678 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(21),
679 fAmpt->GetHINT1(22),
680 fAmpt->GetHINT1(23),
681 fAmpt->GetHINT1(24));
682
683 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(31),
684 fAmpt->GetHINT1(32),
685 fAmpt->GetHINT1(33),
686 fAmpt->GetHINT1(34));
687 // After final state gluon radiation.
688 TLorentzVector* jet3 = new TLorentzVector(fAmpt->GetHINT1(26),
689 fAmpt->GetHINT1(27),
690 fAmpt->GetHINT1(28),
691 fAmpt->GetHINT1(29));
692
693 TLorentzVector* jet4 = new TLorentzVector(fAmpt->GetHINT1(36),
694 fAmpt->GetHINT1(37),
695 fAmpt->GetHINT1(38),
696 fAmpt->GetHINT1(39));
697 fHeader->SetJets(jet1, jet2, jet3, jet4);
698 // Bookkeeping for kinematic bias
699 fHeader->SetTrials(fTrials);
700 // Event Vertex
701 fHeader->SetPrimaryVertex(fVertex);
702 fHeader->SetInteractionTime(fEventTime);
703
704 fCollisionGeometry = fHeader;
705 AddHeader(fHeader);
706}
707
708
709Bool_t AliGenAmpt::CheckTrigger()
710{
711 // Check the kinematic trigger condition
712
713 Bool_t triggered = kFALSE;
714
715 if (fTrigger == 1) {
716 // jet-jet Trigger
717 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(26),
718 fAmpt->GetHINT1(27),
719 fAmpt->GetHINT1(28),
720 fAmpt->GetHINT1(29));
721
722 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(36),
723 fAmpt->GetHINT1(37),
724 fAmpt->GetHINT1(38),
725 fAmpt->GetHINT1(39));
726 Double_t eta1 = jet1->Eta();
727 Double_t eta2 = jet2->Eta();
728 Double_t phi1 = jet1->Phi();
729 Double_t phi2 = jet2->Phi();
730 //printf("\n Trigger: %f %f %f %f", fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
731 if ( (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&
732 phi1 < fPhiMaxJet && phi1 > fPhiMinJet)
733 ||
734 (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&
735 phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
736 )
737 triggered = kTRUE;
738 } else if (fTrigger == 2) {
739 // Gamma Jet
740 Int_t np = fParticles.GetEntriesFast();
741 for (Int_t i = 0; i < np; i++) {
742 TParticle* part = (TParticle*) fParticles.At(i);
743 Int_t kf = part->GetPdgCode();
744 Int_t ksp = part->GetUniqueID();
745 if (kf == 22 && ksp == 40) {
746 Float_t phi = part->Phi();
747 Float_t eta = part->Eta();
748 if (eta < fEtaMaxJet &&
749 eta > fEtaMinJet &&
750 phi < fPhiMaxJet &&
751 phi > fPhiMinJet) {
752 triggered = 1;
753 break;
754 } // check phi,eta within limits
755 } // direct gamma ?
756 } // particle loop
757 } // fTrigger == 2
758 return triggered;
759}