]> 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
Fix to access the efficiency map file
[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);
afe6642c 362 for (Int_t jj = 1; jj < ndecayed; jj++) {
d12123e8 363 TParticle *jp = (TParticle *)arr.At(jj);
364 if (jp->GetFirstMother()!=1)
365 continue;
afe6642c 366
367 TParticle *newp = new(fParticles[np2]) TParticle(jp->GetPdgCode(),
d12123e8 368 0, //1, //to be compatible with Hijing
369 i,
370 -1,
371 -1,
372 -1,
373 jp->Px(),jp->Py(),jp->Pz(),jp->Energy(),
afe6642c 374 jp->Vx(),jp->Vy(),jp->Vz(),jp->T());
375 //take care of the phi
376 //if((kf == 333)||(kf == 313)) {
377 if(IsThisAKnownParticle(iparticle)) {
378 //Printf("=============PANOS===================");
379 //Printf("Phi detected - daughet is: %d",jp->GetPdgCode());
380 newp->SetUniqueID(4);
381 }
382 else newp->SetUniqueID( jp->GetStatusCode() );
d12123e8 383 np2++;
384 } // end of jj->nDecayedParticles
385 iparticle->SetLastDaughter(np2-1);
386 } // end of nDecayedPrticles>1
387 } // end of i->np
388 np = fParticles.GetEntries();
389 if (np!=np2) {
390 AliError(Form("Something is fishy: %d %d\n", np,np2));
d60fac1e 391 }
d12123e8 392 } // end of nLoop->NumberOfNestedLoops
d60fac1e 393 } else {
394 if (fDecay)
395 AliError("No decayer found, but fDecay==kTRUE!");
396 }
397
0119ef9a 398 if (fLHC)
399 Boost();
400
401 Int_t nc = 0;
402 Int_t* newPos = new Int_t[np];
403 Int_t* pSelected = new Int_t[np];
404
405 for (Int_t i = 0; i < np; i++) {
406 newPos[i] = i;
407 pSelected[i] = 0;
408 }
409
410 // Get event vertex
411 //TParticle * iparticle = (TParticle *) fParticles.At(0);
412 fVertex[0] = origin0[0];
413 fVertex[1] = origin0[1];
414 fVertex[2] = origin0[2];
d12123e8 415 //fTime = time0;
0119ef9a 416
417 // First select parent particles
418 for (Int_t i = 0; i < np; i++) {
419 TParticle *iparticle = (TParticle *) fParticles.At(i);
420
421 // Is this a parent particle ?
d12123e8 422 if (Stable(iparticle)) continue; // quit if particle has no daughters
0119ef9a 423 Bool_t selected = kTRUE;
424 Bool_t hasSelectedDaughters = kFALSE;
d12123e8 425 kf = iparticle->GetPdgCode();
426 ks = iparticle->GetStatusCode();
0119ef9a 427 if (kf == 92)
428 continue;
429
430 if (!fSelectAll)
431 selected = KinematicSelection(iparticle, 0) && SelectFlavor(kf);
432 hasSelectedDaughters = DaughtersSelection(iparticle);
433
434 // Put particle on the stack if it is either selected or
435 // it is the mother of at least one seleted particle
436 if (selected || hasSelectedDaughters) {
437 nc++;
438 pSelected[i] = 1;
439 } // selected
440 } // particle loop parents
441
442 // Now select the final state particles
443 fProjectileSpecn = 0;
444 fProjectileSpecp = 0;
445 fTargetSpecn = 0;
446 fTargetSpecp = 0;
447 for (Int_t i = 0; i<np; i++) {
448 TParticle *iparticle = (TParticle *) fParticles.At(i);
449 // Is this a final state particle ?
d12123e8 450 if (!Stable(iparticle)) continue; // quit if particle has daughters
451 Bool_t selected = kTRUE;
452 kf = iparticle->GetPdgCode();
d60fac1e 453 if (kf == 92)
454 continue;
d12123e8 455 ks = iparticle->GetStatusCode();
456 ksp = iparticle->GetUniqueID();
0119ef9a 457
458 // --------------------------------------------------------------------------
459 // Count spectator neutrons and protons
460 if(ksp == 0 || ksp == 1) {
461 if(kf == kNeutron) fProjectileSpecn += 1;
462 if(kf == kProton) fProjectileSpecp += 1;
463 } else if(ksp == 10 || ksp == 11) {
464 if(kf == kNeutron) fTargetSpecn += 1;
465 if(kf == kProton) fTargetSpecp += 1;
466 }
467 // --------------------------------------------------------------------------
468 if (!fSelectAll) {
469 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
470 if (!fSpectators && selected)
471 selected = (ksp != 0 && ksp != 1 && ksp != 10 && ksp != 11);
472 }
473
474 // Put particle on the stack if selected
475 if (selected) {
476 nc++;
477 pSelected[i] = 1;
68a6e988 478 if (0) printf("---> %d %d %d %s\n",i,nc,kf,iparticle->GetName());
0119ef9a 479 } // selected
480 } // particle loop final state
481
0119ef9a 482 // Write particles to stack
483 for (Int_t i = 0; i<np; i++) {
0119ef9a 484 if (pSelected[i]) {
d12123e8 485 TParticle *iparticle = (TParticle *) fParticles.At(i);
486 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
487 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
0119ef9a 488 kf = iparticle->GetPdgCode();
489 ks = iparticle->GetStatusCode();
490 p[0] = iparticle->Px();
491 p[1] = iparticle->Py();
492 p[2] = iparticle->Pz() * sign;
774ceaaf 493 origin[0] = origin0[0]+iparticle->Vx();
494 origin[1] = origin0[1]+iparticle->Vy();
495 origin[2] = origin0[2]+iparticle->Vz();
496 tof = time0 + iparticle->T();
cdfa7be6 497
0119ef9a 498 imo = -1;
499 TParticle* mother = 0;
d12123e8 500 TMCProcess procID = (TMCProcess) iparticle->GetUniqueID();
0119ef9a 501 if (hasMother) {
502 imo = iparticle->GetFirstMother();
503 mother = (TParticle *) fParticles.At(imo);
504 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
d12123e8 505 } else { // if has no mothers then it was created by AMPT
506 if(procID==999)
507 procID = kPPrimary; // reseting to ALIROOT convention
508 else
509 procID = kPNoProcess; // for expectators
0119ef9a 510 } // if has mother
511 Bool_t tFlag = (fTrackIt && !hasDaughter);
d12123e8 512 PushTrack(tFlag,imo,kf,p,origin,polar,tof,procID,nt, 1., ks);
0119ef9a 513 fNprimaries++;
514 KeepTrack(nt);
515 newPos[i] = nt;
516 } // if selected
517 } // particle loop
518 delete[] newPos;
519 delete[] pSelected;
520
521 AliInfo(Form("\n I've put %i particles on the stack \n",nc));
522 if (nc > 0) {
523 jev += nc;
524 if (jev >= fNpart || fNpart == -1) {
525 fKineBias = Float_t(fNpart)/Float_t(fTrials);
526 AliInfo(Form("\n Trials: %i %i %i\n",fTrials, fNpart, jev));
527 break;
528 }
529 }
530 } // event loop
531 MakeHeader();
532 SetHighWaterMark(nt);
533}
534
d12123e8 535Bool_t AliGenAmpt::IsThisAKnownParticle(TParticle *thisGuy)
536{
537 // In order to prevent AMPT to introduce weird particles into the decayer and transporter
538 // blame cperez@cern.ch for this method
539
540 Int_t pdgcode = TMath::Abs( thisGuy->GetPdgCode() );
541
542 Int_t myFavoriteParticles[ 38] = { 3322, 3314, 3312, 3224, 3222, // Xi0 Xi*+- Xi+- Sigma*-+ Sigma-+
543 3214, 3212, 3122, 3114, 3112, // Sigma*0 Sigma0 Lambda0 Sigma*+- Sigma+-
544 2224, 2214, 2212, 2114, 2112, // Delta--++ Delta-+ proton Delta0 neutron
545 1114, 323, 321, 313, 311, // Delta+- K*-+ K-+ K*0 K0
546 213, 211, 11, 22, 111, // rho-+ pi-+ e+- gamma pi0
547 113, 130, 221, 223, 310, // rho0 K_L0 eta omega K_S0
548 331, 333, 3324, 431, 421, // eta' phi Xi*0 Ds-+ D0
549 411, 413, 13 // D-+ D*-+ mu+-
550 };
551
552 Bool_t found = kFALSE;
553 for(Int_t i=0; i!=38; ++i)
554 if( myFavoriteParticles[i] == pdgcode ) {
555 found = kTRUE;
556 break;
557 }
558
559 return found;
560}
561
0119ef9a 562void AliGenAmpt::EvaluateCrossSections()
563{
564 // Glauber Calculation of geometrical x-section
565
566 Float_t xTot = 0.; // barn
567 Float_t xTotHard = 0.; // barn
568 Float_t xPart = 0.; // barn
569 Float_t xPartHard = 0.; // barn
570 Float_t sigmaHard = 0.1; // mbarn
571 Float_t bMin = 0.;
572 Float_t bMax = fAmpt->GetHIPR1(34)+fAmpt->GetHIPR1(35);
573 const Float_t kdib = 0.2;
574 Int_t kMax = Int_t((bMax-bMin)/kdib)+1;
575
576 printf("\n Projectile Radius (fm): %f \n",fAmpt->GetHIPR1(34));
577 printf("\n Target Radius (fm): %f \n",fAmpt->GetHIPR1(35));
578
579 Int_t i;
580 Float_t oldvalue= 0.;
48eb126c 581 Float_t* b = new Float_t[kMax]; memset(b,0,kMax*sizeof(Float_t));
582 Float_t* si1 = new Float_t[kMax]; memset(si1,0,kMax*sizeof(Float_t));
583 Float_t* si2 = new Float_t[kMax]; memset(si2,0,kMax*sizeof(Float_t));
0119ef9a 584 for (i = 0; i < kMax; i++) {
585 Float_t xb = bMin+i*kdib;
586 Float_t ov=fAmpt->Profile(xb);
587 Float_t gb = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fAmpt->GetHINT1(12)*ov));
588 Float_t gbh = 2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
589 xTot+=gb;
590 xTotHard += gbh;
591 printf("profile %f %f %f\n", xb, ov, fAmpt->GetHINT1(12));
592
593 if (xb > fMinImpactParam && xb < fMaxImpactParam) {
594 xPart += gb;
595 xPartHard += gbh;
596 }
597
598 if ((oldvalue) && ((xTot-oldvalue)/oldvalue<0.0001))
599 break;
600 oldvalue = xTot;
601 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
602 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
603 if (i>0) {
604 si1[i] = gb/kdib;
605 si2[i] = gbh/gb;
606 b[i] = xb;
607 }
608 }
609
610 printf("\n Total cross section (barn): %f \n",xTot);
611 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
612 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
613 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
614
615 // Store result as a graph
616 b[0] = 0;
617 si1[0] = 0;
618 si2[0]=si2[1];
619 delete fDsigmaDb;
620 fDsigmaDb = new TGraph(i, b, si1);
621 delete fDnDb;
622 fDnDb = new TGraph(i, b, si2);
623}
624
625Bool_t AliGenAmpt::DaughtersSelection(TParticle* iparticle)
626{
627 // Looks recursively if one of the daughters has been selected
628 //printf("\n Consider daughters %d:",iparticle->GetPdgCode());
629 Int_t imin = -1;
630 Int_t imax = -1;
631 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
632 Bool_t selected = kFALSE;
633 if (hasDaughters) {
634 imin = iparticle->GetFirstDaughter();
635 imax = iparticle->GetLastDaughter();
636 for (Int_t i = imin; i <= imax; i++){
637 TParticle * jparticle = (TParticle *) fParticles.At(i);
638 Int_t ip = jparticle->GetPdgCode();
639 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
640 selected=kTRUE; break;
641 }
642 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
643 }
644 } else {
645 return kFALSE;
646 }
647 return selected;
648}
649
650Bool_t AliGenAmpt::SelectFlavor(Int_t pid)
651{
652 // Select flavor of particle
653 // 0: all
654 // 4: charm and beauty
655 // 5: beauty
656 Bool_t res = 0;
657
658 if (fFlavor == 0) {
659 res = kTRUE;
660 } else {
661 Int_t ifl = TMath::Abs(pid/100);
662 if (ifl > 10) ifl/=10;
663 res = (fFlavor == ifl);
664 }
665
666 // This part if gamma writing is inhibited
667 if (fNoGammas)
668 res = res && (pid != kGamma && pid != kPi0);
669
670 return res;
671}
672
d60fac1e 673Bool_t AliGenAmpt::Stable(TParticle* particle) const
0119ef9a 674{
675 // Return true for a stable particle
676
d60fac1e 677 if (!particle)
678 return kFALSE;
0119ef9a 679 if (particle->GetFirstDaughter() < 0 )
0119ef9a 680 return kTRUE;
d60fac1e 681 return kFALSE;
d12123e8 682
683 /// ADD LIST
684
0119ef9a 685}
686
687void AliGenAmpt::MakeHeader()
688{
689 // Fills the event header, to be called after each event
690
691 fHeader->SetNProduced(fNprimaries);
692 fHeader->SetImpactParameter(fAmpt->GetHINT1(19));
693 fHeader->SetTotalEnergy(fAmpt->GetEATT());
694 fHeader->SetHardScatters(fAmpt->GetJATT());
695 fHeader->SetParticipants(fAmpt->GetNP(), fAmpt->GetNT());
696 fHeader->SetCollisions(fAmpt->GetN0(),
697 fAmpt->GetN01(),
698 fAmpt->GetN10(),
699 fAmpt->GetN11());
700 fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
701 fTargetSpecn,fTargetSpecp);
48beeea0 702 //fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
703 fHeader->SetReactionPlaneAngle(fAmpt->GetReactionPlaneAngle());
0119ef9a 704 //printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
705
706 // 4-momentum vectors of the triggered jets.
707 // Before final state gluon radiation.
708 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(21),
709 fAmpt->GetHINT1(22),
710 fAmpt->GetHINT1(23),
711 fAmpt->GetHINT1(24));
712
713 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(31),
714 fAmpt->GetHINT1(32),
715 fAmpt->GetHINT1(33),
716 fAmpt->GetHINT1(34));
717 // After final state gluon radiation.
718 TLorentzVector* jet3 = new TLorentzVector(fAmpt->GetHINT1(26),
719 fAmpt->GetHINT1(27),
720 fAmpt->GetHINT1(28),
721 fAmpt->GetHINT1(29));
722
723 TLorentzVector* jet4 = new TLorentzVector(fAmpt->GetHINT1(36),
724 fAmpt->GetHINT1(37),
725 fAmpt->GetHINT1(38),
726 fAmpt->GetHINT1(39));
727 fHeader->SetJets(jet1, jet2, jet3, jet4);
728 // Bookkeeping for kinematic bias
729 fHeader->SetTrials(fTrials);
730 // Event Vertex
731 fHeader->SetPrimaryVertex(fVertex);
d12123e8 732 fHeader->SetInteractionTime(fEventTime);
733
0119ef9a 734 fCollisionGeometry = fHeader;
735 AddHeader(fHeader);
736}
737
738
739Bool_t AliGenAmpt::CheckTrigger()
740{
741 // Check the kinematic trigger condition
742
743 Bool_t triggered = kFALSE;
744
745 if (fTrigger == 1) {
746 // jet-jet Trigger
747 TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(26),
748 fAmpt->GetHINT1(27),
749 fAmpt->GetHINT1(28),
750 fAmpt->GetHINT1(29));
751
752 TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(36),
753 fAmpt->GetHINT1(37),
754 fAmpt->GetHINT1(38),
755 fAmpt->GetHINT1(39));
756 Double_t eta1 = jet1->Eta();
757 Double_t eta2 = jet2->Eta();
758 Double_t phi1 = jet1->Phi();
759 Double_t phi2 = jet2->Phi();
760 //printf("\n Trigger: %f %f %f %f", fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
761 if ( (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&
762 phi1 < fPhiMaxJet && phi1 > fPhiMinJet)
763 ||
764 (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&
765 phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
766 )
767 triggered = kTRUE;
768 } else if (fTrigger == 2) {
769 // Gamma Jet
770 Int_t np = fParticles.GetEntriesFast();
771 for (Int_t i = 0; i < np; i++) {
772 TParticle* part = (TParticle*) fParticles.At(i);
773 Int_t kf = part->GetPdgCode();
774 Int_t ksp = part->GetUniqueID();
775 if (kf == 22 && ksp == 40) {
776 Float_t phi = part->Phi();
777 Float_t eta = part->Eta();
778 if (eta < fEtaMaxJet &&
779 eta > fEtaMinJet &&
780 phi < fPhiMaxJet &&
781 phi > fPhiMinJet) {
782 triggered = 1;
783 break;
784 } // check phi,eta within limits
785 } // direct gamma ?
786 } // particle loop
787 } // fTrigger == 2
788 return triggered;
789}