]> 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
During simulation: fill STU region w/ non null time sums
[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(),
d12123e8 41 fDecayer(NULL),
0119ef9a 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),
c6db63eb 77 fIsoft(4),
0119ef9a 78 fNtMax(150),
79 fIpop(1),
80 fXmu(3.2264),
a004b331 81 fAlpha(1./3),
82 fStringA(0.5),
83 fStringB(0.9),
d12123e8 84 fEventTime(0.),
d60fac1e 85 fHeader(new AliGenAmptEventHeader("Ampt")),
86 fDecay(kTRUE)
0119ef9a 87{
88 // Constructor
d60fac1e 89 fEnergyCMS = 2760.;
0119ef9a 90 AliAmptRndm::SetAmptRandom(GetRandom());
91}
92
93AliGenAmpt::AliGenAmpt(Int_t npart)
94 : AliGenMC(npart),
d12123e8 95 fDecayer(NULL),
0119ef9a 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),
0119ef9a 131 fIsoft(1),
132 fNtMax(150),
133 fIpop(1),
134 fXmu(3.2264),
a004b331 135 fAlpha(1./3),
136 fStringA(0.5),
137 fStringB(0.9),
d12123e8 138 fEventTime(0.),
d60fac1e 139 fHeader(new AliGenAmptEventHeader("Ampt")),
140 fDecay(kTRUE)
0119ef9a 141{
d60fac1e 142 // Default PbPb collisions at 2.76 TeV
0119ef9a 143
d60fac1e 144 fEnergyCMS = 2760.;
0119ef9a 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);
a004b331 222 fAmpt->SetAlpha(fAlpha);
223 fAmpt->SetStringFrag(fStringA, fStringB);
0119ef9a 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};
cdfa7be6 240 Float_t time0 = 0.;
0119ef9a 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];
d12123e8 255 //time0 = fTimeOrigin;
0119ef9a 256
257 if(fVertexSmear == kPerEvent) {
258 Vertex();
259 for (j=0; j < 3; j++)
260 origin0[j] = fVertex[j];
d12123e8 261 //time0 = fTime;
0119ef9a 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 }
d60fac1e 283
284 AliDecayer *decayer = 0;
d12123e8 285 //if (gMC)
286 // decayer = gMC->GetDecayer();
287 decayer = fDecayer; //AMPT does not do the strong decays per dafault
d60fac1e 288
289 if (decayer&&fDecay) {
290 TClonesArray arr("TParticle",100);
d12123e8 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));
d60fac1e 362 }
d12123e8 363 } // end of nLoop->NumberOfNestedLoops
d60fac1e 364 } else {
365 if (fDecay)
366 AliError("No decayer found, but fDecay==kTRUE!");
367 }
368
0119ef9a 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];
d12123e8 386 //fTime = time0;
0119ef9a 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 ?
d12123e8 393 if (Stable(iparticle)) continue; // quit if particle has no daughters
0119ef9a 394 Bool_t selected = kTRUE;
395 Bool_t hasSelectedDaughters = kFALSE;
d12123e8 396 kf = iparticle->GetPdgCode();
397 ks = iparticle->GetStatusCode();
0119ef9a 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 ?
d12123e8 421 if (!Stable(iparticle)) continue; // quit if particle has daughters
422 Bool_t selected = kTRUE;
423 kf = iparticle->GetPdgCode();
d60fac1e 424 if (kf == 92)
425 continue;
d12123e8 426 ks = iparticle->GetStatusCode();
427 ksp = iparticle->GetUniqueID();
0119ef9a 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;
68a6e988 449 if (0) printf("---> %d %d %d %s\n",i,nc,kf,iparticle->GetName());
0119ef9a 450 } // selected
451 } // particle loop final state
452
0119ef9a 453 // Write particles to stack
454 for (Int_t i = 0; i<np; i++) {
0119ef9a 455 if (pSelected[i]) {
d12123e8 456 TParticle *iparticle = (TParticle *) fParticles.At(i);
457 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
458 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
0119ef9a 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;
cdfa7be6 467 tof = time0+kconv * iparticle->T();
468
0119ef9a 469 imo = -1;
470 TParticle* mother = 0;
d12123e8 471 TMCProcess procID = (TMCProcess) iparticle->GetUniqueID();
0119ef9a 472 if (hasMother) {
473 imo = iparticle->GetFirstMother();
474 mother = (TParticle *) fParticles.At(imo);
475 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
d12123e8 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
0119ef9a 481 } // if has mother
482 Bool_t tFlag = (fTrackIt && !hasDaughter);
d12123e8 483 PushTrack(tFlag,imo,kf,p,origin,polar,tof,procID,nt, 1., ks);
0119ef9a 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
d12123e8 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
0119ef9a 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.;
48eb126c 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));
0119ef9a 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
d60fac1e 644Bool_t AliGenAmpt::Stable(TParticle* particle) const
0119ef9a 645{
646 // Return true for a stable particle
647
d60fac1e 648 if (!particle)
649 return kFALSE;
0119ef9a 650 if (particle->GetFirstDaughter() < 0 )
0119ef9a 651 return kTRUE;
d60fac1e 652 return kFALSE;
d12123e8 653
654 /// ADD LIST
655
0119ef9a 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);
d12123e8 702 fHeader->SetInteractionTime(fEventTime);
703
0119ef9a 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}