]> 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
Moving required CMake version from 2.8.4 to 2.8.8
[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
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};
cdfa7be6 244 Float_t time0 = 0.;
0119ef9a 245 Float_t p[3];
246 Float_t tof;
247
0119ef9a 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];
d12123e8 256 //time0 = fTimeOrigin;
0119ef9a 257
258 if(fVertexSmear == kPerEvent) {
259 Vertex();
260 for (j=0; j < 3; j++)
261 origin0[j] = fVertex[j];
d12123e8 262 //time0 = fTime;
0119ef9a 263 }
264
265 Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
266
267 while(1) {
48beeea0 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
0119ef9a 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;
48beeea0 282
283
0119ef9a 284 fAmpt->ImportParticles(&fParticles,"All");
285 Int_t np = fParticles.GetEntriesFast();
286 if (np == 0 )
287 continue;
774ceaaf 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 //
0119ef9a 299 if (fTrigger != kNoTrigger) {
300 if (!CheckTrigger())
301 continue;
302 }
d60fac1e 303
304 AliDecayer *decayer = 0;
d12123e8 305 //if (gMC)
306 // decayer = gMC->GetDecayer();
307 decayer = fDecayer; //AMPT does not do the strong decays per dafault
d60fac1e 308
309 if (decayer&&fDecay) {
310 TClonesArray arr("TParticle",100);
d12123e8 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);
afe6642c 361 for (Int_t jj = 1; jj < ndecayed; jj++) {
d12123e8 362 TParticle *jp = (TParticle *)arr.At(jj);
363 if (jp->GetFirstMother()!=1)
364 continue;
afe6642c 365
366 TParticle *newp = new(fParticles[np2]) TParticle(jp->GetPdgCode(),
d12123e8 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(),
afe6642c 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() );
d12123e8 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));
d60fac1e 390 }
d12123e8 391 } // end of nLoop->NumberOfNestedLoops
d60fac1e 392 } else {
393 if (fDecay)
394 AliError("No decayer found, but fDecay==kTRUE!");
395 }
396
0119ef9a 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];
d12123e8 414 //fTime = time0;
0119ef9a 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 ?
d12123e8 421 if (Stable(iparticle)) continue; // quit if particle has no daughters
0119ef9a 422 Bool_t selected = kTRUE;
423 Bool_t hasSelectedDaughters = kFALSE;
d12123e8 424 kf = iparticle->GetPdgCode();
425 ks = iparticle->GetStatusCode();
0119ef9a 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 ?
d12123e8 449 if (!Stable(iparticle)) continue; // quit if particle has daughters
450 Bool_t selected = kTRUE;
451 kf = iparticle->GetPdgCode();
d60fac1e 452 if (kf == 92)
453 continue;
d12123e8 454 ks = iparticle->GetStatusCode();
455 ksp = iparticle->GetUniqueID();
0119ef9a 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;
68a6e988 477 if (0) printf("---> %d %d %d %s\n",i,nc,kf,iparticle->GetName());
0119ef9a 478 } // selected
479 } // particle loop final state
480
0119ef9a 481 // Write particles to stack
482 for (Int_t i = 0; i<np; i++) {
0119ef9a 483 if (pSelected[i]) {
d12123e8 484 TParticle *iparticle = (TParticle *) fParticles.At(i);
485 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
486 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
0119ef9a 487 kf = iparticle->GetPdgCode();
488 ks = iparticle->GetStatusCode();
489 p[0] = iparticle->Px();
490 p[1] = iparticle->Py();
491 p[2] = iparticle->Pz() * sign;
774ceaaf 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();
cdfa7be6 496
0119ef9a 497 imo = -1;
498 TParticle* mother = 0;
d12123e8 499 TMCProcess procID = (TMCProcess) iparticle->GetUniqueID();
0119ef9a 500 if (hasMother) {
501 imo = iparticle->GetFirstMother();
502 mother = (TParticle *) fParticles.At(imo);
503 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
d12123e8 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
0119ef9a 509 } // if has mother
510 Bool_t tFlag = (fTrackIt && !hasDaughter);
d12123e8 511 PushTrack(tFlag,imo,kf,p,origin,polar,tof,procID,nt, 1., ks);
0119ef9a 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
d12123e8 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
0119ef9a 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.;
48eb126c 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));
0119ef9a 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
d60fac1e 672Bool_t AliGenAmpt::Stable(TParticle* particle) const
0119ef9a 673{
674 // Return true for a stable particle
675
d60fac1e 676 if (!particle)
677 return kFALSE;
0119ef9a 678 if (particle->GetFirstDaughter() < 0 )
0119ef9a 679 return kTRUE;
d60fac1e 680 return kFALSE;
d12123e8 681
682 /// ADD LIST
683
0119ef9a 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);
48beeea0 701 //fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
702 fHeader->SetReactionPlaneAngle(fAmpt->GetReactionPlaneAngle());
0119ef9a 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);
d12123e8 731 fHeader->SetInteractionTime(fEventTime);
732
0119ef9a 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}