]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA6/AliGenPythia.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix for #97283: AliGenPythia : Trigger on Pythia particles for any process
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
CommitLineData
8d2cd130 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
7cdba479 16/* $Id$ */
8d2cd130 17
18//
19// Generator using the TPythia interface (via AliPythia)
20// to generate pp collisions.
21// Using SetNuclei() also nuclear modifications to the structure functions
22// can be taken into account. This makes, of course, only sense for the
23// generation of the products of hard processes (heavy flavor, jets ...)
24//
25// andreas.morsch@cern.ch
26//
27
37b09b91 28#include <TClonesArray.h>
8d2cd130 29#include <TDatabasePDG.h>
30#include <TParticle.h>
31#include <TPDGCode.h>
1058d9df 32#include <TObjArray.h>
8d2cd130 33#include <TSystem.h>
34#include <TTree.h>
8d2cd130 35#include "AliConst.h"
36#include "AliDecayerPythia.h"
37#include "AliGenPythia.h"
cd07c39b 38#include "AliFastGlauber.h"
5fa4b20b 39#include "AliHeader.h"
8d2cd130 40#include "AliGenPythiaEventHeader.h"
41#include "AliPythia.h"
7cdba479 42#include "AliPythiaRndm.h"
8d2cd130 43#include "AliRun.h"
7ea3ea5b 44#include "AliStack.h"
45#include "AliRunLoader.h"
5d12ce38 46#include "AliMC.h"
c93255fe 47#include "AliLog.h"
e2d34d35 48#include "PyquenCommon.h"
800be8b7 49#include "AliLog.h"
8d2cd130 50
014a9521 51ClassImp(AliGenPythia)
8d2cd130 52
e8a8adcd 53
54AliGenPythia::AliGenPythia():
55 AliGenMC(),
56 fProcess(kPyCharm),
efe3b1cd 57 fItune(-1),
e8a8adcd 58 fStrucFunc(kCTEQ5L),
e8a8adcd 59 fKineBias(0.),
60 fTrials(0),
61 fTrialsRun(0),
62 fQ(0.),
63 fX1(0.),
64 fX2(0.),
65 fEventTime(0.),
66 fInteractionRate(0.),
67 fTimeWindow(0.),
68 fCurSubEvent(0),
69 fEventsTime(0),
70 fNev(0),
71 fFlavorSelect(0),
72 fXsection(0.),
73 fPythia(0),
74 fPtHardMin(0.),
75 fPtHardMax(1.e4),
76 fYHardMin(-1.e10),
77 fYHardMax(1.e10),
78 fGinit(1),
79 fGfinal(1),
80 fHadronisation(1),
03358a32 81 fPatchOmegaDalitz(0),
e8a8adcd 82 fNpartons(0),
83 fReadFromFile(0),
84 fQuench(0),
cd07c39b 85 fQhat(0.),
86 fLength(0.),
66b8652c 87 fpyquenT(1.),
88 fpyquenTau(0.1),
89 fpyquenNf(0),
90 fpyquenEloss(0),
91 fpyquenAngle(0),
e6fe9b82 92 fImpact(0.),
e8a8adcd 93 fPtKick(1.),
94 fFullEvent(kTRUE),
95 fDecayer(new AliDecayerPythia()),
96 fDebugEventFirst(-1),
97 fDebugEventLast(-1),
98 fEtMinJet(0.),
99 fEtMaxJet(1.e4),
100 fEtaMinJet(-20.),
101 fEtaMaxJet(20.),
102 fPhiMinJet(0.),
103 fPhiMaxJet(2.* TMath::Pi()),
104 fJetReconstruction(kCell),
105 fEtaMinGamma(-20.),
106 fEtaMaxGamma(20.),
107 fPhiMinGamma(0.),
108 fPhiMaxGamma(2. * TMath::Pi()),
109 fPycellEtaMax(2.),
110 fPycellNEta(274),
111 fPycellNPhi(432),
112 fPycellThreshold(0.),
113 fPycellEtSeed(4.),
114 fPycellMinEtJet(10.),
115 fPycellMaxRadius(1.),
116 fStackFillOpt(kFlavorSelection),
117 fFeedDownOpt(kTRUE),
118 fFragmentation(kTRUE),
119 fSetNuclei(kFALSE),
120 fNewMIS(kFALSE),
121 fHFoff(kFALSE),
20e47f08 122 fNucPdf(0),
e8a8adcd 123 fTriggerParticle(0),
124 fTriggerEta(0.9),
2591bd0e 125 fTriggerMinPt(-1),
126 fTriggerMaxPt(1000),
700b9416 127 fTriggerMultiplicity(0),
128 fTriggerMultiplicityEta(0),
38112f3f 129 fTriggerMultiplicityPtMin(0),
e8a8adcd 130 fCountMode(kCountAll),
131 fHeader(0),
132 fRL(0),
777e69b0 133 fkFileName(0),
9fd8e520 134 fFragPhotonInCalo(kFALSE),
43e3b80a 135 fHadronInCalo(kFALSE) ,
ec2c406e 136 fPi0InCalo(kFALSE) ,
90a236ce 137 fPhotonInCalo(kFALSE),
9dfe63b3 138 fEleInEMCAL(kFALSE),
9fd8e520 139 fCheckEMCAL(kFALSE),
140 fCheckPHOS(kFALSE),
90a236ce 141 fCheckPHOSeta(kFALSE),
43e3b80a 142 fTriggerParticleMinPt(0),
90a236ce 143 fPhotonMinPt(0),
9dfe63b3 144 fElectronMinPt(0),
9fd8e520 145 fPHOSMinPhi(219.),
146 fPHOSMaxPhi(321.),
147 fPHOSEta(0.13),
148 fEMCALMinPhi(79.),
149 fEMCALMaxPhi(191.),
800be8b7 150 fEMCALEta(0.71),
151 fkTuneForDiff(0),
152 fProcDiff(0)
8d2cd130 153{
154// Default Constructor
e7c989e4 155 fEnergyCMS = 5500.;
7cdba479 156 if (!AliPythiaRndm::GetPythiaRandom())
e8a8adcd 157 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 158}
159
160AliGenPythia::AliGenPythia(Int_t npart)
e8a8adcd 161 :AliGenMC(npart),
162 fProcess(kPyCharm),
efe3b1cd 163 fItune(-1),
e8a8adcd 164 fStrucFunc(kCTEQ5L),
e8a8adcd 165 fKineBias(0.),
166 fTrials(0),
167 fTrialsRun(0),
168 fQ(0.),
169 fX1(0.),
170 fX2(0.),
171 fEventTime(0.),
172 fInteractionRate(0.),
173 fTimeWindow(0.),
174 fCurSubEvent(0),
175 fEventsTime(0),
176 fNev(0),
177 fFlavorSelect(0),
178 fXsection(0.),
179 fPythia(0),
180 fPtHardMin(0.),
181 fPtHardMax(1.e4),
182 fYHardMin(-1.e10),
183 fYHardMax(1.e10),
184 fGinit(kTRUE),
185 fGfinal(kTRUE),
186 fHadronisation(kTRUE),
03358a32 187 fPatchOmegaDalitz(0),
e8a8adcd 188 fNpartons(0),
189 fReadFromFile(kFALSE),
190 fQuench(kFALSE),
cd07c39b 191 fQhat(0.),
192 fLength(0.),
66b8652c 193 fpyquenT(1.),
194 fpyquenTau(0.1),
195 fpyquenNf(0),
196 fpyquenEloss(0),
197 fpyquenAngle(0),
e6fe9b82 198 fImpact(0.),
e8a8adcd 199 fPtKick(1.),
200 fFullEvent(kTRUE),
201 fDecayer(new AliDecayerPythia()),
202 fDebugEventFirst(-1),
203 fDebugEventLast(-1),
204 fEtMinJet(0.),
205 fEtMaxJet(1.e4),
206 fEtaMinJet(-20.),
207 fEtaMaxJet(20.),
208 fPhiMinJet(0.),
209 fPhiMaxJet(2.* TMath::Pi()),
210 fJetReconstruction(kCell),
211 fEtaMinGamma(-20.),
212 fEtaMaxGamma(20.),
213 fPhiMinGamma(0.),
214 fPhiMaxGamma(2. * TMath::Pi()),
215 fPycellEtaMax(2.),
216 fPycellNEta(274),
217 fPycellNPhi(432),
218 fPycellThreshold(0.),
219 fPycellEtSeed(4.),
220 fPycellMinEtJet(10.),
221 fPycellMaxRadius(1.),
222 fStackFillOpt(kFlavorSelection),
223 fFeedDownOpt(kTRUE),
224 fFragmentation(kTRUE),
225 fSetNuclei(kFALSE),
226 fNewMIS(kFALSE),
227 fHFoff(kFALSE),
20e47f08 228 fNucPdf(0),
e8a8adcd 229 fTriggerParticle(0),
2591bd0e 230 fTriggerEta(0.9),
231 fTriggerMinPt(-1),
232 fTriggerMaxPt(1000),
700b9416 233 fTriggerMultiplicity(0),
234 fTriggerMultiplicityEta(0),
38112f3f 235 fTriggerMultiplicityPtMin(0),
e8a8adcd 236 fCountMode(kCountAll),
237 fHeader(0),
238 fRL(0),
777e69b0 239 fkFileName(0),
9fd8e520 240 fFragPhotonInCalo(kFALSE),
43e3b80a 241 fHadronInCalo(kFALSE) ,
ec2c406e 242 fPi0InCalo(kFALSE) ,
90a236ce 243 fPhotonInCalo(kFALSE),
9dfe63b3 244 fEleInEMCAL(kFALSE),
9fd8e520 245 fCheckEMCAL(kFALSE),
246 fCheckPHOS(kFALSE),
90a236ce 247 fCheckPHOSeta(kFALSE),
43e3b80a 248 fTriggerParticleMinPt(0),
90a236ce 249 fPhotonMinPt(0),
9dfe63b3 250 fElectronMinPt(0),
9fd8e520 251 fPHOSMinPhi(219.),
252 fPHOSMaxPhi(321.),
253 fPHOSEta(0.13),
254 fEMCALMinPhi(79.),
255 fEMCALMaxPhi(191.),
800be8b7 256 fEMCALEta(0.71),
257 fkTuneForDiff(0),
258 fProcDiff(0)
8d2cd130 259{
260// default charm production at 5. 5 TeV
261// semimuonic decay
262// structure function GRVHO
263//
e7c989e4 264 fEnergyCMS = 5500.;
8d2cd130 265 fName = "Pythia";
266 fTitle= "Particle Generator using PYTHIA";
8d2cd130 267 SetForceDecay();
8d2cd130 268 // Set random number generator
7cdba479 269 if (!AliPythiaRndm::GetPythiaRandom())
270 AliPythiaRndm::SetPythiaRandom(GetRandom());
76d6ba9a 271 }
8d2cd130 272
8d2cd130 273AliGenPythia::~AliGenPythia()
274{
275// Destructor
9d7108a7 276 if(fEventsTime) delete fEventsTime;
277}
278
279void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
280{
281// Generate pileup using user specified rate
282 fInteractionRate = rate;
283 fTimeWindow = timewindow;
284 GeneratePileup();
285}
286
287void AliGenPythia::GeneratePileup()
288{
289// Generate sub events time for pileup
290 fEventsTime = 0;
291 if(fInteractionRate == 0.) {
292 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
293 return;
294 }
295
296 Int_t npart = NumberParticles();
297 if(npart < 0) {
298 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
299 return;
300 }
301
302 if(fEventsTime) delete fEventsTime;
303 fEventsTime = new TArrayF(npart);
304 TArrayF &array = *fEventsTime;
305 for(Int_t ipart = 0; ipart < npart; ipart++)
306 array[ipart] = 0.;
307
308 Float_t eventtime = 0.;
309 while(1)
310 {
311 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
312 if(eventtime > fTimeWindow) break;
313 array.Set(array.GetSize()+1);
314 array[array.GetSize()-1] = eventtime;
315 }
316
317 eventtime = 0.;
318 while(1)
319 {
320 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
321 if(TMath::Abs(eventtime) > fTimeWindow) break;
322 array.Set(array.GetSize()+1);
323 array[array.GetSize()-1] = eventtime;
324 }
325
326 SetNumberParticles(fEventsTime->GetSize());
8d2cd130 327}
328
592f8307 329void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
330 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
331{
332// Set pycell parameters
333 fPycellEtaMax = etamax;
334 fPycellNEta = neta;
335 fPycellNPhi = nphi;
336 fPycellThreshold = thresh;
337 fPycellEtSeed = etseed;
338 fPycellMinEtJet = minet;
339 fPycellMaxRadius = r;
340}
341
342
343
8d2cd130 344void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
345{
346 // Set a range of event numbers, for which a table
347 // of generated particle will be printed
348 fDebugEventFirst = eventFirst;
349 fDebugEventLast = eventLast;
350 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
351}
352
353void AliGenPythia::Init()
354{
355// Initialisation
356
357 SetMC(AliPythia::Instance());
b88f5cea 358 fPythia=(AliPythia*) fMCEvGen;
e2bddf81 359
8d2cd130 360//
361 fParentWeight=1./Float_t(fNpart);
362//
8d2cd130 363
364
365 fPythia->SetCKIN(3,fPtHardMin);
366 fPythia->SetCKIN(4,fPtHardMax);
367 fPythia->SetCKIN(7,fYHardMin);
368 fPythia->SetCKIN(8,fYHardMax);
369
20e47f08 370 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
8d2cd130 371 // Fragmentation?
372 if (fFragmentation) {
373 fPythia->SetMSTP(111,1);
374 } else {
375 fPythia->SetMSTP(111,0);
376 }
377
378
379// initial state radiation
380 fPythia->SetMSTP(61,fGinit);
381// final state radiation
382 fPythia->SetMSTP(71,fGfinal);
383// pt - kick
384 if (fPtKick > 0.) {
385 fPythia->SetMSTP(91,1);
688950ef 386 fPythia->SetPARP(91,fPtKick);
387 fPythia->SetPARP(93, 4. * fPtKick);
8d2cd130 388 } else {
389 fPythia->SetMSTP(91,0);
390 }
391
5fa4b20b 392
393 if (fReadFromFile) {
777e69b0 394 fRL = AliRunLoader::Open(fkFileName, "Partons");
5fa4b20b 395 fRL->LoadKinematics();
396 fRL->LoadHeader();
397 } else {
398 fRL = 0x0;
399 }
fdea4387 400 //
efe3b1cd 401 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
fdea4387 402 // Forward Paramters to the AliPythia object
403 fDecayer->SetForceDecay(fForceDecay);
beac474c 404// Switch off Heavy Flavors on request
405 if (fHFoff) {
51387949 406 // Maximum number of quark flavours used in pdf
beac474c 407 fPythia->SetMSTP(58, 3);
51387949 408 // Maximum number of flavors that can be used in showers
beac474c 409 fPythia->SetMSTJ(45, 3);
51387949 410 // Switch off g->QQbar splitting in decay table
411 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
beac474c 412 }
fdea4387 413
51387949 414 fDecayer->Init();
415
8d2cd130 416
417// Parent and Children Selection
418 switch (fProcess)
419 {
65f2626c 420 case kPyOldUEQ2ordered:
421 case kPyOldUEQ2ordered2:
422 case kPyOldPopcorn:
423 break;
8d2cd130 424 case kPyCharm:
425 case kPyCharmUnforced:
adf4d898 426 case kPyCharmPbPbMNR:
aabc7187 427 case kPyCharmpPbMNR:
e0e89f40 428 case kPyCharmppMNR:
429 case kPyCharmppMNRwmi:
8d2cd130 430 fParentSelect[0] = 411;
431 fParentSelect[1] = 421;
432 fParentSelect[2] = 431;
433 fParentSelect[3] = 4122;
9ccc3d0e 434 fParentSelect[4] = 4232;
435 fParentSelect[5] = 4132;
436 fParentSelect[6] = 4332;
8d2cd130 437 fFlavorSelect = 4;
438 break;
adf4d898 439 case kPyD0PbPbMNR:
440 case kPyD0pPbMNR:
441 case kPyD0ppMNR:
8d2cd130 442 fParentSelect[0] = 421;
443 fFlavorSelect = 4;
444 break;
90d7b703 445 case kPyDPlusPbPbMNR:
446 case kPyDPluspPbMNR:
447 case kPyDPlusppMNR:
448 fParentSelect[0] = 411;
449 fFlavorSelect = 4;
450 break;
e0e89f40 451 case kPyDPlusStrangePbPbMNR:
452 case kPyDPlusStrangepPbMNR:
453 case kPyDPlusStrangeppMNR:
454 fParentSelect[0] = 431;
455 fFlavorSelect = 4;
456 break;
68504d42 457 case kPyLambdacppMNR:
458 fParentSelect[0] = 4122;
459 fFlavorSelect = 4;
460 break;
8d2cd130 461 case kPyBeauty:
9dfe63b3 462 case kPyBeautyJets:
adf4d898 463 case kPyBeautyPbPbMNR:
464 case kPyBeautypPbMNR:
465 case kPyBeautyppMNR:
e0e89f40 466 case kPyBeautyppMNRwmi:
8d2cd130 467 fParentSelect[0]= 511;
468 fParentSelect[1]= 521;
469 fParentSelect[2]= 531;
470 fParentSelect[3]= 5122;
471 fParentSelect[4]= 5132;
472 fParentSelect[5]= 5232;
473 fParentSelect[6]= 5332;
474 fFlavorSelect = 5;
475 break;
476 case kPyBeautyUnforced:
477 fParentSelect[0] = 511;
478 fParentSelect[1] = 521;
479 fParentSelect[2] = 531;
480 fParentSelect[3] = 5122;
481 fParentSelect[4] = 5132;
482 fParentSelect[5] = 5232;
483 fParentSelect[6] = 5332;
484 fFlavorSelect = 5;
485 break;
486 case kPyJpsiChi:
487 case kPyJpsi:
488 fParentSelect[0] = 443;
489 break;
f529e69b 490 case kPyMbDefault:
0bd3d7c5 491 case kPyMbAtlasTuneMC09:
8d2cd130 492 case kPyMb:
04081a91 493 case kPyMbWithDirectPhoton:
8d2cd130 494 case kPyMbNonDiffr:
d7de4a67 495 case kPyMbMSEL1:
8d2cd130 496 case kPyJets:
497 case kPyDirectGamma:
e33acb67 498 case kPyLhwgMb:
8d2cd130 499 break;
589380c6 500 case kPyW:
0f6ee828 501 case kPyZ:
589380c6 502 break;
8d2cd130 503 }
504//
592f8307 505//
506// JetFinder for Trigger
507//
508// Configure detector (EMCAL like)
509//
d7de4a67 510 fPythia->SetPARU(51, fPycellEtaMax);
511 fPythia->SetMSTU(51, fPycellNEta);
512 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 513//
514// Configure Jet Finder
515//
d7de4a67 516 fPythia->SetPARU(58, fPycellThreshold);
517 fPythia->SetPARU(52, fPycellEtSeed);
518 fPythia->SetPARU(53, fPycellMinEtJet);
519 fPythia->SetPARU(54, fPycellMaxRadius);
520 fPythia->SetMSTU(54, 2);
592f8307 521//
8d2cd130 522// This counts the total number of calls to Pyevnt() per run.
523 fTrialsRun = 0;
524 fQ = 0.;
525 fX1 = 0.;
526 fX2 = 0.;
527 fNev = 0 ;
528//
1d568bc2 529//
530//
8d2cd130 531 AliGenMC::Init();
1d568bc2 532//
533//
534//
535 if (fSetNuclei) {
536 fDyBoost = 0;
537 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
538 }
d7de4a67 539
cd07c39b 540 fPythia->SetPARJ(200, 0.0);
541 fPythia->SetPARJ(199, 0.0);
542 fPythia->SetPARJ(198, 0.0);
543 fPythia->SetPARJ(197, 0.0);
544
545 if (fQuench == 1) {
5fa4b20b 546 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 547 }
3a709cfa 548
66b8652c 549 if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
550
b976f7d8 551 if (fQuench == 3) {
552 // Nestor's change of the splittings
553 fPythia->SetPARJ(200, 0.8);
554 fPythia->SetMSTJ(41, 1); // QCD radiation only
555 fPythia->SetMSTJ(42, 2); // angular ordering
556 fPythia->SetMSTJ(44, 2); // option to run alpha_s
557 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
558 fPythia->SetMSTJ(50, 0); // No coherence in first branching
559 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
cd07c39b 560 } else if (fQuench == 4) {
561 // Armesto-Cunqueiro-Salgado change of the splittings.
562 AliFastGlauber* glauber = AliFastGlauber::Instance();
563 glauber->Init(2);
564 //read and store transverse almonds corresponding to differnt
565 //impact parameters.
cd07c39b 566 glauber->SetCentralityClass(0.,0.1);
567 fPythia->SetPARJ(200, 1.);
568 fPythia->SetPARJ(198, fQhat);
569 fPythia->SetPARJ(199, fLength);
cd07c39b 570 fPythia->SetMSTJ(42, 2); // angular ordering
571 fPythia->SetMSTJ(44, 2); // option to run alpha_s
cd07c39b 572 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
b976f7d8 573 }
8d2cd130 574}
575
576void AliGenPythia::Generate()
577{
578// Generate one event
19ef8cf0 579 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
8d2cd130 580 fDecayer->ForceDecay();
581
13cca24a 582 Double_t polar[3] = {0,0,0};
583 Double_t origin[3] = {0,0,0};
584 Double_t p[4];
8d2cd130 585// converts from mm/c to s
586 const Float_t kconv=0.001/2.999792458e8;
587//
588 Int_t nt=0;
589 Int_t jev=0;
590 Int_t j, kf;
591 fTrials=0;
f913ec4f 592 fEventTime = 0.;
593
2590ccf9 594
8d2cd130 595
596 // Set collision vertex position
2590ccf9 597 if (fVertexSmear == kPerEvent) Vertex();
598
8d2cd130 599// event loop
600 while(1)
601 {
32d6ef7d 602//
5fa4b20b 603// Produce event
32d6ef7d 604//
32d6ef7d 605//
5fa4b20b 606// Switch hadronisation off
607//
608 fPythia->SetMSTJ(1, 0);
cd07c39b 609
610 if (fQuench ==4){
611 Double_t bimp;
612 // Quenching comes through medium-modified splitting functions.
613 AliFastGlauber::Instance()->GetRandomBHard(bimp);
e6fe9b82 614 fPythia->SetPARJ(197, bimp);
615 fImpact = bimp;
6c43eccb 616 fPythia->Qpygin0();
cd07c39b 617 }
32d6ef7d 618//
5fa4b20b 619// Either produce new event or read partons from file
620//
621 if (!fReadFromFile) {
beac474c 622 if (!fNewMIS) {
623 fPythia->Pyevnt();
624 } else {
625 fPythia->Pyevnw();
626 }
5fa4b20b 627 fNpartons = fPythia->GetN();
628 } else {
33c3c91a 629 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
630 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
5fa4b20b 631 fPythia->SetN(0);
632 LoadEvent(fRL->Stack(), 0 , 1);
633 fPythia->Pyedit(21);
634 }
635
32d6ef7d 636//
637// Run quenching routine
638//
5fa4b20b 639 if (fQuench == 1) {
640 fPythia->Quench();
641 } else if (fQuench == 2){
642 fPythia->Pyquen(208., 0, 0.);
b976f7d8 643 } else if (fQuench == 3) {
644 // Quenching is via multiplicative correction of the splittings
5fa4b20b 645 }
b976f7d8 646
32d6ef7d 647//
5fa4b20b 648// Switch hadronisation on
32d6ef7d 649//
20e47f08 650 if (fHadronisation) {
651 fPythia->SetMSTJ(1, 1);
5fa4b20b 652//
653// .. and perform hadronisation
aea21c57 654// printf("Calling hadronisation %d\n", fPythia->GetN());
03358a32 655
656 if (fPatchOmegaDalitz) {
657 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
658 fPythia->Pyexec();
659 fPythia->DalitzDecays();
660 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
661 }
662 fPythia->Pyexec();
20e47f08 663 }
8d2cd130 664 fTrials++;
8507138f 665 fPythia->ImportParticles(&fParticles,"All");
03358a32 666
df275cfa 667 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
e1cf8975 668 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
f64640b5 669
8d2cd130 670//
671//
672//
673 Int_t i;
674
dbd64db6 675 fNprimaries = 0;
8507138f 676 Int_t np = fParticles.GetEntriesFast();
5fa4b20b 677
7c21f297 678 if (np == 0) continue;
8d2cd130 679//
2590ccf9 680
8d2cd130 681//
682 Int_t* pParent = new Int_t[np];
683 Int_t* pSelected = new Int_t[np];
684 Int_t* trackIt = new Int_t[np];
5fa4b20b 685 for (i = 0; i < np; i++) {
8d2cd130 686 pParent[i] = -1;
687 pSelected[i] = 0;
688 trackIt[i] = 0;
689 }
690
691 Int_t nc = 0; // Total n. of selected particles
692 Int_t nParents = 0; // Selected parents
693 Int_t nTkbles = 0; // Trackable particles
f529e69b 694 if (fProcess != kPyMbDefault &&
695 fProcess != kPyMb &&
6d2ec66d 696 fProcess != kPyMbAtlasTuneMC09 &&
04081a91 697 fProcess != kPyMbWithDirectPhoton &&
f529e69b 698 fProcess != kPyJets &&
8d2cd130 699 fProcess != kPyDirectGamma &&
589380c6 700 fProcess != kPyMbNonDiffr &&
d7de4a67 701 fProcess != kPyMbMSEL1 &&
f529e69b 702 fProcess != kPyW &&
703 fProcess != kPyZ &&
704 fProcess != kPyCharmppMNRwmi &&
9dfe63b3 705 fProcess != kPyBeautyppMNRwmi &&
706 fProcess != kPyBeautyJets) {
8d2cd130 707
5fa4b20b 708 for (i = 0; i < np; i++) {
8507138f 709 TParticle* iparticle = (TParticle *) fParticles.At(i);
8d2cd130 710 Int_t ks = iparticle->GetStatusCode();
711 kf = CheckPDGCode(iparticle->GetPdgCode());
712// No initial state partons
713 if (ks==21) continue;
714//
715// Heavy Flavor Selection
716//
717 // quark ?
718 kf = TMath::Abs(kf);
719 Int_t kfl = kf;
9ff6c04c 720 // Resonance
f913ec4f 721
9ff6c04c 722 if (kfl > 100000) kfl %= 100000;
183a5ca9 723 if (kfl > 10000) kfl %= 10000;
8d2cd130 724 // meson ?
725 if (kfl > 10) kfl/=100;
726 // baryon
727 if (kfl > 10) kfl/=10;
8d2cd130 728 Int_t ipa = iparticle->GetFirstMother()-1;
729 Int_t kfMo = 0;
f913ec4f 730//
731// Establish mother daughter relation between heavy quarks and mesons
732//
733 if (kf >= fFlavorSelect && kf <= 6) {
734 Int_t idau = iparticle->GetFirstDaughter() - 1;
735 if (idau > -1) {
8507138f 736 TParticle* daughter = (TParticle *) fParticles.At(idau);
f913ec4f 737 Int_t pdgD = daughter->GetPdgCode();
738 if (pdgD == 91 || pdgD == 92) {
739 Int_t jmin = daughter->GetFirstDaughter() - 1;
740 Int_t jmax = daughter->GetLastDaughter() - 1;
2ab330c9 741 for (Int_t jp = jmin; jp <= jmax; jp++)
742 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
f913ec4f 743 } // is string or cluster
744 } // has daughter
745 } // heavy quark
8d2cd130 746
f913ec4f 747
8d2cd130 748 if (ipa > -1) {
8507138f 749 TParticle * mother = (TParticle *) fParticles.At(ipa);
8d2cd130 750 kfMo = TMath::Abs(mother->GetPdgCode());
751 }
f913ec4f 752
8d2cd130 753 // What to keep in Stack?
754 Bool_t flavorOK = kFALSE;
755 Bool_t selectOK = kFALSE;
756 if (fFeedDownOpt) {
32d6ef7d 757 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 758 } else {
32d6ef7d 759 if (kfl > fFlavorSelect) {
760 nc = -1;
761 break;
762 }
763 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 764 }
765 switch (fStackFillOpt) {
06956108 766 case kHeavyFlavor:
8d2cd130 767 case kFlavorSelection:
32d6ef7d 768 selectOK = kTRUE;
769 break;
8d2cd130 770 case kParentSelection:
32d6ef7d 771 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
772 break;
8d2cd130 773 }
774 if (flavorOK && selectOK) {
775//
776// Heavy flavor hadron or quark
777//
778// Kinematic seletion on final state heavy flavor mesons
779 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
780 {
9ff6c04c 781 continue;
8d2cd130 782 }
783 pSelected[i] = 1;
784 if (ParentSelected(kf)) ++nParents; // Update parent count
785// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
786 } else {
787// Kinematic seletion on decay products
788 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 789 && !KinematicSelection(iparticle, 1))
8d2cd130 790 {
9ff6c04c 791 continue;
8d2cd130 792 }
793//
794// Decay products
795// Select if mother was selected and is not tracked
796
797 if (pSelected[ipa] &&
798 !trackIt[ipa] && // mother will be tracked ?
799 kfMo != 5 && // mother is b-quark, don't store fragments
800 kfMo != 4 && // mother is c-quark, don't store fragments
801 kf != 92) // don't store string
802 {
803//
804// Semi-stable or de-selected: diselect decay products:
805//
806//
807 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
808 {
809 Int_t ipF = iparticle->GetFirstDaughter();
810 Int_t ipL = iparticle->GetLastDaughter();
811 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
812 }
813// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
814 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
815 }
816 }
817 if (pSelected[i] == -1) pSelected[i] = 0;
818 if (!pSelected[i]) continue;
819 // Count quarks only if you did not include fragmentation
820 if (fFragmentation && kf <= 10) continue;
9ff6c04c 821
8d2cd130 822 nc++;
823// Decision on tracking
824 trackIt[i] = 0;
825//
826// Track final state particle
827 if (ks == 1) trackIt[i] = 1;
828// Track semi-stable particles
d25cfd65 829 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 830// Track particles selected by process if undecayed.
831 if (fForceDecay == kNoDecay) {
832 if (ParentSelected(kf)) trackIt[i] = 1;
833 } else {
834 if (ParentSelected(kf)) trackIt[i] = 0;
835 }
836 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
837//
838//
839
840 } // particle selection loop
841 if (nc > 0) {
842 for (i = 0; i<np; i++) {
843 if (!pSelected[i]) continue;
8507138f 844 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 845 kf = CheckPDGCode(iparticle->GetPdgCode());
846 Int_t ks = iparticle->GetStatusCode();
847 p[0] = iparticle->Px();
848 p[1] = iparticle->Py();
849 p[2] = iparticle->Pz();
a920faf9 850 p[3] = iparticle->Energy();
851
2590ccf9 852 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
853 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
854 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
855
21391258 856 Float_t tof = fTime + kconv*iparticle->T();
8d2cd130 857 Int_t ipa = iparticle->GetFirstMother()-1;
858 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 859
860 PushTrack(fTrackIt*trackIt[i], iparent, kf,
861 p[0], p[1], p[2], p[3],
862 origin[0], origin[1], origin[2], tof,
863 polar[0], polar[1], polar[2],
864 kPPrimary, nt, 1., ks);
8d2cd130 865 pParent[i] = nt;
dbd64db6 866 KeepTrack(nt);
867 fNprimaries++;
642f15cf 868 } // PushTrack loop
8d2cd130 869 }
870 } else {
871 nc = GenerateMB();
872 } // mb ?
f913ec4f 873
874 GetSubEventTime();
8d2cd130 875
234f6d32 876 delete[] pParent;
877 delete[] pSelected;
878 delete[] trackIt;
8d2cd130 879
880 if (nc > 0) {
881 switch (fCountMode) {
882 case kCountAll:
883 // printf(" Count all \n");
884 jev += nc;
885 break;
886 case kCountParents:
887 // printf(" Count parents \n");
888 jev += nParents;
889 break;
890 case kCountTrackables:
891 // printf(" Count trackable \n");
892 jev += nTkbles;
893 break;
894 }
895 if (jev >= fNpart || fNpart == -1) {
896 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 897
8d2cd130 898 fQ += fPythia->GetVINT(51);
899 fX1 += fPythia->GetVINT(41);
900 fX2 += fPythia->GetVINT(42);
901 fTrialsRun += fTrials;
902 fNev++;
903 MakeHeader();
904 break;
905 }
906 }
907 } // event loop
908 SetHighWaterMark(nt);
909// adjust weight due to kinematic selection
b88f5cea 910// AdjustWeights();
8d2cd130 911// get cross-section
912 fXsection=fPythia->GetPARI(1);
913}
914
915Int_t AliGenPythia::GenerateMB()
916{
917//
918// Min Bias selection and other global selections
919//
920 Int_t i, kf, nt, iparent;
921 Int_t nc = 0;
13cca24a 922 Double_t p[4];
923 Double_t polar[3] = {0,0,0};
924 Double_t origin[3] = {0,0,0};
8d2cd130 925// converts from mm/c to s
926 const Float_t kconv=0.001/2.999792458e8;
927
e0e89f40 928
929
8507138f 930 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
aea21c57 931
5fa4b20b 932
e0e89f40 933
8d2cd130 934 Int_t* pParent = new Int_t[np];
935 for (i=0; i< np; i++) pParent[i] = -1;
2f405d65 936 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8507138f 937 TParticle* jet1 = (TParticle *) fParticles.At(6);
938 TParticle* jet2 = (TParticle *) fParticles.At(7);
234f6d32 939 if (!CheckTrigger(jet1, jet2)) {
940 delete [] pParent;
941 return 0;
942 }
8d2cd130 943 }
e0e89f40 944
ac01c7c6 945 // Select events with fragmentation photon or pi0 or hadrons going to PHOS or EMCAL,
946 // implemented primaryly for kPyJets, but extended to any kind of process.
947 if ((fFragPhotonInCalo || fPi0InCalo || fHadronInCalo) &&
43e3b80a 948 (fCheckPHOS || fCheckEMCAL) ) {
949
950 Bool_t ok = kFALSE;
951
952 for (i=0; i< np; i++) {
953
954 TParticle* iparticle = (TParticle *) fParticles.At(i);
955
956 Int_t pdg = iparticle->GetPdgCode();
957 Int_t status = iparticle->GetStatusCode();
958 ok = kFALSE;
959
960 if (fFragPhotonInCalo && pdg == 22 && status == 1)
961 {
962 Int_t imother = iparticle->GetFirstMother() - 1;
963 TParticle* pmother = (TParticle *) fParticles.At(imother);
964
965 if(pmother->GetStatusCode() != 11) ok = kTRUE ; // No photon from hadron decay
ec2c406e 966 }
43e3b80a 967 else if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
968 {
969 //printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
970 ok = kTRUE;
971 }
972 else if (fHadronInCalo && status == 1)
973 {
974 if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
975 // (in case neutral mesons were declared stable)
976 ok = kTRUE;
1ee02229 977 }
43e3b80a 978
979 if(ok && iparticle->Pt() > fTriggerParticleMinPt)
980 {
981 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
982 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
983
984 if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
985 (fCheckPHOS && IsInPHOS (phi,eta)) )
986 {
987 ok =kTRUE;
988 AliDebug(1,Form("Selected trigger pdg %d, status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
989 break;
990 }
991 }
992 else ok = kFALSE;
993 }
994
995 if(!ok) {
996 delete[] pParent;
997 return 0;
ec2c406e 998 }
43e3b80a 999 }
9dfe63b3 1000
43e3b80a 1001 // Select beauty jets with electron in EMCAL
9dfe63b3 1002 if (fProcess == kPyBeautyJets && fEleInEMCAL) {
1003
1004 Bool_t ok = kFALSE;
1005
1006 Int_t pdg = 11; //electron
1007
70574ff8 1008 Float_t pt = 0.;
1009 Float_t eta = 0.;
1010 Float_t phi = 0.;
9dfe63b3 1011 for (i=0; i< np; i++) {
1012 TParticle* iparticle = (TParticle *) fParticles.At(i);
1013 if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
1014 iparticle->Pt() > fElectronMinPt){
1015 pt = iparticle->Pt();
1016 phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1017 eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
1018 if(IsInEMCAL(phi,eta))
1019 ok =kTRUE;
1020 }
1021 }
92847124 1022 if(!ok) {
1023 delete[] pParent;
9dfe63b3 1024 return 0;
92847124 1025 }
1026
9dfe63b3 1027 AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
1028 }
800be8b7 1029 // Check for diffraction
1030 if(fkTuneForDiff) {
c5dfa3e4 1031 if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
800be8b7 1032 if(!CheckDiffraction()) {
1033 delete [] pParent;
1034 return 0;
1035 }
1036 }
1037 }
1038
700b9416 1039 // Check for minimum multiplicity
1040 if (fTriggerMultiplicity > 0) {
1041 Int_t multiplicity = 0;
1042 for (i = 0; i < np; i++) {
8507138f 1043 TParticle * iparticle = (TParticle *) fParticles.At(i);
700b9416 1044
1045 Int_t statusCode = iparticle->GetStatusCode();
1046
1047 // Initial state particle
e296848e 1048 if (statusCode != 1)
700b9416 1049 continue;
38112f3f 1050 // eta cut
700b9416 1051 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1052 continue;
38112f3f 1053 // pt cut
1054 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1055 continue;
1056
700b9416 1057 TParticlePDG* pdgPart = iparticle->GetPDG();
1058 if (pdgPart && pdgPart->Charge() == 0)
1059 continue;
1060
1061 ++multiplicity;
1062 }
1063
1064 if (multiplicity < fTriggerMultiplicity) {
1065 delete [] pParent;
1066 return 0;
1067 }
38112f3f 1068 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
700b9416 1069 }
90a236ce 1070
1071 // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
1072 if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
1073
1074 Bool_t okd = kFALSE;
1075
1076 Int_t pdg = 22;
1077 Int_t iphcand = -1;
1078 for (i=0; i< np; i++) {
8507138f 1079 TParticle* iparticle = (TParticle *) fParticles.At(i);
90a236ce 1080 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1081 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
1082
1083 if(iparticle->GetStatusCode() == 1
1084 && iparticle->GetPdgCode() == pdg
1085 && iparticle->Pt() > fPhotonMinPt
04a26b0a 1086 && eta < fPHOSEta){
90a236ce 1087
1088 // first check if the photon is in PHOS phi
1089 if(IsInPHOS(phi,eta)){
1090 okd = kTRUE;
1091 break;
1092 }
1093 if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
1094
1095 }
1096 }
1097
1098 if(!okd && iphcand != -1) // execute rotation in phi
1099 RotatePhi(iphcand,okd);
1100
b3305127 1101 if(!okd) {
1102 delete [] pParent;
1103 return 0;
1104 }
90a236ce 1105 }
1106
7ce1321b 1107 if (fTriggerParticle) {
1108 Bool_t triggered = kFALSE;
1109 for (i = 0; i < np; i++) {
8507138f 1110 TParticle * iparticle = (TParticle *) fParticles.At(i);
7ce1321b 1111 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 1112 if (kf != fTriggerParticle) continue;
7ce1321b 1113 if (iparticle->Pt() == 0.) continue;
1114 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
2591bd0e 1115 if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
7ce1321b 1116 triggered = kTRUE;
1117 break;
1118 }
234f6d32 1119 if (!triggered) {
1120 delete [] pParent;
1121 return 0;
1122 }
7ce1321b 1123 }
e0e89f40 1124
1125
1126 // Check if there is a ccbar or bbbar pair with at least one of the two
1127 // in fYMin < y < fYMax
2f405d65 1128
9dfe63b3 1129 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
9ccc3d0e 1130 TParticle *partCheck;
1131 TParticle *mother;
e0e89f40 1132 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 1133 Bool_t theChild=kFALSE;
5d76923e 1134 Bool_t triggered=kFALSE;
9ccc3d0e 1135 Float_t y;
1136 Int_t pdg,mpdg,mpdgUpperFamily;
e0e89f40 1137 for(i=0; i<np; i++) {
9ccc3d0e 1138 partCheck = (TParticle*)fParticles.At(i);
1139 pdg = partCheck->GetPdgCode();
1140 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1141 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1142 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1143 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1144 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1145 }
5d76923e 1146 if(fTriggerParticle) {
1147 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1148 }
9ccc3d0e 1149 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1150 Int_t mi = partCheck->GetFirstMother() - 1;
1151 if(mi<0) continue;
1152 mother = (TParticle*)fParticles.At(mi);
1153 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 1154 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 1155 if ( ParentSelected(mpdg) ||
1156 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1157 if (KinematicSelection(partCheck,1)) {
1158 theChild=kTRUE;
1159 }
1160 }
1161 }
e0e89f40 1162 }
9ccc3d0e 1163 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
234f6d32 1164 delete[] pParent;
e0e89f40 1165 return 0;
1166 }
9ccc3d0e 1167 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1168 delete[] pParent;
1169 return 0;
1170 }
5d76923e 1171 if(fTriggerParticle && !triggered) { // particle requested is not produced
1172 delete[] pParent;
1173 return 0;
1174 }
9ccc3d0e 1175
e0e89f40 1176 }
aea21c57 1177
0f6ee828 1178 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
f529e69b 1179 if ( (fProcess == kPyW ||
1180 fProcess == kPyZ ||
1181 fProcess == kPyMbDefault ||
1182 fProcess == kPyMb ||
6d2ec66d 1183 fProcess == kPyMbAtlasTuneMC09 ||
04081a91 1184 fProcess == kPyMbWithDirectPhoton ||
f529e69b 1185 fProcess == kPyMbNonDiffr)
0f6ee828 1186 && (fCutOnChild == 1) ) {
1187 if ( !CheckKinematicsOnChild() ) {
234f6d32 1188 delete[] pParent;
0f6ee828 1189 return 0;
1190 }
aea21c57 1191 }
1192
1193
f913ec4f 1194 for (i = 0; i < np; i++) {
8d2cd130 1195 Int_t trackIt = 0;
8507138f 1196 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 1197 kf = CheckPDGCode(iparticle->GetPdgCode());
1198 Int_t ks = iparticle->GetStatusCode();
1199 Int_t km = iparticle->GetFirstMother();
06956108 1200 if (
1201 (((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
1202 ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
1203 )
1204 {
1205 nc++;
8d2cd130 1206 if (ks == 1) trackIt = 1;
1207 Int_t ipa = iparticle->GetFirstMother()-1;
1208
1209 iparent = (ipa > -1) ? pParent[ipa] : -1;
1210
1211//
1212// store track information
1213 p[0] = iparticle->Px();
1214 p[1] = iparticle->Py();
1215 p[2] = iparticle->Pz();
a920faf9 1216 p[3] = iparticle->Energy();
1406f599 1217
a920faf9 1218
2590ccf9 1219 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1220 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1221 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1222
21391258 1223 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
a920faf9 1224
1225 PushTrack(fTrackIt*trackIt, iparent, kf,
1226 p[0], p[1], p[2], p[3],
1227 origin[0], origin[1], origin[2], tof,
1228 polar[0], polar[1], polar[2],
1229 kPPrimary, nt, 1., ks);
dbd64db6 1230 fNprimaries++;
8d2cd130 1231 KeepTrack(nt);
1232 pParent[i] = nt;
f913ec4f 1233 SetHighWaterMark(nt);
1234
8d2cd130 1235 } // select particle
1236 } // particle loop
1237
234f6d32 1238 delete[] pParent;
e0e89f40 1239
f913ec4f 1240 return 1;
8d2cd130 1241}
1242
1243
1244void AliGenPythia::FinishRun()
1245{
1246// Print x-section summary
1247 fPythia->Pystat(1);
2779fc64 1248
1249 if (fNev > 0.) {
1250 fQ /= fNev;
1251 fX1 /= fNev;
1252 fX2 /= fNev;
1253 }
1254
8d2cd130 1255 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1256 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1257}
1258
7184e472 1259void AliGenPythia::AdjustWeights() const
8d2cd130 1260{
1261// Adjust the weights after generation of all events
1262//
e2bddf81 1263 if (gAlice) {
1264 TParticle *part;
1265 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1266 for (Int_t i=0; i<ntrack; i++) {
1267 part= gAlice->GetMCApp()->Particle(i);
1268 part->SetWeight(part->GetWeight()*fKineBias);
1269 }
8d2cd130 1270 }
1271}
1272
20e47f08 1273void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
8d2cd130 1274{
1275// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1276
1a626d4e 1277 fAProjectile = a1;
1278 fATarget = a2;
66f02a7f 1279 fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
1d568bc2 1280 fSetNuclei = kTRUE;
8d2cd130 1281}
1282
1283
1284void AliGenPythia::MakeHeader()
1285{
7184e472 1286//
1287// Make header for the simulated event
1288//
183a5ca9 1289 if (gAlice) {
1290 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1291 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1292 }
1293
8d2cd130 1294// Builds the event header, to be called after each event
e5c87a3d 1295 if (fHeader) delete fHeader;
1296 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1297//
1298// Event type
800be8b7 1299 if(fProcDiff>0){
1300 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1301 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1302 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1303 }
1304 else
e5c87a3d 1305 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1306//
1307// Number of trials
e5c87a3d 1308 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1309//
1310// Event Vertex
d25cfd65 1311 fHeader->SetPrimaryVertex(fVertex);
21391258 1312 fHeader->SetInteractionTime(fTime+fEventTime);
dbd64db6 1313//
1314// Number of primaries
1315 fHeader->SetNProduced(fNprimaries);
8d2cd130 1316//
1317// Jets that have triggered
f913ec4f 1318
9dfe63b3 1319 //Need to store jets for b-jet studies too!
2f405d65 1320 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
8d2cd130 1321 {
1322 Int_t ntrig, njet;
1323 Float_t jets[4][10];
1324 GetJets(njet, ntrig, jets);
9ff6c04c 1325
8d2cd130 1326
1327 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1328 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1329 jets[3][i]);
1330 }
1331 }
5fa4b20b 1332//
1333// Copy relevant information from external header, if present.
1334//
1335 Float_t uqJet[4];
1336
1337 if (fRL) {
1338 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1339 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1340 {
1341 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1342
1343
1344 exHeader->TriggerJet(i, uqJet);
1345 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1346 }
1347 }
1348//
1349// Store quenching parameters
1350//
1351 if (fQuench){
b6262a45 1352 Double_t z[4] = {0.};
1353 Double_t xp = 0.;
1354 Double_t yp = 0.;
1355
7c21f297 1356 if (fQuench == 1) {
1357 // Pythia::Quench()
1358 fPythia->GetQuenchingParameters(xp, yp, z);
1bab4b79 1359 } else if (fQuench == 2){
7c21f297 1360 // Pyquen
1361 Double_t r1 = PARIMP.rb1;
1362 Double_t r2 = PARIMP.rb2;
1363 Double_t b = PARIMP.b1;
1364 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1365 Double_t phi = PARIMP.psib1;
1366 xp = r * TMath::Cos(phi);
1367 yp = r * TMath::Sin(phi);
1368
1bab4b79 1369 } else if (fQuench == 4) {
1370 // QPythia
5831e75f 1371 Double_t xy[2];
e9719084 1372 Double_t i0i1[2];
1373 AliFastGlauber::Instance()->GetSavedXY(xy);
1374 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
5831e75f 1375 xp = xy[0];
1376 yp = xy[1];
e6fe9b82 1377 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
7c21f297 1378 }
1bab4b79 1379
7c21f297 1380 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1381 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1bab4b79 1382 }
beac474c 1383//
1384// Store pt^hard
1385 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 1386//
cf57b268 1387// Pass header
5fa4b20b 1388//
cf57b268 1389 AddHeader(fHeader);
4c4eac97 1390 fHeader = 0x0;
8d2cd130 1391}
cf57b268 1392
c2fc9d31 1393Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
8d2cd130 1394{
1395// Check the kinematic trigger condition
1396//
1397 Double_t eta[2];
1398 eta[0] = jet1->Eta();
1399 eta[1] = jet2->Eta();
1400 Double_t phi[2];
1401 phi[0] = jet1->Phi();
1402 phi[1] = jet2->Phi();
1403 Int_t pdg[2];
1404 pdg[0] = jet1->GetPdgCode();
1405 pdg[1] = jet2->GetPdgCode();
1406 Bool_t triggered = kFALSE;
1407
2f405d65 1408 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8d2cd130 1409 Int_t njets = 0;
1410 Int_t ntrig = 0;
1411 Float_t jets[4][10];
1412//
1413// Use Pythia clustering on parton level to determine jet axis
1414//
1415 GetJets(njets, ntrig, jets);
1416
76d6ba9a 1417 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1418//
1419 } else {
1420 Int_t ij = 0;
1421 Int_t ig = 1;
1422 if (pdg[0] == kGamma) {
1423 ij = 1;
1424 ig = 0;
1425 }
1426 //Check eta range first...
1427 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1428 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1429 {
1430 //Eta is okay, now check phi range
1431 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1432 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1433 {
1434 triggered = kTRUE;
1435 }
1436 }
1437 }
1438 return triggered;
1439}
aea21c57 1440
1441
aea21c57 1442
7184e472 1443Bool_t AliGenPythia::CheckKinematicsOnChild(){
1444//
1445//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1446//
aea21c57 1447 Bool_t checking = kFALSE;
1448 Int_t j, kcode, ks, km;
1449 Int_t nPartAcc = 0; //number of particles in the acceptance range
1450 Int_t numberOfAcceptedParticles = 1;
1451 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1452 Int_t npart = fParticles.GetEntriesFast();
aea21c57 1453
0f6ee828 1454 for (j = 0; j<npart; j++) {
8507138f 1455 TParticle * jparticle = (TParticle *) fParticles.At(j);
aea21c57 1456 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1457 ks = jparticle->GetStatusCode();
1458 km = jparticle->GetFirstMother();
1459
1460 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1461 nPartAcc++;
1462 }
0f6ee828 1463 if( numberOfAcceptedParticles <= nPartAcc){
1464 checking = kTRUE;
1465 break;
1466 }
aea21c57 1467 }
0f6ee828 1468
aea21c57 1469 return checking;
aea21c57 1470}
1471
5fa4b20b 1472void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1473{
1058d9df 1474 //
1475 // Load event into Pythia Common Block
1476 //
1477
1478 Int_t npart = stack -> GetNprimary();
1479 Int_t n0 = 0;
1480
1481 if (!flag) {
1482 (fPythia->GetPyjets())->N = npart;
1483 } else {
1484 n0 = (fPythia->GetPyjets())->N;
1485 (fPythia->GetPyjets())->N = n0 + npart;
1486 }
1487
1488
1489 for (Int_t part = 0; part < npart; part++) {
1490 TParticle *mPart = stack->Particle(part);
32d6ef7d 1491
1058d9df 1492 Int_t kf = mPart->GetPdgCode();
1493 Int_t ks = mPart->GetStatusCode();
1494 Int_t idf = mPart->GetFirstDaughter();
1495 Int_t idl = mPart->GetLastDaughter();
1496
1497 if (reHadr) {
1498 if (ks == 11 || ks == 12) {
1499 ks -= 10;
1500 idf = -1;
1501 idl = -1;
1502 }
32d6ef7d 1503 }
1504
1058d9df 1505 Float_t px = mPart->Px();
1506 Float_t py = mPart->Py();
1507 Float_t pz = mPart->Pz();
1508 Float_t e = mPart->Energy();
1509 Float_t m = mPart->GetCalcMass();
32d6ef7d 1510
1058d9df 1511
1512 (fPythia->GetPyjets())->P[0][part+n0] = px;
1513 (fPythia->GetPyjets())->P[1][part+n0] = py;
1514 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1515 (fPythia->GetPyjets())->P[3][part+n0] = e;
1516 (fPythia->GetPyjets())->P[4][part+n0] = m;
1517
1518 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1519 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1520 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1521 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1522 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1523 }
1524}
1525
c2fc9d31 1526void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1058d9df 1527{
1528 //
1529 // Load event into Pythia Common Block
1530 //
1531
1532 Int_t npart = stack -> GetEntries();
1533 Int_t n0 = 0;
1534
1535 if (!flag) {
1536 (fPythia->GetPyjets())->N = npart;
1537 } else {
1538 n0 = (fPythia->GetPyjets())->N;
1539 (fPythia->GetPyjets())->N = n0 + npart;
1540 }
1541
1542
1543 for (Int_t part = 0; part < npart; part++) {
1544 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
92847124 1545 if (!mPart) continue;
1546
1058d9df 1547 Int_t kf = mPart->GetPdgCode();
1548 Int_t ks = mPart->GetStatusCode();
1549 Int_t idf = mPart->GetFirstDaughter();
1550 Int_t idl = mPart->GetLastDaughter();
1551
1552 if (reHadr) {
92847124 1553 if (ks == 11 || ks == 12) {
1554 ks -= 10;
1555 idf = -1;
1556 idl = -1;
1557 }
8d2cd130 1558 }
1058d9df 1559
1560 Float_t px = mPart->Px();
1561 Float_t py = mPart->Py();
1562 Float_t pz = mPart->Pz();
1563 Float_t e = mPart->Energy();
1564 Float_t m = mPart->GetCalcMass();
1565
1566
1567 (fPythia->GetPyjets())->P[0][part+n0] = px;
1568 (fPythia->GetPyjets())->P[1][part+n0] = py;
1569 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1570 (fPythia->GetPyjets())->P[3][part+n0] = e;
1571 (fPythia->GetPyjets())->P[4][part+n0] = m;
1572
1573 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1574 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1575 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1576 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1577 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1578 }
8d2cd130 1579}
1580
5fa4b20b 1581
014a9521 1582void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1583{
1584//
1585// Calls the Pythia jet finding algorithm to find jets in the current event
1586//
1587//
8d2cd130 1588//
1589// Save jets
1590 Int_t n = fPythia->GetN();
1591
1592//
1593// Run Jet Finder
1594 fPythia->Pycell(njets);
1595 Int_t i;
1596 for (i = 0; i < njets; i++) {
1597 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1598 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1599 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1600 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1601
1602 jets[0][i] = px;
1603 jets[1][i] = py;
1604 jets[2][i] = pz;
1605 jets[3][i] = e;
1606 }
1607}
1608
1609
1610
1611void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1612{
1613//
1614// Calls the Pythia clustering algorithm to find jets in the current event
1615//
1616 Int_t n = fPythia->GetN();
1617 nJets = 0;
1618 nJetsTrig = 0;
1619 if (fJetReconstruction == kCluster) {
1620//
1621// Configure cluster algorithm
1622//
1623 fPythia->SetPARU(43, 2.);
1624 fPythia->SetMSTU(41, 1);
1625//
1626// Call cluster algorithm
1627//
1628 fPythia->Pyclus(nJets);
1629//
1630// Loading jets from common block
1631//
1632 } else {
592f8307 1633
8d2cd130 1634//
1635// Run Jet Finder
1636 fPythia->Pycell(nJets);
1637 }
1638
1639 Int_t i;
1640 for (i = 0; i < nJets; i++) {
1641 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1642 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1643 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1644 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1645 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1646 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1647 Float_t theta = TMath::ATan2(pt,pz);
1648 Float_t et = e * TMath::Sin(theta);
1649 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1650 if (
1651 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1652 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1653 et > fEtMinJet && et < fEtMaxJet
1654 )
1655 {
1656 jets[0][nJetsTrig] = px;
1657 jets[1][nJetsTrig] = py;
1658 jets[2][nJetsTrig] = pz;
1659 jets[3][nJetsTrig] = e;
1660 nJetsTrig++;
5fa4b20b 1661// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1662 } else {
1663// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1664 }
1665 }
1666}
1667
f913ec4f 1668void AliGenPythia::GetSubEventTime()
1669{
1670 // Calculates time of the next subevent
9d7108a7 1671 fEventTime = 0.;
1672 if (fEventsTime) {
1673 TArrayF &array = *fEventsTime;
1674 fEventTime = array[fCurSubEvent++];
1675 }
1676 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1677 return;
f913ec4f 1678}
8d2cd130 1679
c2fc9d31 1680Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
ec2c406e 1681{
1682 // Is particle in EMCAL acceptance?
ec2c406e 1683 // phi in degrees, etamin=-etamax
9fd8e520 1684 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1685 eta < fEMCALEta )
ec2c406e 1686 return kTRUE;
1687 else
1688 return kFALSE;
1689}
1690
c2fc9d31 1691Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta) const
ec2c406e 1692{
1693 // Is particle in PHOS acceptance?
1694 // Acceptance slightly larger considered.
1695 // phi in degrees, etamin=-etamax
9fd8e520 1696 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1697 eta < fPHOSEta )
ec2c406e 1698 return kTRUE;
1699 else
1700 return kFALSE;
1701}
1702
90a236ce 1703void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1704{
1705 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1706 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1707 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
800be8b7 1708 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
90a236ce 1709
1710 //calculate deltaphi
8507138f 1711 TParticle* ph = (TParticle *) fParticles.At(iphcand);
90a236ce 1712 Double_t phphi = ph->Phi();
1713 Double_t deltaphi = phiPHOS - phphi;
1714
1715
1716
1717 //loop for all particles and produce the phi rotation
8507138f 1718 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
90a236ce 1719 Double_t oldphi, newphi;
777e69b0 1720 Double_t newVx, newVy, r, vZ, time;
1721 Double_t newPx, newPy, pt, pz, e;
90a236ce 1722 for(Int_t i=0; i< np; i++) {
8507138f 1723 TParticle* iparticle = (TParticle *) fParticles.At(i);
90a236ce 1724 oldphi = iparticle->Phi();
1725 newphi = oldphi + deltaphi;
1726 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1727 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1728
777e69b0 1729 r = iparticle->R();
1730 newVx = r * TMath::Cos(newphi);
1731 newVy = r * TMath::Sin(newphi);
1732 vZ = iparticle->Vz(); // don't transform
90a236ce 1733 time = iparticle->T(); // don't transform
1734
1735 pt = iparticle->Pt();
777e69b0 1736 newPx = pt * TMath::Cos(newphi);
1737 newPy = pt * TMath::Sin(newphi);
1738 pz = iparticle->Pz(); // don't transform
1739 e = iparticle->Energy(); // don't transform
90a236ce 1740
1741 // apply rotation
777e69b0 1742 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1743 iparticle->SetMomentum(newPx, newPy, pz, e);
90a236ce 1744
1745 } //end particle loop
1746
1747 // now let's check that we put correctly the candidate photon in PHOS
1748 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1749 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1750 if(IsInPHOS(phi,eta))
1751 okdd = kTRUE;
1752}
ec2c406e 1753
1754
589380c6 1755
c5dfa3e4 1756
800be8b7 1757Bool_t AliGenPythia::CheckDiffraction()
1758{
1759 // use this method only with Perugia-0 tune!
1760
1761 // printf("AAA\n");
1762
1763 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1764
1765 Int_t iPart1=-1;
1766 Int_t iPart2=-1;
1767
1768 Double_t y1 = 1e10;
1769 Double_t y2 = -1e10;
1770
1771 const Int_t kNstable=20;
1772 const Int_t pdgStable[20] = {
1773 22, // Photon
1774 11, // Electron
1775 12, // Electron Neutrino
1776 13, // Muon
1777 14, // Muon Neutrino
1778 15, // Tau
1779 16, // Tau Neutrino
1780 211, // Pion
1781 321, // Kaon
1782 311, // K0
1783 130, // K0s
1784 310, // K0l
1785 2212, // Proton
1786 2112, // Neutron
1787 3122, // Lambda_0
1788 3112, // Sigma Minus
1789 3222, // Sigma Plus
1790 3312, // Xsi Minus
1791 3322, // Xsi0
1792 3334 // Omega
1793 };
1794
1795 for (Int_t i = 0; i < np; i++) {
1796 TParticle * part = (TParticle *) fParticles.At(i);
1797
1798 Int_t statusCode = part->GetStatusCode();
1799
1800 // Initial state particle
1801 if (statusCode != 1)
1802 continue;
1803
1804 Int_t pdg = TMath::Abs(part->GetPdgCode());
1805 Bool_t isStable = kFALSE;
1806 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1807 if (pdg == pdgStable[i1]) {
1808 isStable = kTRUE;
1809 break;
1810 }
1811 }
1812 if(!isStable)
1813 continue;
1814
1815 Double_t y = part->Y();
1816
1817 if (y < y1)
1818 {
1819 y1 = y;
1820 iPart1 = i;
1821 }
1822 if (y > y2)
1823 {
1824 y2 = y;
1825 iPart2 = i;
1826 }
1827 }
1828
1829 if(iPart1<0 || iPart2<0) return kFALSE;
1830
1831 y1=TMath::Abs(y1);
1832 y2=TMath::Abs(y2);
1833
1834 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1835 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1836
1837 Int_t pdg1 = part1->GetPdgCode();
1838 Int_t pdg2 = part2->GetPdgCode();
1839
1840
1841 Int_t iPart = -1;
1842 if (pdg1 == 2212 && pdg2 == 2212)
1843 {
1844 if(y1 > y2)
1845 iPart = iPart1;
1846 else if(y1 < y2)
1847 iPart = iPart2;
1848 else {
1849 iPart = iPart1;
1850 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1851 }
1852 }
1853 else if (pdg1 == 2212)
1854 iPart = iPart1;
1855 else if (pdg2 == 2212)
1856 iPart = iPart2;
1857
1858
1859
1860
1861
1862 Double_t M=-1.;
1863 if(iPart>0) {
1864 TParticle * part = (TParticle *) fParticles.At(iPart);
1865 Double_t E= part->Energy();
1866 Double_t P= part->P();
1867 M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
1868 }
1869
c5dfa3e4 1870 Double_t Mmin, Mmax, wSD, wDD, wND;
1871 if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1872
1873 if(M>-1 && M<Mmin) return kFALSE;
1874 if(M>Mmax) M=-1;
1875
1876 Int_t procType=fPythia->GetMSTI(1);
1877 Int_t proc0=2;
1878 if(procType== 94) proc0=1;
1879 if(procType== 92 || procType== 93) proc0=0;
1880
1881 Int_t proc=2;
1882 if(M>0) proc=0;
1883 else if(proc0==1) proc=1;
1884
1885 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1886 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
800be8b7 1887
1888
1889 // if(proc==1 || proc==2) return kFALSE;
1890
c5dfa3e4 1891 if(proc!=0) {
1892 if(proc0!=0) fProcDiff = procType;
1893 else fProcDiff = 95;
1894 return kTRUE;
1895 }
1896
1897 if(wSD<0) AliError("wSD<0 ! \n");
1898
1899 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1900
1901 // printf("iPart = %d\n", iPart);
1902
1903 if(iPart==iPart1) fProcDiff=93;
1904 else if(iPart==iPart2) fProcDiff=92;
1905 else {
1906 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1907
800be8b7 1908 }
1909
c5dfa3e4 1910 return kTRUE;
1911}
1912
1913
1914
1915Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1916 Double_t &wSD, Double_t &wDD, Double_t &wND)
1917{
1918
1919 // 900 GeV
1920 if(TMath::Abs(fEnergyCMS-900)<1 ){
1921
1922const Int_t nbin=400;
1923Double_t bin[]={
19241.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
19254.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
19267.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
192710.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
192813.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
192915.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
193018.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
193121.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
193224.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
193327.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
193430.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
193533.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
193636.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
193739.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
193842.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
193945.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
194048.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
194151.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
194254.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
194357.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
194460.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
194563.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
194666.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
194769.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
194872.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
194975.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
195078.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
195181.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
195284.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
195387.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
195490.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
195593.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
195696.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
195799.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1958102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1959105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1960108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1961111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1962114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1963117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1964120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1965123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1966126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
1967129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
1968132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
1969135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
1970138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
1971141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
1972144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
1973147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
1974150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
1975153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
1976156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
1977159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
1978162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
1979165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
1980168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
1981171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
1982174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
1983177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
1984180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
1985183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
1986186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
1987189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
1988192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
1989195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
1990198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1991Double_t w[]={
19921.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
19930.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
19940.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
19950.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
19960.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
19970.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
19980.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
19990.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
20000.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
20010.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
20020.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
20030.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
20040.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
20050.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
20060.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
20070.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
20080.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
20090.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
20100.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
20110.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
20120.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
20130.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
20140.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
20150.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
20160.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
20170.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
20180.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
20190.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
20200.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
20210.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
20220.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
20230.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
20240.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
20250.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
20260.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
20270.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
20280.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
20290.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
20300.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
20310.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
20320.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
20330.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
20340.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
20350.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
20360.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
20370.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
20380.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
20390.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
20400.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
20410.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
20420.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
20430.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
20440.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
20450.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
20460.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
20470.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
20480.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
20490.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
20500.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
20510.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
20520.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
20530.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
20540.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
20550.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
20560.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
20570.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
20580.386112, 0.364314, 0.434375, 0.334629};
2059wDD = 0.379611;
2060wND = 0.496961;
2061wSD = -1;
2062
2063 Mmin = bin[0];
2064 Mmax = bin[nbin];
2065 if(M<Mmin || M>Mmax) return kTRUE;
2066
7873275c 2067 Int_t ibin=nbin-1;
93a60586 2068 for(Int_t i=1; i<=nbin; i++)
2ff57420 2069 if(M<=bin[i]) {
2070 ibin=i-1;
2071 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
800be8b7 2072 break;
2073 }
c5dfa3e4 2074 wSD=w[ibin];
2075 return kTRUE;
2076 }
2077 else if(TMath::Abs(fEnergyCMS-2760)<1 ){
800be8b7 2078
c5dfa3e4 2079const Int_t nbin=400;
2080Double_t bin[]={
20811.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
20824.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
20837.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
208410.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
208513.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
208615.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
208718.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
208821.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
208924.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
209027.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
209130.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
209233.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
209336.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
209439.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
209542.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
209645.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
209748.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
209851.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
209954.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
210057.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
210160.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
210263.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
210366.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
210469.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
210572.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
210675.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
210778.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
210881.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
210984.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
211087.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
211190.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
211293.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
211396.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
211499.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2115102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2116105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2117108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2118111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2119114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2120117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2121120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2122123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2123126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2124129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2125132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2126135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2127138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2128141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2129144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2130147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2131150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2132153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2133156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2134159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2135162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2136165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2137168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2138171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2139174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2140177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2141180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2142183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2143186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2144189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2145192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2146195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2147198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2148Double_t w[]={
21491.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
21500.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
21510.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
21520.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
21530.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
21540.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
21550.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
21560.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
21570.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
21580.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
21590.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
21600.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
21610.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
21620.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
21630.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
21640.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
21650.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
21660.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
21670.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
21680.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
21690.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
21700.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
21710.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
21720.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
21730.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
21740.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
21750.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
21760.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
21770.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
21780.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
21790.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
21800.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
21810.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
21820.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
21830.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
21840.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
21850.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
21860.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
21870.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
21880.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
21890.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
21900.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
21910.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
21920.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
21930.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
21940.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
21950.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
21960.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
21970.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
21980.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
21990.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
22000.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
22010.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
22020.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
22030.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
22040.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
22050.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
22060.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
22070.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
22080.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
22090.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
22100.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
22110.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
22120.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
22130.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
22140.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
22150.201712, 0.242225, 0.255565, 0.258738};
2216wDD = 0.512813;
2217wND = 0.518820;
2218wSD = -1;
2219
2220 Mmin = bin[0];
2221 Mmax = bin[nbin];
2222 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2223
c5dfa3e4 2224 Int_t ibin=nbin-1;
2225 for(Int_t i=1; i<=nbin; i++)
2226 if(M<=bin[i]) {
2227 ibin=i-1;
2228 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2229 break;
2230 }
2231 wSD=w[ibin];
2232 return kTRUE;
2233 }
800be8b7 2234
800be8b7 2235
c5dfa3e4 2236 else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2237const Int_t nbin=400;
2238Double_t bin[]={
22391.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
22404.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
22417.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
224210.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
224313.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
224415.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
224518.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
224621.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
224724.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
224827.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
224930.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
225033.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
225136.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
225239.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
225342.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
225445.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
225548.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
225651.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
225754.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
225857.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
225960.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
226063.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
226166.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
226269.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
226372.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
226475.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
226578.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
226681.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
226784.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
226887.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
226990.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
227093.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
227196.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
227299.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2273102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2274105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2275108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2276111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2277114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2278117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2279120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2280123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2281126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2282129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2283132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2284135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2285138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2286141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2287144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2288147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2289150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2290153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2291156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2292159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2293162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2294165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2295168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2296171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2297174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2298177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2299180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2300183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2301186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2302189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2303192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2304195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2305198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2306Double_t w[]={
23071.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
23080.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
23090.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
23100.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
23110.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
23120.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
23130.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
23140.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
23150.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
23160.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
23170.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
23180.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
23190.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
23200.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
23210.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
23220.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
23230.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
23240.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
23250.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
23260.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
23270.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
23280.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
23290.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
23300.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
23310.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
23320.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
23330.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
23340.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
23350.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
23360.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
23370.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
23380.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
23390.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
23400.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
23410.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
23420.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
23430.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
23440.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
23450.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
23460.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
23470.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
23480.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
23490.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
23500.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
23510.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
23520.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
23530.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
23540.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
23550.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
23560.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
23570.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
23580.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
23590.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
23600.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
23610.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
23620.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
23630.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
23640.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
23650.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
23660.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
23670.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
23680.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
23690.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
23700.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
23710.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
23720.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
23730.175006, 0.223482, 0.202706, 0.218108};
2374wDD = 0.207705;
2375wND = 0.289628;
2376wSD = -1;
2377
2378 Mmin = bin[0];
2379 Mmax = bin[nbin];
2380
2381 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2382
c5dfa3e4 2383 Int_t ibin=nbin-1;
2384 for(Int_t i=1; i<=nbin; i++)
2385 if(M<=bin[i]) {
2386 ibin=i-1;
2387 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2388 break;
2389 }
2390 wSD=w[ibin];
800be8b7 2391 return kTRUE;
c5dfa3e4 2392 }
2393
2394 return kFALSE;
800be8b7 2395}
06956108 2396
2397Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2398{
2399// Check if this is a heavy flavor decay product
2400 TParticle * part = (TParticle *) fParticles.At(ipart);
2401 Int_t mpdg = TMath::Abs(part->GetPdgCode());
2402 Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2403 //
2404 // Light hadron
2405 if (mfl >= 4 && mfl < 6) return kTRUE;
2406
2407 Int_t imo = part->GetFirstMother()-1;
2408 TParticle* pm = part;
2409 //
2410 // Heavy flavor hadron produced by generator
2411 while (imo > -1) {
2412 pm = (TParticle*)fParticles.At(imo);
2413 mpdg = TMath::Abs(pm->GetPdgCode());
2414 mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2415 if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2416 imo = pm->GetFirstMother()-1;
2417 }
2418 return kFALSE;
2419}