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