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