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