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