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