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