Update timestamp for new data points simulation
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia6 / 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"
8d2cd130 49
014a9521 50ClassImp(AliGenPythia)
8d2cd130 51
e8a8adcd 52
53AliGenPythia::AliGenPythia():
54 AliGenMC(),
55 fProcess(kPyCharm),
efe3b1cd 56 fItune(-1),
e8a8adcd 57 fStrucFunc(kCTEQ5L),
e8a8adcd 58 fKineBias(0.),
59 fTrials(0),
60 fTrialsRun(0),
61 fQ(0.),
62 fX1(0.),
63 fX2(0.),
64 fEventTime(0.),
65 fInteractionRate(0.),
66 fTimeWindow(0.),
67 fCurSubEvent(0),
68 fEventsTime(0),
69 fNev(0),
70 fFlavorSelect(0),
71 fXsection(0.),
72 fPythia(0),
1cdfb96f 73 fWeightPower(0.),
e8a8adcd 74 fPtHardMin(0.),
75 fPtHardMax(1.e4),
76 fYHardMin(-1.e10),
77 fYHardMax(1.e10),
78 fGinit(1),
79 fGfinal(1),
036662e5 80 fCRoff(0),
e8a8adcd 81 fHadronisation(1),
03358a32 82 fPatchOmegaDalitz(0),
d7d0c637 83 fDecayerExodus(0),
e8a8adcd 84 fNpartons(0),
85 fReadFromFile(0),
64da86aa 86 fReadLHEF(0),
e8a8adcd 87 fQuench(0),
cd07c39b 88 fQhat(0.),
89 fLength(0.),
66b8652c 90 fpyquenT(1.),
91 fpyquenTau(0.1),
92 fpyquenNf(0),
93 fpyquenEloss(0),
94 fpyquenAngle(0),
e6fe9b82 95 fImpact(0.),
e8a8adcd 96 fPtKick(1.),
97 fFullEvent(kTRUE),
98 fDecayer(new AliDecayerPythia()),
99 fDebugEventFirst(-1),
100 fDebugEventLast(-1),
101 fEtMinJet(0.),
102 fEtMaxJet(1.e4),
103 fEtaMinJet(-20.),
104 fEtaMaxJet(20.),
105 fPhiMinJet(0.),
106 fPhiMaxJet(2.* TMath::Pi()),
107 fJetReconstruction(kCell),
108 fEtaMinGamma(-20.),
109 fEtaMaxGamma(20.),
110 fPhiMinGamma(0.),
111 fPhiMaxGamma(2. * TMath::Pi()),
d8850d4f 112 fUseYCutHQ(kFALSE),
113 fYMinHQ(-20.),
114 fYMaxHQ(20.),
e8a8adcd 115 fPycellEtaMax(2.),
116 fPycellNEta(274),
117 fPycellNPhi(432),
118 fPycellThreshold(0.),
119 fPycellEtSeed(4.),
120 fPycellMinEtJet(10.),
121 fPycellMaxRadius(1.),
122 fStackFillOpt(kFlavorSelection),
123 fFeedDownOpt(kTRUE),
124 fFragmentation(kTRUE),
125 fSetNuclei(kFALSE),
fb355e51 126 fUseNuclearPDF(kFALSE),
127 fUseLorentzBoost(kTRUE),
e8a8adcd 128 fNewMIS(kFALSE),
129 fHFoff(kFALSE),
20e47f08 130 fNucPdf(0),
e8a8adcd 131 fTriggerParticle(0),
132 fTriggerEta(0.9),
4b5ec27d 133 fTriggerY(999.),
1f0caa6a 134 fTriggerEtaMin(0.9),
2591bd0e 135 fTriggerMinPt(-1),
136 fTriggerMaxPt(1000),
700b9416 137 fTriggerMultiplicity(0),
138 fTriggerMultiplicityEta(0),
440c2873 139 fTriggerMultiplicityEtaMin(0),
140 fTriggerMultiplicityEtaMax(0),
38112f3f 141 fTriggerMultiplicityPtMin(0),
e8a8adcd 142 fCountMode(kCountAll),
143 fHeader(0),
144 fRL(0),
777e69b0 145 fkFileName(0),
64da86aa 146 fkNameLHEF(0),
9fd8e520 147 fFragPhotonInCalo(kFALSE),
43e3b80a 148 fHadronInCalo(kFALSE) ,
ec2c406e 149 fPi0InCalo(kFALSE) ,
40fe59d4 150 fEtaInCalo(kFALSE) ,
151 fPhotonInCalo(kFALSE), // not in use
152 fDecayPhotonInCalo(kFALSE),
153 fForceNeutralMeson2PhotonDecay(kFALSE),
154 fEleInCalo(kFALSE),
155 fEleInEMCAL(kFALSE), // not in use
156 fCheckBarrel(kFALSE),
9fd8e520 157 fCheckEMCAL(kFALSE),
158 fCheckPHOS(kFALSE),
90a236ce 159 fCheckPHOSeta(kFALSE),
40fe59d4 160 fPHOSRotateCandidate(-1),
43e3b80a 161 fTriggerParticleMinPt(0),
40fe59d4 162 fPhotonMinPt(0), // not in use
163 fElectronMinPt(0), // not in use
9fd8e520 164 fPHOSMinPhi(219.),
165 fPHOSMaxPhi(321.),
166 fPHOSEta(0.13),
167 fEMCALMinPhi(79.),
168 fEMCALMaxPhi(191.),
800be8b7 169 fEMCALEta(0.71),
170 fkTuneForDiff(0),
171 fProcDiff(0)
8d2cd130 172{
173// Default Constructor
e7c989e4 174 fEnergyCMS = 5500.;
7cdba479 175 if (!AliPythiaRndm::GetPythiaRandom())
e8a8adcd 176 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 177}
178
179AliGenPythia::AliGenPythia(Int_t npart)
e8a8adcd 180 :AliGenMC(npart),
181 fProcess(kPyCharm),
efe3b1cd 182 fItune(-1),
e8a8adcd 183 fStrucFunc(kCTEQ5L),
e8a8adcd 184 fKineBias(0.),
185 fTrials(0),
186 fTrialsRun(0),
187 fQ(0.),
188 fX1(0.),
189 fX2(0.),
190 fEventTime(0.),
191 fInteractionRate(0.),
192 fTimeWindow(0.),
193 fCurSubEvent(0),
194 fEventsTime(0),
195 fNev(0),
196 fFlavorSelect(0),
197 fXsection(0.),
198 fPythia(0),
1cdfb96f 199 fWeightPower(0.),
e8a8adcd 200 fPtHardMin(0.),
201 fPtHardMax(1.e4),
202 fYHardMin(-1.e10),
203 fYHardMax(1.e10),
204 fGinit(kTRUE),
205 fGfinal(kTRUE),
036662e5 206 fCRoff(kFALSE),
e8a8adcd 207 fHadronisation(kTRUE),
d7d0c637 208 fPatchOmegaDalitz(0),
209 fDecayerExodus(0),
e8a8adcd 210 fNpartons(0),
211 fReadFromFile(kFALSE),
64da86aa 212 fReadLHEF(0),
e8a8adcd 213 fQuench(kFALSE),
cd07c39b 214 fQhat(0.),
215 fLength(0.),
66b8652c 216 fpyquenT(1.),
217 fpyquenTau(0.1),
218 fpyquenNf(0),
219 fpyquenEloss(0),
220 fpyquenAngle(0),
e6fe9b82 221 fImpact(0.),
e8a8adcd 222 fPtKick(1.),
223 fFullEvent(kTRUE),
224 fDecayer(new AliDecayerPythia()),
225 fDebugEventFirst(-1),
226 fDebugEventLast(-1),
227 fEtMinJet(0.),
228 fEtMaxJet(1.e4),
229 fEtaMinJet(-20.),
230 fEtaMaxJet(20.),
231 fPhiMinJet(0.),
232 fPhiMaxJet(2.* TMath::Pi()),
233 fJetReconstruction(kCell),
234 fEtaMinGamma(-20.),
235 fEtaMaxGamma(20.),
236 fPhiMinGamma(0.),
237 fPhiMaxGamma(2. * TMath::Pi()),
d8850d4f 238 fUseYCutHQ(kFALSE),
239 fYMinHQ(-20.),
240 fYMaxHQ(20.),
e8a8adcd 241 fPycellEtaMax(2.),
242 fPycellNEta(274),
243 fPycellNPhi(432),
244 fPycellThreshold(0.),
245 fPycellEtSeed(4.),
246 fPycellMinEtJet(10.),
247 fPycellMaxRadius(1.),
248 fStackFillOpt(kFlavorSelection),
249 fFeedDownOpt(kTRUE),
250 fFragmentation(kTRUE),
251 fSetNuclei(kFALSE),
fb355e51 252 fUseNuclearPDF(kFALSE),
253 fUseLorentzBoost(kTRUE),
e8a8adcd 254 fNewMIS(kFALSE),
255 fHFoff(kFALSE),
20e47f08 256 fNucPdf(0),
e8a8adcd 257 fTriggerParticle(0),
2591bd0e 258 fTriggerEta(0.9),
4b5ec27d 259 fTriggerY(999.),
1f0caa6a 260 fTriggerEtaMin(0.9),
2591bd0e 261 fTriggerMinPt(-1),
262 fTriggerMaxPt(1000),
700b9416 263 fTriggerMultiplicity(0),
264 fTriggerMultiplicityEta(0),
440c2873 265 fTriggerMultiplicityEtaMin(0),
266 fTriggerMultiplicityEtaMax(0),
38112f3f 267 fTriggerMultiplicityPtMin(0),
e8a8adcd 268 fCountMode(kCountAll),
269 fHeader(0),
270 fRL(0),
777e69b0 271 fkFileName(0),
64da86aa 272 fkNameLHEF(0),
9fd8e520 273 fFragPhotonInCalo(kFALSE),
43e3b80a 274 fHadronInCalo(kFALSE) ,
ec2c406e 275 fPi0InCalo(kFALSE) ,
40fe59d4 276 fEtaInCalo(kFALSE) ,
277 fPhotonInCalo(kFALSE), // not in use
278 fDecayPhotonInCalo(kFALSE),
279 fForceNeutralMeson2PhotonDecay(kFALSE),
280 fEleInCalo(kFALSE),
281 fEleInEMCAL(kFALSE), // not in use
282 fCheckBarrel(kFALSE),
9fd8e520 283 fCheckEMCAL(kFALSE),
284 fCheckPHOS(kFALSE),
90a236ce 285 fCheckPHOSeta(kFALSE),
40fe59d4 286 fPHOSRotateCandidate(-1),
43e3b80a 287 fTriggerParticleMinPt(0),
40fe59d4 288 fPhotonMinPt(0), // not in use
289 fElectronMinPt(0), // not in use
9fd8e520 290 fPHOSMinPhi(219.),
291 fPHOSMaxPhi(321.),
292 fPHOSEta(0.13),
293 fEMCALMinPhi(79.),
294 fEMCALMaxPhi(191.),
800be8b7 295 fEMCALEta(0.71),
296 fkTuneForDiff(0),
297 fProcDiff(0)
8d2cd130 298{
299// default charm production at 5. 5 TeV
300// semimuonic decay
301// structure function GRVHO
302//
e7c989e4 303 fEnergyCMS = 5500.;
8d2cd130 304 fName = "Pythia";
305 fTitle= "Particle Generator using PYTHIA";
8d2cd130 306 SetForceDecay();
8d2cd130 307 // Set random number generator
7cdba479 308 if (!AliPythiaRndm::GetPythiaRandom())
309 AliPythiaRndm::SetPythiaRandom(GetRandom());
76d6ba9a 310 }
8d2cd130 311
8d2cd130 312AliGenPythia::~AliGenPythia()
313{
314// Destructor
9d7108a7 315 if(fEventsTime) delete fEventsTime;
316}
317
318void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
319{
320// Generate pileup using user specified rate
321 fInteractionRate = rate;
322 fTimeWindow = timewindow;
323 GeneratePileup();
324}
325
326void AliGenPythia::GeneratePileup()
327{
328// Generate sub events time for pileup
329 fEventsTime = 0;
330 if(fInteractionRate == 0.) {
331 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
332 return;
333 }
334
335 Int_t npart = NumberParticles();
336 if(npart < 0) {
337 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
338 return;
339 }
340
341 if(fEventsTime) delete fEventsTime;
342 fEventsTime = new TArrayF(npart);
343 TArrayF &array = *fEventsTime;
344 for(Int_t ipart = 0; ipart < npart; ipart++)
345 array[ipart] = 0.;
346
347 Float_t eventtime = 0.;
348 while(1)
349 {
350 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
351 if(eventtime > fTimeWindow) break;
352 array.Set(array.GetSize()+1);
353 array[array.GetSize()-1] = eventtime;
354 }
355
356 eventtime = 0.;
357 while(1)
358 {
359 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
360 if(TMath::Abs(eventtime) > fTimeWindow) break;
361 array.Set(array.GetSize()+1);
362 array[array.GetSize()-1] = eventtime;
363 }
364
365 SetNumberParticles(fEventsTime->GetSize());
8d2cd130 366}
367
592f8307 368void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
369 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
370{
371// Set pycell parameters
372 fPycellEtaMax = etamax;
373 fPycellNEta = neta;
374 fPycellNPhi = nphi;
375 fPycellThreshold = thresh;
376 fPycellEtSeed = etseed;
377 fPycellMinEtJet = minet;
378 fPycellMaxRadius = r;
379}
380
381
382
8d2cd130 383void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
384{
385 // Set a range of event numbers, for which a table
386 // of generated particle will be printed
387 fDebugEventFirst = eventFirst;
388 fDebugEventLast = eventLast;
389 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
390}
391
392void AliGenPythia::Init()
393{
394// Initialisation
395
396 SetMC(AliPythia::Instance());
b88f5cea 397 fPythia=(AliPythia*) fMCEvGen;
e2bddf81 398
8d2cd130 399//
400 fParentWeight=1./Float_t(fNpart);
401//
1cdfb96f 402 if (fWeightPower != 0)
403 fPythia->SetWeightPower(fWeightPower);
8d2cd130 404 fPythia->SetCKIN(3,fPtHardMin);
405 fPythia->SetCKIN(4,fPtHardMax);
406 fPythia->SetCKIN(7,fYHardMin);
407 fPythia->SetCKIN(8,fYHardMax);
408
e28ccdaf 409 if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
410
fb355e51 411 if(fUseNuclearPDF)
412 fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
8d2cd130 413 // Fragmentation?
414 if (fFragmentation) {
415 fPythia->SetMSTP(111,1);
416 } else {
417 fPythia->SetMSTP(111,0);
418 }
419
420
421// initial state radiation
422 fPythia->SetMSTP(61,fGinit);
423// final state radiation
424 fPythia->SetMSTP(71,fGfinal);
036662e5 425 //color reconnection strength
426 if(fCRoff==1)fPythia->SetMSTP(95,0);
8d2cd130 427// pt - kick
428 if (fPtKick > 0.) {
429 fPythia->SetMSTP(91,1);
688950ef 430 fPythia->SetPARP(91,fPtKick);
431 fPythia->SetPARP(93, 4. * fPtKick);
8d2cd130 432 } else {
433 fPythia->SetMSTP(91,0);
434 }
435
64da86aa 436 if (fReadLHEF) fPythia->OpenFortranFile(97, const_cast<char*>(fkNameLHEF));
437
5fa4b20b 438 if (fReadFromFile) {
777e69b0 439 fRL = AliRunLoader::Open(fkFileName, "Partons");
5fa4b20b 440 fRL->LoadKinematics();
441 fRL->LoadHeader();
442 } else {
443 fRL = 0x0;
444 }
fdea4387 445 //
efe3b1cd 446 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
fdea4387 447 // Forward Paramters to the AliPythia object
448 fDecayer->SetForceDecay(fForceDecay);
beac474c 449// Switch off Heavy Flavors on request
450 if (fHFoff) {
51387949 451 // Maximum number of quark flavours used in pdf
beac474c 452 fPythia->SetMSTP(58, 3);
51387949 453 // Maximum number of flavors that can be used in showers
beac474c 454 fPythia->SetMSTJ(45, 3);
51387949 455 // Switch off g->QQbar splitting in decay table
456 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
beac474c 457 }
fdea4387 458
51387949 459 fDecayer->Init();
460
8d2cd130 461
462// Parent and Children Selection
463 switch (fProcess)
464 {
65f2626c 465 case kPyOldUEQ2ordered:
466 case kPyOldUEQ2ordered2:
467 case kPyOldPopcorn:
468 break;
8d2cd130 469 case kPyCharm:
470 case kPyCharmUnforced:
adf4d898 471 case kPyCharmPbPbMNR:
aabc7187 472 case kPyCharmpPbMNR:
e0e89f40 473 case kPyCharmppMNR:
474 case kPyCharmppMNRwmi:
64da86aa 475 case kPyCharmPWHG:
8d2cd130 476 fParentSelect[0] = 411;
477 fParentSelect[1] = 421;
478 fParentSelect[2] = 431;
479 fParentSelect[3] = 4122;
9ccc3d0e 480 fParentSelect[4] = 4232;
481 fParentSelect[5] = 4132;
482 fParentSelect[6] = 4332;
8d2cd130 483 fFlavorSelect = 4;
484 break;
adf4d898 485 case kPyD0PbPbMNR:
486 case kPyD0pPbMNR:
487 case kPyD0ppMNR:
8d2cd130 488 fParentSelect[0] = 421;
489 fFlavorSelect = 4;
490 break;
90d7b703 491 case kPyDPlusPbPbMNR:
492 case kPyDPluspPbMNR:
493 case kPyDPlusppMNR:
494 fParentSelect[0] = 411;
495 fFlavorSelect = 4;
496 break;
e0e89f40 497 case kPyDPlusStrangePbPbMNR:
498 case kPyDPlusStrangepPbMNR:
499 case kPyDPlusStrangeppMNR:
500 fParentSelect[0] = 431;
501 fFlavorSelect = 4;
502 break;
68504d42 503 case kPyLambdacppMNR:
504 fParentSelect[0] = 4122;
505 fFlavorSelect = 4;
506 break;
8d2cd130 507 case kPyBeauty:
9dfe63b3 508 case kPyBeautyJets:
adf4d898 509 case kPyBeautyPbPbMNR:
510 case kPyBeautypPbMNR:
511 case kPyBeautyppMNR:
e0e89f40 512 case kPyBeautyppMNRwmi:
64da86aa 513 case kPyBeautyPWHG:
8d2cd130 514 fParentSelect[0]= 511;
515 fParentSelect[1]= 521;
516 fParentSelect[2]= 531;
517 fParentSelect[3]= 5122;
518 fParentSelect[4]= 5132;
519 fParentSelect[5]= 5232;
520 fParentSelect[6]= 5332;
521 fFlavorSelect = 5;
522 break;
523 case kPyBeautyUnforced:
524 fParentSelect[0] = 511;
525 fParentSelect[1] = 521;
526 fParentSelect[2] = 531;
527 fParentSelect[3] = 5122;
528 fParentSelect[4] = 5132;
529 fParentSelect[5] = 5232;
530 fParentSelect[6] = 5332;
531 fFlavorSelect = 5;
532 break;
533 case kPyJpsiChi:
534 case kPyJpsi:
535 fParentSelect[0] = 443;
536 break;
f529e69b 537 case kPyMbDefault:
0bd3d7c5 538 case kPyMbAtlasTuneMC09:
8d2cd130 539 case kPyMb:
04081a91 540 case kPyMbWithDirectPhoton:
8d2cd130 541 case kPyMbNonDiffr:
d7de4a67 542 case kPyMbMSEL1:
8d2cd130 543 case kPyJets:
64da86aa 544 case kPyJetsPWHG:
8d2cd130 545 case kPyDirectGamma:
e33acb67 546 case kPyLhwgMb:
8d2cd130 547 break;
df607629 548 case kPyWPWHG:
589380c6 549 case kPyW:
0f6ee828 550 case kPyZ:
df607629 551 case kPyZgamma:
9a8774a1 552 case kPyMBRSingleDiffraction:
553 case kPyMBRDoubleDiffraction:
554 case kPyMBRCentralDiffraction:
589380c6 555 break;
8d2cd130 556 }
557//
592f8307 558//
559// JetFinder for Trigger
560//
561// Configure detector (EMCAL like)
562//
d7de4a67 563 fPythia->SetPARU(51, fPycellEtaMax);
564 fPythia->SetMSTU(51, fPycellNEta);
565 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 566//
567// Configure Jet Finder
568//
d7de4a67 569 fPythia->SetPARU(58, fPycellThreshold);
570 fPythia->SetPARU(52, fPycellEtSeed);
571 fPythia->SetPARU(53, fPycellMinEtJet);
572 fPythia->SetPARU(54, fPycellMaxRadius);
573 fPythia->SetMSTU(54, 2);
592f8307 574//
8d2cd130 575// This counts the total number of calls to Pyevnt() per run.
576 fTrialsRun = 0;
577 fQ = 0.;
578 fX1 = 0.;
579 fX2 = 0.;
580 fNev = 0 ;
581//
1d568bc2 582//
583//
8d2cd130 584 AliGenMC::Init();
fb355e51 585
586 // Reset Lorentz boost if demanded
587 if(!fUseLorentzBoost) {
588 fDyBoost = 0;
589 Warning("Init","Demand to discard Lorentz boost.\n");
590 }
1d568bc2 591//
592//
593//
594 if (fSetNuclei) {
fb355e51 595 fDyBoost = 0;
596 Warning("Init","Deprecated function SetNuclei() used (nPDFs + no boost). Use SetProjectile + SetTarget + SetUseNuclearPDF + SetUseLorentzBoost instead.\n");
1d568bc2 597 }
cd07c39b 598 fPythia->SetPARJ(200, 0.0);
599 fPythia->SetPARJ(199, 0.0);
600 fPythia->SetPARJ(198, 0.0);
601 fPythia->SetPARJ(197, 0.0);
602
603 if (fQuench == 1) {
5fa4b20b 604 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 605 }
3a709cfa 606
66b8652c 607 if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
608
b976f7d8 609 if (fQuench == 3) {
610 // Nestor's change of the splittings
611 fPythia->SetPARJ(200, 0.8);
612 fPythia->SetMSTJ(41, 1); // QCD radiation only
613 fPythia->SetMSTJ(42, 2); // angular ordering
614 fPythia->SetMSTJ(44, 2); // option to run alpha_s
615 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
616 fPythia->SetMSTJ(50, 0); // No coherence in first branching
617 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
cd07c39b 618 } else if (fQuench == 4) {
619 // Armesto-Cunqueiro-Salgado change of the splittings.
620 AliFastGlauber* glauber = AliFastGlauber::Instance();
621 glauber->Init(2);
622 //read and store transverse almonds corresponding to differnt
623 //impact parameters.
cd07c39b 624 glauber->SetCentralityClass(0.,0.1);
625 fPythia->SetPARJ(200, 1.);
626 fPythia->SetPARJ(198, fQhat);
627 fPythia->SetPARJ(199, fLength);
cd07c39b 628 fPythia->SetMSTJ(42, 2); // angular ordering
629 fPythia->SetMSTJ(44, 2); // option to run alpha_s
cd07c39b 630 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
b976f7d8 631 }
df607629 632
68c5eace 633 if ( AliLog::GetDebugLevel("","AliGenPythia") >= 1 ) {
634 fPythia->Pystat(4);
635 fPythia->Pystat(5);
636 }
8d2cd130 637}
638
639void AliGenPythia::Generate()
640{
641// Generate one event
19ef8cf0 642 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
8d2cd130 643 fDecayer->ForceDecay();
644
13cca24a 645 Double_t polar[3] = {0,0,0};
646 Double_t origin[3] = {0,0,0};
647 Double_t p[4];
8d2cd130 648// converts from mm/c to s
649 const Float_t kconv=0.001/2.999792458e8;
650//
651 Int_t nt=0;
652 Int_t jev=0;
653 Int_t j, kf;
654 fTrials=0;
f913ec4f 655 fEventTime = 0.;
656
2590ccf9 657
8d2cd130 658
659 // Set collision vertex position
2590ccf9 660 if (fVertexSmear == kPerEvent) Vertex();
661
8d2cd130 662// event loop
663 while(1)
664 {
32d6ef7d 665//
5fa4b20b 666// Produce event
32d6ef7d 667//
32d6ef7d 668//
5fa4b20b 669// Switch hadronisation off
670//
671 fPythia->SetMSTJ(1, 0);
cd07c39b 672
673 if (fQuench ==4){
674 Double_t bimp;
675 // Quenching comes through medium-modified splitting functions.
676 AliFastGlauber::Instance()->GetRandomBHard(bimp);
e6fe9b82 677 fPythia->SetPARJ(197, bimp);
678 fImpact = bimp;
6c43eccb 679 fPythia->Qpygin0();
cd07c39b 680 }
32d6ef7d 681//
5fa4b20b 682// Either produce new event or read partons from file
683//
684 if (!fReadFromFile) {
beac474c 685 if (!fNewMIS) {
686 fPythia->Pyevnt();
687 } else {
688 fPythia->Pyevnw();
689 }
5fa4b20b 690 fNpartons = fPythia->GetN();
691 } else {
33c3c91a 692 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
693 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
5fa4b20b 694 fPythia->SetN(0);
695 LoadEvent(fRL->Stack(), 0 , 1);
696 fPythia->Pyedit(21);
697 }
698
32d6ef7d 699//
700// Run quenching routine
701//
5fa4b20b 702 if (fQuench == 1) {
703 fPythia->Quench();
704 } else if (fQuench == 2){
705 fPythia->Pyquen(208., 0, 0.);
b976f7d8 706 } else if (fQuench == 3) {
707 // Quenching is via multiplicative correction of the splittings
5fa4b20b 708 }
b976f7d8 709
32d6ef7d 710//
5fa4b20b 711// Switch hadronisation on
32d6ef7d 712//
20e47f08 713 if (fHadronisation) {
714 fPythia->SetMSTJ(1, 1);
5fa4b20b 715//
716// .. and perform hadronisation
aea21c57 717// printf("Calling hadronisation %d\n", fPythia->GetN());
03358a32 718
719 if (fPatchOmegaDalitz) {
720 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
721 fPythia->Pyexec();
722 fPythia->DalitzDecays();
723 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
724 }
d7d0c637 725
726 else if (fDecayerExodus) {
727
728 fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
729 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
730 fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 0);
731 fPythia->Pyexec();
732 fPythia->OmegaDalitz();
733 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
734 fPythia->PizeroDalitz();
735 fPythia->PhiDalitz();
736 fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 1);
737 fPythia->EtaDalitz();
738 fPythia->EtaprimeDalitz();
739 fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
740 fPythia->RhoDirect();
741 fPythia->OmegaDirect();
742 fPythia->PhiDirect();
743 fPythia->JPsiDirect();
744 }
745
746 fPythia->Pyexec();
747 }
8d2cd130 748 fTrials++;
8507138f 749 fPythia->ImportParticles(&fParticles,"All");
03358a32 750
df275cfa 751 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
e1cf8975 752 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
f64640b5 753
8d2cd130 754//
755//
756//
757 Int_t i;
758
dbd64db6 759 fNprimaries = 0;
8507138f 760 Int_t np = fParticles.GetEntriesFast();
5fa4b20b 761
7c21f297 762 if (np == 0) continue;
8d2cd130 763//
2590ccf9 764
8d2cd130 765//
766 Int_t* pParent = new Int_t[np];
767 Int_t* pSelected = new Int_t[np];
768 Int_t* trackIt = new Int_t[np];
5fa4b20b 769 for (i = 0; i < np; i++) {
8d2cd130 770 pParent[i] = -1;
771 pSelected[i] = 0;
772 trackIt[i] = 0;
773 }
774
775 Int_t nc = 0; // Total n. of selected particles
776 Int_t nParents = 0; // Selected parents
777 Int_t nTkbles = 0; // Trackable particles
f529e69b 778 if (fProcess != kPyMbDefault &&
779 fProcess != kPyMb &&
6d2ec66d 780 fProcess != kPyMbAtlasTuneMC09 &&
04081a91 781 fProcess != kPyMbWithDirectPhoton &&
f529e69b 782 fProcess != kPyJets &&
8d2cd130 783 fProcess != kPyDirectGamma &&
589380c6 784 fProcess != kPyMbNonDiffr &&
d7de4a67 785 fProcess != kPyMbMSEL1 &&
f529e69b 786 fProcess != kPyW &&
787 fProcess != kPyZ &&
1f0caa6a 788 fProcess != kPyZgamma &&
f529e69b 789 fProcess != kPyCharmppMNRwmi &&
9dfe63b3 790 fProcess != kPyBeautyppMNRwmi &&
64da86aa 791 fProcess != kPyBeautyJets &&
1f0caa6a 792 fProcess != kPyWPWHG &&
793 fProcess != kPyJetsPWHG &&
794 fProcess != kPyCharmPWHG &&
795 fProcess != kPyBeautyPWHG) {
8d2cd130 796
5fa4b20b 797 for (i = 0; i < np; i++) {
8507138f 798 TParticle* iparticle = (TParticle *) fParticles.At(i);
8d2cd130 799 Int_t ks = iparticle->GetStatusCode();
800 kf = CheckPDGCode(iparticle->GetPdgCode());
801// No initial state partons
802 if (ks==21) continue;
803//
804// Heavy Flavor Selection
805//
806 // quark ?
807 kf = TMath::Abs(kf);
808 Int_t kfl = kf;
9ff6c04c 809 // Resonance
f913ec4f 810
9ff6c04c 811 if (kfl > 100000) kfl %= 100000;
183a5ca9 812 if (kfl > 10000) kfl %= 10000;
8d2cd130 813 // meson ?
814 if (kfl > 10) kfl/=100;
815 // baryon
816 if (kfl > 10) kfl/=10;
8d2cd130 817 Int_t ipa = iparticle->GetFirstMother()-1;
818 Int_t kfMo = 0;
f913ec4f 819//
820// Establish mother daughter relation between heavy quarks and mesons
821//
822 if (kf >= fFlavorSelect && kf <= 6) {
823 Int_t idau = iparticle->GetFirstDaughter() - 1;
824 if (idau > -1) {
8507138f 825 TParticle* daughter = (TParticle *) fParticles.At(idau);
f913ec4f 826 Int_t pdgD = daughter->GetPdgCode();
827 if (pdgD == 91 || pdgD == 92) {
828 Int_t jmin = daughter->GetFirstDaughter() - 1;
829 Int_t jmax = daughter->GetLastDaughter() - 1;
2ab330c9 830 for (Int_t jp = jmin; jp <= jmax; jp++)
831 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
f913ec4f 832 } // is string or cluster
833 } // has daughter
834 } // heavy quark
8d2cd130 835
f913ec4f 836
8d2cd130 837 if (ipa > -1) {
8507138f 838 TParticle * mother = (TParticle *) fParticles.At(ipa);
8d2cd130 839 kfMo = TMath::Abs(mother->GetPdgCode());
840 }
f913ec4f 841
8d2cd130 842 // What to keep in Stack?
843 Bool_t flavorOK = kFALSE;
844 Bool_t selectOK = kFALSE;
845 if (fFeedDownOpt) {
32d6ef7d 846 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 847 } else {
32d6ef7d 848 if (kfl > fFlavorSelect) {
849 nc = -1;
850 break;
851 }
852 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 853 }
854 switch (fStackFillOpt) {
06956108 855 case kHeavyFlavor:
8d2cd130 856 case kFlavorSelection:
32d6ef7d 857 selectOK = kTRUE;
858 break;
8d2cd130 859 case kParentSelection:
32d6ef7d 860 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
861 break;
8d2cd130 862 }
863 if (flavorOK && selectOK) {
864//
865// Heavy flavor hadron or quark
866//
867// Kinematic seletion on final state heavy flavor mesons
868 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
869 {
9ff6c04c 870 continue;
8d2cd130 871 }
872 pSelected[i] = 1;
873 if (ParentSelected(kf)) ++nParents; // Update parent count
874// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
875 } else {
876// Kinematic seletion on decay products
877 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 878 && !KinematicSelection(iparticle, 1))
8d2cd130 879 {
9ff6c04c 880 continue;
8d2cd130 881 }
882//
883// Decay products
884// Select if mother was selected and is not tracked
885
886 if (pSelected[ipa] &&
887 !trackIt[ipa] && // mother will be tracked ?
888 kfMo != 5 && // mother is b-quark, don't store fragments
889 kfMo != 4 && // mother is c-quark, don't store fragments
890 kf != 92) // don't store string
891 {
892//
893// Semi-stable or de-selected: diselect decay products:
894//
895//
896 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
897 {
898 Int_t ipF = iparticle->GetFirstDaughter();
899 Int_t ipL = iparticle->GetLastDaughter();
900 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
901 }
902// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
903 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
904 }
905 }
906 if (pSelected[i] == -1) pSelected[i] = 0;
907 if (!pSelected[i]) continue;
908 // Count quarks only if you did not include fragmentation
909 if (fFragmentation && kf <= 10) continue;
9ff6c04c 910
8d2cd130 911 nc++;
912// Decision on tracking
913 trackIt[i] = 0;
914//
915// Track final state particle
916 if (ks == 1) trackIt[i] = 1;
917// Track semi-stable particles
d25cfd65 918 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 919// Track particles selected by process if undecayed.
920 if (fForceDecay == kNoDecay) {
921 if (ParentSelected(kf)) trackIt[i] = 1;
922 } else {
923 if (ParentSelected(kf)) trackIt[i] = 0;
924 }
925 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
926//
927//
928
929 } // particle selection loop
930 if (nc > 0) {
931 for (i = 0; i<np; i++) {
932 if (!pSelected[i]) continue;
8507138f 933 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 934 kf = CheckPDGCode(iparticle->GetPdgCode());
935 Int_t ks = iparticle->GetStatusCode();
936 p[0] = iparticle->Px();
937 p[1] = iparticle->Py();
938 p[2] = iparticle->Pz();
a920faf9 939 p[3] = iparticle->Energy();
940
2590ccf9 941 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
942 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
943 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
944
21391258 945 Float_t tof = fTime + kconv*iparticle->T();
8d2cd130 946 Int_t ipa = iparticle->GetFirstMother()-1;
947 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 948
949 PushTrack(fTrackIt*trackIt[i], iparent, kf,
950 p[0], p[1], p[2], p[3],
951 origin[0], origin[1], origin[2], tof,
952 polar[0], polar[1], polar[2],
953 kPPrimary, nt, 1., ks);
8d2cd130 954 pParent[i] = nt;
dbd64db6 955 KeepTrack(nt);
956 fNprimaries++;
642f15cf 957 } // PushTrack loop
8d2cd130 958 }
959 } else {
960 nc = GenerateMB();
961 } // mb ?
f913ec4f 962
963 GetSubEventTime();
8d2cd130 964
234f6d32 965 delete[] pParent;
966 delete[] pSelected;
967 delete[] trackIt;
8d2cd130 968
969 if (nc > 0) {
970 switch (fCountMode) {
971 case kCountAll:
972 // printf(" Count all \n");
973 jev += nc;
974 break;
975 case kCountParents:
976 // printf(" Count parents \n");
977 jev += nParents;
978 break;
979 case kCountTrackables:
980 // printf(" Count trackable \n");
981 jev += nTkbles;
982 break;
983 }
984 if (jev >= fNpart || fNpart == -1) {
985 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 986
8d2cd130 987 fQ += fPythia->GetVINT(51);
988 fX1 += fPythia->GetVINT(41);
989 fX2 += fPythia->GetVINT(42);
990 fTrialsRun += fTrials;
991 fNev++;
992 MakeHeader();
993 break;
994 }
995 }
996 } // event loop
997 SetHighWaterMark(nt);
998// adjust weight due to kinematic selection
b88f5cea 999// AdjustWeights();
8d2cd130 1000// get cross-section
1001 fXsection=fPythia->GetPARI(1);
1002}
1003
1004Int_t AliGenPythia::GenerateMB()
1005{
1006//
1007// Min Bias selection and other global selections
1008//
1009 Int_t i, kf, nt, iparent;
1010 Int_t nc = 0;
13cca24a 1011 Double_t p[4];
1012 Double_t polar[3] = {0,0,0};
1013 Double_t origin[3] = {0,0,0};
8d2cd130 1014// converts from mm/c to s
1015 const Float_t kconv=0.001/2.999792458e8;
1016
e0e89f40 1017
8507138f 1018 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
aea21c57 1019
5fa4b20b 1020
e0e89f40 1021
8d2cd130 1022 Int_t* pParent = new Int_t[np];
1023 for (i=0; i< np; i++) pParent[i] = -1;
64da86aa 1024 if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
0e9e66f0 1025 && fEtMinJet > 0.) {
bed3df2c 1026 TParticle* jet1 = (TParticle *) fParticles.At(6);
8507138f 1027 TParticle* jet2 = (TParticle *) fParticles.At(7);
bed3df2c 1028
1029 if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
234f6d32 1030 delete [] pParent;
1031 return 0;
1032 }
8d2cd130 1033 }
e0e89f40 1034
40fe59d4 1035
1036 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
ac01c7c6 1037 // implemented primaryly for kPyJets, but extended to any kind of process.
40fe59d4 1038 if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
1039 (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
1040 Bool_t ok = TriggerOnSelectedParticles(np);
43e3b80a 1041
1042 if(!ok) {
1043 delete[] pParent;
1044 return 0;
ec2c406e 1045 }
43e3b80a 1046 }
9dfe63b3 1047
800be8b7 1048 // Check for diffraction
1049 if(fkTuneForDiff) {
c5dfa3e4 1050 if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
800be8b7 1051 if(!CheckDiffraction()) {
1052 delete [] pParent;
1053 return 0;
1054 }
1055 }
1056 }
1057
700b9416 1058 // Check for minimum multiplicity
1059 if (fTriggerMultiplicity > 0) {
1060 Int_t multiplicity = 0;
1061 for (i = 0; i < np; i++) {
8507138f 1062 TParticle * iparticle = (TParticle *) fParticles.At(i);
700b9416 1063
1064 Int_t statusCode = iparticle->GetStatusCode();
1065
1066 // Initial state particle
e296848e 1067 if (statusCode != 1)
700b9416 1068 continue;
38112f3f 1069 // eta cut
700b9416 1070 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1071 continue;
440c2873 1072 //multiplicity check for a given eta range
1073 if ((fTriggerMultiplicityEtaMin != fTriggerMultiplicityEtaMax) &&
1074 (iparticle->Eta() < fTriggerMultiplicityEtaMin || iparticle->Eta() > fTriggerMultiplicityEtaMax))
1075 continue;
38112f3f 1076 // pt cut
1077 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1078 continue;
1079
700b9416 1080 TParticlePDG* pdgPart = iparticle->GetPDG();
1081 if (pdgPart && pdgPart->Charge() == 0)
1082 continue;
1083
1084 ++multiplicity;
1085 }
1086
1087 if (multiplicity < fTriggerMultiplicity) {
1088 delete [] pParent;
1089 return 0;
1090 }
38112f3f 1091 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
700b9416 1092 }
90a236ce 1093
90a236ce 1094
7ce1321b 1095 if (fTriggerParticle) {
1096 Bool_t triggered = kFALSE;
1097 for (i = 0; i < np; i++) {
8507138f 1098 TParticle * iparticle = (TParticle *) fParticles.At(i);
7ce1321b 1099 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 1100 if (kf != fTriggerParticle) continue;
7ce1321b 1101 if (iparticle->Pt() == 0.) continue;
4b5ec27d 1102 if (TMath::Abs(iparticle->Y()) > fTriggerY) continue;
1f0caa6a 1103 if (fTriggerEtaMin == fTriggerEta) {
1104 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1105 } else {
1106 if (iparticle->Eta() > fTriggerEta || iparticle->Eta() < fTriggerEtaMin) continue;
1107 }
2591bd0e 1108 if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
7ce1321b 1109 triggered = kTRUE;
1110 break;
1111 }
234f6d32 1112 if (!triggered) {
1113 delete [] pParent;
1114 return 0;
1115 }
7ce1321b 1116 }
e0e89f40 1117
1118
1119 // Check if there is a ccbar or bbbar pair with at least one of the two
1120 // in fYMin < y < fYMax
2f405d65 1121
9dfe63b3 1122 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
9ccc3d0e 1123 TParticle *partCheck;
1124 TParticle *mother;
e0e89f40 1125 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 1126 Bool_t theChild=kFALSE;
5d76923e 1127 Bool_t triggered=kFALSE;
9ccc3d0e 1128 Float_t y;
1f0caa6a 1129 Int_t pdg, mpdg, mpdgUpperFamily;
1130 for(i = 0; i < np; i++) {
9ccc3d0e 1131 partCheck = (TParticle*)fParticles.At(i);
1132 pdg = partCheck->GetPdgCode();
1133 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1134 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1135 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1136 (partCheck->Energy()-partCheck->Pz()+1.e-13));
d8850d4f 1137 if(fUseYCutHQ && y>fYMinHQ && y<fYMaxHQ) inYcut=kTRUE;
1138 if(!fUseYCutHQ && y>fYMin && y<fYMax) inYcut=kTRUE;
9ccc3d0e 1139 }
5d76923e 1140 if(fTriggerParticle) {
1141 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1142 }
9ccc3d0e 1143 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1144 Int_t mi = partCheck->GetFirstMother() - 1;
1145 if(mi<0) continue;
1146 mother = (TParticle*)fParticles.At(mi);
1147 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 1148 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 1149 if ( ParentSelected(mpdg) ||
1150 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1151 if (KinematicSelection(partCheck,1)) {
1152 theChild=kTRUE;
1153 }
1154 }
1155 }
e0e89f40 1156 }
9ccc3d0e 1157 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
234f6d32 1158 delete[] pParent;
e0e89f40 1159 return 0;
1160 }
9ccc3d0e 1161 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1162 delete[] pParent;
1163 return 0;
1164 }
5d76923e 1165 if(fTriggerParticle && !triggered) { // particle requested is not produced
1166 delete[] pParent;
1167 return 0;
1168 }
9ccc3d0e 1169
e0e89f40 1170 }
aea21c57 1171
0f6ee828 1172 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
e136a92a 1173 if ( (
1f0caa6a 1174 fProcess == kPyWPWHG ||
1175 fProcess == kPyW ||
f529e69b 1176 fProcess == kPyZ ||
1f0caa6a 1177 fProcess == kPyZgamma ||
f529e69b 1178 fProcess == kPyMbDefault ||
1179 fProcess == kPyMb ||
6d2ec66d 1180 fProcess == kPyMbAtlasTuneMC09 ||
04081a91 1181 fProcess == kPyMbWithDirectPhoton ||
f529e69b 1182 fProcess == kPyMbNonDiffr)
0f6ee828 1183 && (fCutOnChild == 1) ) {
1184 if ( !CheckKinematicsOnChild() ) {
234f6d32 1185 delete[] pParent;
0f6ee828 1186 return 0;
1187 }
aea21c57 1188 }
1f0caa6a 1189
aea21c57 1190
f913ec4f 1191 for (i = 0; i < np; i++) {
8d2cd130 1192 Int_t trackIt = 0;
8507138f 1193 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 1194 kf = CheckPDGCode(iparticle->GetPdgCode());
1195 Int_t ks = iparticle->GetStatusCode();
1196 Int_t km = iparticle->GetFirstMother();
06956108 1197 if (
1198 (((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
64da86aa 1199 ((fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) && ks == 21 && km == 0 && i>1)
06956108 1200 )
1201 {
1202 nc++;
8d2cd130 1203 if (ks == 1) trackIt = 1;
1204 Int_t ipa = iparticle->GetFirstMother()-1;
1205
1206 iparent = (ipa > -1) ? pParent[ipa] : -1;
1207
1208//
1209// store track information
1210 p[0] = iparticle->Px();
1211 p[1] = iparticle->Py();
1212 p[2] = iparticle->Pz();
a920faf9 1213 p[3] = iparticle->Energy();
1406f599 1214
a920faf9 1215
2590ccf9 1216 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1217 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1218 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1219
21391258 1220 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
a920faf9 1221
1222 PushTrack(fTrackIt*trackIt, iparent, kf,
1223 p[0], p[1], p[2], p[3],
1224 origin[0], origin[1], origin[2], tof,
1225 polar[0], polar[1], polar[2],
1226 kPPrimary, nt, 1., ks);
dbd64db6 1227 fNprimaries++;
8d2cd130 1228 KeepTrack(nt);
1229 pParent[i] = nt;
f913ec4f 1230 SetHighWaterMark(nt);
1231
8d2cd130 1232 } // select particle
1233 } // particle loop
1234
234f6d32 1235 delete[] pParent;
e0e89f40 1236
f913ec4f 1237 return 1;
8d2cd130 1238}
1239
1240
1241void AliGenPythia::FinishRun()
1242{
1243// Print x-section summary
1244 fPythia->Pystat(1);
2779fc64 1245
1246 if (fNev > 0.) {
1247 fQ /= fNev;
1248 fX1 /= fNev;
1249 fX2 /= fNev;
1250 }
1251
8d2cd130 1252 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1253 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1254}
1255
7184e472 1256void AliGenPythia::AdjustWeights() const
8d2cd130 1257{
1258// Adjust the weights after generation of all events
1259//
e2bddf81 1260 if (gAlice) {
1261 TParticle *part;
1262 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1263 for (Int_t i=0; i<ntrack; i++) {
1264 part= gAlice->GetMCApp()->Particle(i);
1265 part->SetWeight(part->GetWeight()*fKineBias);
1266 }
8d2cd130 1267 }
1268}
1269
20e47f08 1270void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
8d2cd130 1271{
1272// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1273
fb355e51 1274 fAProjectile = a1;
1275 fATarget = a2;
1276 fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
1277 fUseNuclearPDF = kTRUE;
1278 fSetNuclei = kTRUE;
8d2cd130 1279}
1280
1281
1282void AliGenPythia::MakeHeader()
1283{
7184e472 1284//
1285// Make header for the simulated event
1286//
183a5ca9 1287 if (gAlice) {
1288 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1289 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1290 }
1291
8d2cd130 1292// Builds the event header, to be called after each event
e5c87a3d 1293 if (fHeader) delete fHeader;
1294 fHeader = new AliGenPythiaEventHeader("Pythia");
d8850d4f 1295 fHeader->SetTitle(GetTitle());
8d2cd130 1296//
1297// Event type
800be8b7 1298 if(fProcDiff>0){
1299 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1300 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1301 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1302 }
1303 else
e5c87a3d 1304 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1305//
1306// Number of trials
e5c87a3d 1307 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1308//
1309// Event Vertex
d25cfd65 1310 fHeader->SetPrimaryVertex(fVertex);
21391258 1311 fHeader->SetInteractionTime(fTime+fEventTime);
dbd64db6 1312//
1313// Number of primaries
1314 fHeader->SetNProduced(fNprimaries);
8d2cd130 1315//
a0027228 1316// Event weight
1317 fHeader->SetEventWeight(fPythia->GetVINT(97));
1318//
8d2cd130 1319// Jets that have triggered
f913ec4f 1320
9dfe63b3 1321 //Need to store jets for b-jet studies too!
64da86aa 1322 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
8d2cd130 1323 {
1324 Int_t ntrig, njet;
1325 Float_t jets[4][10];
1326 GetJets(njet, ntrig, jets);
9ff6c04c 1327
8d2cd130 1328
1329 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1330 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1331 jets[3][i]);
1332 }
1333 }
5fa4b20b 1334//
1335// Copy relevant information from external header, if present.
1336//
1337 Float_t uqJet[4];
1338
1339 if (fRL) {
1340 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1341 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1342 {
1343 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1344
1345
1346 exHeader->TriggerJet(i, uqJet);
1347 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1348 }
1349 }
1350//
1351// Store quenching parameters
1352//
1353 if (fQuench){
b6262a45 1354 Double_t z[4] = {0.};
1355 Double_t xp = 0.;
1356 Double_t yp = 0.;
1357
7c21f297 1358 if (fQuench == 1) {
1359 // Pythia::Quench()
1360 fPythia->GetQuenchingParameters(xp, yp, z);
1bab4b79 1361 } else if (fQuench == 2){
7c21f297 1362 // Pyquen
1363 Double_t r1 = PARIMP.rb1;
1364 Double_t r2 = PARIMP.rb2;
1365 Double_t b = PARIMP.b1;
1366 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1367 Double_t phi = PARIMP.psib1;
1368 xp = r * TMath::Cos(phi);
1369 yp = r * TMath::Sin(phi);
1370
1bab4b79 1371 } else if (fQuench == 4) {
1372 // QPythia
5831e75f 1373 Double_t xy[2];
e9719084 1374 Double_t i0i1[2];
1375 AliFastGlauber::Instance()->GetSavedXY(xy);
1376 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
5831e75f 1377 xp = xy[0];
1378 yp = xy[1];
e6fe9b82 1379 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
7c21f297 1380 }
1bab4b79 1381
7c21f297 1382 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1383 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1bab4b79 1384 }
beac474c 1385//
39ea8b7f 1386// Store pt^hard and cross-section
beac474c 1387 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
39ea8b7f 1388 ((AliGenPythiaEventHeader*) fHeader)->SetXsection(fPythia->GetPARI(1));
64da86aa 1389
1390//
1391// Store Event Weight
1cdfb96f 1392 ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7)*fPythia->GetPARI(10));
1393 // PARI(7) is 1 or -1, for weighted generation with accept/reject, e.g. POWHEG
1394 // PARI(10) is a weight associated with reweighted generation, using Pyevwt
5fa4b20b 1395//
cf57b268 1396// Pass header
5fa4b20b 1397//
cf57b268 1398 AddHeader(fHeader);
4c4eac97 1399 fHeader = 0x0;
8d2cd130 1400}
cf57b268 1401
c2fc9d31 1402Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
8d2cd130 1403{
1404// Check the kinematic trigger condition
1405//
1406 Double_t eta[2];
1407 eta[0] = jet1->Eta();
1408 eta[1] = jet2->Eta();
1409 Double_t phi[2];
1410 phi[0] = jet1->Phi();
1411 phi[1] = jet2->Phi();
1412 Int_t pdg[2];
1413 pdg[0] = jet1->GetPdgCode();
1414 pdg[1] = jet2->GetPdgCode();
1415 Bool_t triggered = kFALSE;
1416
64da86aa 1417 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) {
8d2cd130 1418 Int_t njets = 0;
1419 Int_t ntrig = 0;
1420 Float_t jets[4][10];
1421//
1422// Use Pythia clustering on parton level to determine jet axis
1423//
1424 GetJets(njets, ntrig, jets);
1425
76d6ba9a 1426 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1427//
1428 } else {
1429 Int_t ij = 0;
1430 Int_t ig = 1;
1431 if (pdg[0] == kGamma) {
1432 ij = 1;
1433 ig = 0;
1434 }
1435 //Check eta range first...
1436 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1437 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1438 {
1439 //Eta is okay, now check phi range
1440 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1441 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1442 {
1443 triggered = kTRUE;
1444 }
1445 }
1446 }
1447 return triggered;
1448}
aea21c57 1449
1450
aea21c57 1451
7184e472 1452Bool_t AliGenPythia::CheckKinematicsOnChild(){
1453//
1454//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1455//
aea21c57 1456 Bool_t checking = kFALSE;
1457 Int_t j, kcode, ks, km;
1458 Int_t nPartAcc = 0; //number of particles in the acceptance range
1459 Int_t numberOfAcceptedParticles = 1;
1460 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1461 Int_t npart = fParticles.GetEntriesFast();
aea21c57 1462
0f6ee828 1463 for (j = 0; j<npart; j++) {
8507138f 1464 TParticle * jparticle = (TParticle *) fParticles.At(j);
aea21c57 1465 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1466 ks = jparticle->GetStatusCode();
1467 km = jparticle->GetFirstMother();
1468
1469 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1470 nPartAcc++;
1471 }
0f6ee828 1472 if( numberOfAcceptedParticles <= nPartAcc){
1473 checking = kTRUE;
1474 break;
1475 }
aea21c57 1476 }
0f6ee828 1477
aea21c57 1478 return checking;
aea21c57 1479}
1480
5fa4b20b 1481void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1482{
1058d9df 1483 //
1484 // Load event into Pythia Common Block
1485 //
1486
1487 Int_t npart = stack -> GetNprimary();
1488 Int_t n0 = 0;
1489
1490 if (!flag) {
1491 (fPythia->GetPyjets())->N = npart;
1492 } else {
1493 n0 = (fPythia->GetPyjets())->N;
1494 (fPythia->GetPyjets())->N = n0 + npart;
1495 }
1496
1497
1498 for (Int_t part = 0; part < npart; part++) {
1499 TParticle *mPart = stack->Particle(part);
32d6ef7d 1500
1058d9df 1501 Int_t kf = mPart->GetPdgCode();
1502 Int_t ks = mPart->GetStatusCode();
1503 Int_t idf = mPart->GetFirstDaughter();
1504 Int_t idl = mPart->GetLastDaughter();
1505
1506 if (reHadr) {
1507 if (ks == 11 || ks == 12) {
1508 ks -= 10;
1509 idf = -1;
1510 idl = -1;
1511 }
32d6ef7d 1512 }
1513
1058d9df 1514 Float_t px = mPart->Px();
1515 Float_t py = mPart->Py();
1516 Float_t pz = mPart->Pz();
1517 Float_t e = mPart->Energy();
1518 Float_t m = mPart->GetCalcMass();
32d6ef7d 1519
1058d9df 1520
1521 (fPythia->GetPyjets())->P[0][part+n0] = px;
1522 (fPythia->GetPyjets())->P[1][part+n0] = py;
1523 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1524 (fPythia->GetPyjets())->P[3][part+n0] = e;
1525 (fPythia->GetPyjets())->P[4][part+n0] = m;
1526
1527 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1528 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1529 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1530 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1531 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1532 }
1533}
1534
c2fc9d31 1535void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1058d9df 1536{
1537 //
1538 // Load event into Pythia Common Block
1539 //
1540
1541 Int_t npart = stack -> GetEntries();
1542 Int_t n0 = 0;
1543
1544 if (!flag) {
1545 (fPythia->GetPyjets())->N = npart;
1546 } else {
1547 n0 = (fPythia->GetPyjets())->N;
1548 (fPythia->GetPyjets())->N = n0 + npart;
1549 }
1550
1551
1552 for (Int_t part = 0; part < npart; part++) {
1553 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
92847124 1554 if (!mPart) continue;
1555
1058d9df 1556 Int_t kf = mPart->GetPdgCode();
1557 Int_t ks = mPart->GetStatusCode();
1558 Int_t idf = mPart->GetFirstDaughter();
1559 Int_t idl = mPart->GetLastDaughter();
1560
1561 if (reHadr) {
92847124 1562 if (ks == 11 || ks == 12) {
1563 ks -= 10;
1564 idf = -1;
1565 idl = -1;
1566 }
8d2cd130 1567 }
1058d9df 1568
1569 Float_t px = mPart->Px();
1570 Float_t py = mPart->Py();
1571 Float_t pz = mPart->Pz();
1572 Float_t e = mPart->Energy();
1573 Float_t m = mPart->GetCalcMass();
1574
1575
1576 (fPythia->GetPyjets())->P[0][part+n0] = px;
1577 (fPythia->GetPyjets())->P[1][part+n0] = py;
1578 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1579 (fPythia->GetPyjets())->P[3][part+n0] = e;
1580 (fPythia->GetPyjets())->P[4][part+n0] = m;
1581
1582 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1583 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1584 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1585 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1586 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1587 }
8d2cd130 1588}
1589
5fa4b20b 1590
014a9521 1591void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1592{
1593//
1594// Calls the Pythia jet finding algorithm to find jets in the current event
1595//
1596//
8d2cd130 1597//
1598// Save jets
1599 Int_t n = fPythia->GetN();
1600
1601//
1602// Run Jet Finder
1603 fPythia->Pycell(njets);
1604 Int_t i;
1605 for (i = 0; i < njets; i++) {
1606 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1607 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1608 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1609 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1610
1611 jets[0][i] = px;
1612 jets[1][i] = py;
1613 jets[2][i] = pz;
1614 jets[3][i] = e;
1615 }
1616}
1617
1618
1619
1620void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1621{
1622//
1623// Calls the Pythia clustering algorithm to find jets in the current event
1624//
1625 Int_t n = fPythia->GetN();
1626 nJets = 0;
1627 nJetsTrig = 0;
1628 if (fJetReconstruction == kCluster) {
1629//
1630// Configure cluster algorithm
1631//
1632 fPythia->SetPARU(43, 2.);
1633 fPythia->SetMSTU(41, 1);
1634//
1635// Call cluster algorithm
1636//
1637 fPythia->Pyclus(nJets);
1638//
1639// Loading jets from common block
1640//
1641 } else {
592f8307 1642
8d2cd130 1643//
1644// Run Jet Finder
1645 fPythia->Pycell(nJets);
1646 }
1647
1648 Int_t i;
1649 for (i = 0; i < nJets; i++) {
1650 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1651 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1652 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1653 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1654 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1655 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1656 Float_t theta = TMath::ATan2(pt,pz);
1657 Float_t et = e * TMath::Sin(theta);
1658 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1659 if (
1660 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1661 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1662 et > fEtMinJet && et < fEtMaxJet
1663 )
1664 {
1665 jets[0][nJetsTrig] = px;
1666 jets[1][nJetsTrig] = py;
1667 jets[2][nJetsTrig] = pz;
1668 jets[3][nJetsTrig] = e;
1669 nJetsTrig++;
5fa4b20b 1670// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1671 } else {
1672// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1673 }
1674 }
1675}
1676
f913ec4f 1677void AliGenPythia::GetSubEventTime()
1678{
1679 // Calculates time of the next subevent
9d7108a7 1680 fEventTime = 0.;
1681 if (fEventsTime) {
1682 TArrayF &array = *fEventsTime;
1683 fEventTime = array[fCurSubEvent++];
1684 }
1685 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1686 return;
f913ec4f 1687}
8d2cd130 1688
e81f2bac 1689Bool_t AliGenPythia::IsInBarrel(Float_t eta) const
40fe59d4 1690{
1691 // Is particle in Central Barrel acceptance?
1692 // etamin=-etamax
1693 if( eta < fTriggerEta )
1694 return kTRUE;
1695 else
1696 return kFALSE;
1697}
1698
e81f2bac 1699Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
ec2c406e 1700{
1701 // Is particle in EMCAL acceptance?
ec2c406e 1702 // phi in degrees, etamin=-etamax
9fd8e520 1703 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1704 eta < fEMCALEta )
ec2c406e 1705 return kTRUE;
1706 else
1707 return kFALSE;
1708}
1709
e81f2bac 1710Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
ec2c406e 1711{
1712 // Is particle in PHOS acceptance?
1713 // Acceptance slightly larger considered.
1714 // phi in degrees, etamin=-etamax
40fe59d4 1715 // iparticle is the index of the particle to be checked, for PHOS rotation case
1716
9fd8e520 1717 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1718 eta < fPHOSEta )
ec2c406e 1719 return kTRUE;
1720 else
40fe59d4 1721 {
1722 if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1723
ec2c406e 1724 return kFALSE;
40fe59d4 1725 }
ec2c406e 1726}
1727
40fe59d4 1728void AliGenPythia::RotatePhi(Bool_t& okdd)
90a236ce 1729{
40fe59d4 1730 //Rotate event in phi to enhance events in PHOS acceptance
1731
1732 if(fPHOSRotateCandidate < 0) return ;
1733
90a236ce 1734 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1735 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1736 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
800be8b7 1737 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
90a236ce 1738
1739 //calculate deltaphi
40fe59d4 1740 TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
90a236ce 1741 Double_t phphi = ph->Phi();
1742 Double_t deltaphi = phiPHOS - phphi;
40fe59d4 1743
90a236ce 1744
1745
1746 //loop for all particles and produce the phi rotation
8507138f 1747 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
90a236ce 1748 Double_t oldphi, newphi;
777e69b0 1749 Double_t newVx, newVy, r, vZ, time;
1750 Double_t newPx, newPy, pt, pz, e;
90a236ce 1751 for(Int_t i=0; i< np; i++) {
40fe59d4 1752 TParticle* iparticle = (TParticle *) fParticles.At(i);
1753 oldphi = iparticle->Phi();
1754 newphi = oldphi + deltaphi;
1755 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1756 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1757
1758 r = iparticle->R();
1759 newVx = r * TMath::Cos(newphi);
1760 newVy = r * TMath::Sin(newphi);
1761 vZ = iparticle->Vz(); // don't transform
1762 time = iparticle->T(); // don't transform
1763
1764 pt = iparticle->Pt();
1765 newPx = pt * TMath::Cos(newphi);
1766 newPy = pt * TMath::Sin(newphi);
1767 pz = iparticle->Pz(); // don't transform
1768 e = iparticle->Energy(); // don't transform
1769
1770 // apply rotation
1771 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1772 iparticle->SetMomentum(newPx, newPy, pz, e);
1773
90a236ce 1774 } //end particle loop
1775
40fe59d4 1776 // now let's check that we put correctly the candidate photon in PHOS
1777 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1778 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1779 if(IsInPHOS(phi,eta,-1))
1780 okdd = kTRUE;
1781
1782 // reset the value for next event
1783 fPHOSRotateCandidate = -1;
1784
90a236ce 1785}
ec2c406e 1786
1787
800be8b7 1788Bool_t AliGenPythia::CheckDiffraction()
1789{
1790 // use this method only with Perugia-0 tune!
1791
1792 // printf("AAA\n");
1793
1794 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1795
1796 Int_t iPart1=-1;
1797 Int_t iPart2=-1;
1798
1799 Double_t y1 = 1e10;
1800 Double_t y2 = -1e10;
1801
1802 const Int_t kNstable=20;
1803 const Int_t pdgStable[20] = {
1804 22, // Photon
1805 11, // Electron
1806 12, // Electron Neutrino
1807 13, // Muon
1808 14, // Muon Neutrino
1809 15, // Tau
1810 16, // Tau Neutrino
1811 211, // Pion
1812 321, // Kaon
1813 311, // K0
1814 130, // K0s
1815 310, // K0l
1816 2212, // Proton
1817 2112, // Neutron
1818 3122, // Lambda_0
1819 3112, // Sigma Minus
1820 3222, // Sigma Plus
1821 3312, // Xsi Minus
1822 3322, // Xsi0
1823 3334 // Omega
1824 };
1825
1826 for (Int_t i = 0; i < np; i++) {
1827 TParticle * part = (TParticle *) fParticles.At(i);
1828
1829 Int_t statusCode = part->GetStatusCode();
1830
1831 // Initial state particle
1832 if (statusCode != 1)
1833 continue;
1834
1835 Int_t pdg = TMath::Abs(part->GetPdgCode());
1836 Bool_t isStable = kFALSE;
1837 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1838 if (pdg == pdgStable[i1]) {
1839 isStable = kTRUE;
1840 break;
1841 }
1842 }
1843 if(!isStable)
1844 continue;
1845
1846 Double_t y = part->Y();
1847
1848 if (y < y1)
1849 {
1850 y1 = y;
1851 iPart1 = i;
1852 }
1853 if (y > y2)
1854 {
1855 y2 = y;
1856 iPart2 = i;
1857 }
1858 }
1859
1860 if(iPart1<0 || iPart2<0) return kFALSE;
1861
1862 y1=TMath::Abs(y1);
1863 y2=TMath::Abs(y2);
1864
1865 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1866 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1867
1868 Int_t pdg1 = part1->GetPdgCode();
1869 Int_t pdg2 = part2->GetPdgCode();
1870
1871
1872 Int_t iPart = -1;
1873 if (pdg1 == 2212 && pdg2 == 2212)
1874 {
1875 if(y1 > y2)
1876 iPart = iPart1;
1877 else if(y1 < y2)
1878 iPart = iPart2;
1879 else {
1880 iPart = iPart1;
1881 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1882 }
1883 }
1884 else if (pdg1 == 2212)
1885 iPart = iPart1;
1886 else if (pdg2 == 2212)
1887 iPart = iPart2;
1888
1889
1890
1891
1892
1893 Double_t M=-1.;
1894 if(iPart>0) {
1895 TParticle * part = (TParticle *) fParticles.At(iPart);
1896 Double_t E= part->Energy();
1897 Double_t P= part->P();
4f813e90 1898 Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
1899 if(M2<0) return kFALSE;
1900 M= TMath::Sqrt(M2);
800be8b7 1901 }
1902
c5dfa3e4 1903 Double_t Mmin, Mmax, wSD, wDD, wND;
1904 if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1905
1906 if(M>-1 && M<Mmin) return kFALSE;
1907 if(M>Mmax) M=-1;
1908
1909 Int_t procType=fPythia->GetMSTI(1);
1910 Int_t proc0=2;
1911 if(procType== 94) proc0=1;
1912 if(procType== 92 || procType== 93) proc0=0;
1913
1914 Int_t proc=2;
1915 if(M>0) proc=0;
1916 else if(proc0==1) proc=1;
1917
1918 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1919 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
800be8b7 1920
1921
1922 // if(proc==1 || proc==2) return kFALSE;
1923
c5dfa3e4 1924 if(proc!=0) {
1925 if(proc0!=0) fProcDiff = procType;
1926 else fProcDiff = 95;
1927 return kTRUE;
1928 }
1929
1930 if(wSD<0) AliError("wSD<0 ! \n");
1931
1932 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1933
1934 // printf("iPart = %d\n", iPart);
1935
1936 if(iPart==iPart1) fProcDiff=93;
1937 else if(iPart==iPart2) fProcDiff=92;
1938 else {
1939 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1940
800be8b7 1941 }
1942
c5dfa3e4 1943 return kTRUE;
1944}
1945
1946
1947
1948Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1949 Double_t &wSD, Double_t &wDD, Double_t &wND)
1950{
1951
1952 // 900 GeV
1953 if(TMath::Abs(fEnergyCMS-900)<1 ){
1954
1955const Int_t nbin=400;
1956Double_t bin[]={
19571.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
19584.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
19597.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
196010.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
196113.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
196215.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
196318.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
196421.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
196524.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
196627.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
196730.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
196833.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
196936.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
197039.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
197142.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
197245.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
197348.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
197451.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
197554.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
197657.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
197760.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
197863.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
197966.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
198069.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
198172.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
198275.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
198378.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
198481.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
198584.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
198687.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
198790.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
198893.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
198996.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
199099.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1991102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1992105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1993108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1994111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1995114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1996117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1997120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1998123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1999126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2000129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2001132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2002135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2003138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2004141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2005144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2006147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2007150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2008153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2009156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2010159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2011162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2012165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2013168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2014171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2015174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2016177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2017180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2018183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2019186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2020189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2021192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2022195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2023198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2024Double_t w[]={
20251.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
20260.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
20270.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
20280.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
20290.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
20300.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
20310.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
20320.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
20330.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
20340.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
20350.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
20360.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
20370.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
20380.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
20390.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
20400.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
20410.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
20420.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
20430.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
20440.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
20450.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
20460.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
20470.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
20480.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
20490.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
20500.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
20510.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
20520.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
20530.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
20540.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
20550.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
20560.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
20570.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
20580.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
20590.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
20600.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
20610.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
20620.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
20630.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
20640.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
20650.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
20660.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
20670.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
20680.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
20690.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
20700.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
20710.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
20720.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
20730.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
20740.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
20750.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
20760.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
20770.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
20780.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
20790.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
20800.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
20810.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
20820.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
20830.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
20840.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
20850.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
20860.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
20870.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
20880.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
20890.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
20900.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
20910.386112, 0.364314, 0.434375, 0.334629};
2092wDD = 0.379611;
2093wND = 0.496961;
2094wSD = -1;
2095
2096 Mmin = bin[0];
2097 Mmax = bin[nbin];
2098 if(M<Mmin || M>Mmax) return kTRUE;
2099
7873275c 2100 Int_t ibin=nbin-1;
93a60586 2101 for(Int_t i=1; i<=nbin; i++)
2ff57420 2102 if(M<=bin[i]) {
2103 ibin=i-1;
2104 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
800be8b7 2105 break;
2106 }
c5dfa3e4 2107 wSD=w[ibin];
2108 return kTRUE;
2109 }
2110 else if(TMath::Abs(fEnergyCMS-2760)<1 ){
800be8b7 2111
c5dfa3e4 2112const Int_t nbin=400;
2113Double_t bin[]={
21141.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
21154.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
21167.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
211710.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
211813.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
211915.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
212018.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
212121.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
212224.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
212327.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
212430.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
212533.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
212636.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
212739.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
212842.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
212945.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
213048.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
213151.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
213254.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
213357.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
213460.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
213563.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
213666.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
213769.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
213872.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
213975.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
214078.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
214181.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
214284.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
214387.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
214490.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
214593.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
214696.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
214799.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2148102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2149105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2150108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2151111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2152114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2153117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2154120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2155123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2156126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2157129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2158132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2159135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2160138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2161141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2162144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2163147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2164150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2165153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2166156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2167159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2168162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2169165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2170168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2171171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2172174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2173177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2174180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2175183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2176186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2177189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2178192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2179195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2180198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2181Double_t w[]={
21821.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
21830.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
21840.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
21850.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
21860.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
21870.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
21880.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
21890.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
21900.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
21910.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
21920.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
21930.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
21940.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
21950.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
21960.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
21970.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
21980.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
21990.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
22000.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
22010.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
22020.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
22030.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
22040.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
22050.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
22060.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
22070.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
22080.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
22090.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
22100.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
22110.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
22120.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
22130.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
22140.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
22150.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
22160.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
22170.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
22180.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
22190.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
22200.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
22210.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
22220.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
22230.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
22240.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
22250.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
22260.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
22270.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
22280.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
22290.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
22300.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
22310.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
22320.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
22330.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
22340.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
22350.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
22360.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
22370.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
22380.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
22390.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
22400.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
22410.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
22420.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
22430.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
22440.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
22450.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
22460.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
22470.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
22480.201712, 0.242225, 0.255565, 0.258738};
2249wDD = 0.512813;
2250wND = 0.518820;
2251wSD = -1;
2252
2253 Mmin = bin[0];
2254 Mmax = bin[nbin];
2255 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2256
c5dfa3e4 2257 Int_t ibin=nbin-1;
2258 for(Int_t i=1; i<=nbin; i++)
2259 if(M<=bin[i]) {
2260 ibin=i-1;
2261 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2262 break;
2263 }
2264 wSD=w[ibin];
2265 return kTRUE;
2266 }
800be8b7 2267
800be8b7 2268
c5dfa3e4 2269 else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2270const Int_t nbin=400;
2271Double_t bin[]={
22721.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
22734.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
22747.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
227510.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
227613.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
227715.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
227818.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
227921.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
228024.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
228127.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
228230.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
228333.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
228436.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
228539.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
228642.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
228745.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
228848.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
228951.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
229054.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
229157.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
229260.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
229363.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
229466.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
229569.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
229672.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
229775.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
229878.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
229981.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
230084.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
230187.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
230290.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
230393.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
230496.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
230599.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2306102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2307105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2308108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2309111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2310114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2311117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2312120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2313123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2314126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2315129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2316132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2317135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2318138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2319141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2320144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2321147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2322150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2323153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2324156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2325159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2326162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2327165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2328168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2329171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2330174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2331177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2332180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2333183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2334186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2335189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2336192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2337195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2338198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2339Double_t w[]={
23401.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
23410.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
23420.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
23430.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
23440.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
23450.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
23460.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
23470.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
23480.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
23490.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
23500.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
23510.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
23520.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
23530.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
23540.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
23550.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
23560.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
23570.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
23580.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
23590.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
23600.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
23610.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
23620.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
23630.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
23640.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
23650.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
23660.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
23670.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
23680.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
23690.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
23700.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
23710.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
23720.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
23730.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
23740.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
23750.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
23760.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
23770.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
23780.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
23790.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
23800.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
23810.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
23820.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
23830.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
23840.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
23850.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
23860.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
23870.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
23880.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
23890.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
23900.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
23910.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
23920.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
23930.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
23940.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
23950.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
23960.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
23970.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
23980.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
23990.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
24000.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
24010.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
24020.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
24030.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
24040.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
24050.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
24060.175006, 0.223482, 0.202706, 0.218108};
2407wDD = 0.207705;
2408wND = 0.289628;
2409wSD = -1;
2410
2411 Mmin = bin[0];
2412 Mmax = bin[nbin];
2413
2414 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2415
c5dfa3e4 2416 Int_t ibin=nbin-1;
2417 for(Int_t i=1; i<=nbin; i++)
2418 if(M<=bin[i]) {
2419 ibin=i-1;
2420 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2421 break;
2422 }
2423 wSD=w[ibin];
800be8b7 2424 return kTRUE;
c5dfa3e4 2425 }
2426
2427 return kFALSE;
800be8b7 2428}
06956108 2429
2430Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2431{
2432// Check if this is a heavy flavor decay product
2433 TParticle * part = (TParticle *) fParticles.At(ipart);
2434 Int_t mpdg = TMath::Abs(part->GetPdgCode());
2435 Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2436 //
2437 // Light hadron
2438 if (mfl >= 4 && mfl < 6) return kTRUE;
2439
2440 Int_t imo = part->GetFirstMother()-1;
2441 TParticle* pm = part;
2442 //
2443 // Heavy flavor hadron produced by generator
2444 while (imo > -1) {
2445 pm = (TParticle*)fParticles.At(imo);
2446 mpdg = TMath::Abs(pm->GetPdgCode());
2447 mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2448 if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2449 imo = pm->GetFirstMother()-1;
2450 }
2451 return kFALSE;
2452}
40fe59d4 2453
e81f2bac 2454Bool_t AliGenPythia::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
40fe59d4 2455{
2456 // check the eta/phi correspond to the detectors acceptance
2457 // iparticle is the index of the particle to be checked, for PHOS rotation case
2458 if (fCheckPHOS && IsInPHOS (phi,eta,iparticle)) return kTRUE;
2459 else if(fCheckEMCAL && IsInEMCAL (phi,eta)) return kTRUE;
2460 else if(fCheckBarrel && IsInBarrel( eta)) return kTRUE;
2461 else return kFALSE;
2462}
2463
e81f2bac 2464Bool_t AliGenPythia::TriggerOnSelectedParticles(Int_t np)
40fe59d4 2465{
2466 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2467 // implemented primaryly for kPyJets, but extended to any kind of process.
2468
2469 //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2470 // fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
2471
2472 Bool_t ok = kFALSE;
2473 for (Int_t i=0; i< np; i++) {
2474
2475 TParticle* iparticle = (TParticle *) fParticles.At(i);
2476
2477 Int_t pdg = iparticle->GetPdgCode();
2478 Int_t status = iparticle->GetStatusCode();
2479 Int_t imother = iparticle->GetFirstMother() - 1;
2480
2481 TParticle* pmother = 0x0;
2482 Int_t momStatus = -1;
2483 Int_t momPdg = -1;
2484 if(imother > 0 ){
2485 pmother = (TParticle *) fParticles.At(imother);
2486 momStatus = pmother->GetStatusCode();
2487 momPdg = pmother->GetPdgCode();
2488 }
2489
2490 ok = kFALSE;
2491
2492 //
2493 // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2494 //
2495 // Hadron
2496 if (fHadronInCalo && status == 1)
2497 {
2498 if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
2499 // (in case neutral mesons were declared stable)
2500 ok = kTRUE;
2501 }
2502 //Electron
2503 else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2504 {
2505 ok = kTRUE;
2506 }
2507 //Fragmentation photon
2508 else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2509 {
2510 if(momStatus != 11) ok = kTRUE ; // No photon from hadron decay
2511 }
2512 // Decay photon
2513 else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2514 {
2515 if( momStatus == 11)
2516 {
2517 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2518 // pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2519 ok = kTRUE ; // photon from hadron decay
2520
2521 //In case only decays from pi0 or eta requested
2522 if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2523 if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2524 }
2525
2526 }
2527 // Pi0 or Eta particle
2528 else if ((fPi0InCalo || fEtaInCalo))
2529 {
2530 if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2531
2532 if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2533 {
2534 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2535 ok = kTRUE;
2536 }
2537 else if (fEtaInCalo && pdg == 221)
2538 {
2539 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2540 ok = kTRUE;
2541 }
2542
2543 }// pi0 or eta
2544
2545 //
2546 // Check that the selected particle is in the calorimeter acceptance
2547 //
2548 if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2549 {
2550 //Just check if the selected particle falls in the acceptance
2551 if(!fForceNeutralMeson2PhotonDecay )
2552 {
2553 //printf("\t Check acceptance! \n");
2554 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2555 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2556
2557 if(CheckDetectorAcceptance(phi,eta,i))
2558 {
2559 ok =kTRUE;
2560 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));
2561 //printf("\t Accept \n");
2562 break;
2563 }
2564 else ok = kFALSE;
2565 }
2566 //Mesons have several decay modes, select only those decaying into 2 photons
2567 else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2568 {
2569 // In case we want the pi0/eta trigger,
2570 // check the decay mode (2 photons)
2571
2572 //printf("\t Force decay 2 gamma\n");
2573
2574 Int_t ndaughters = iparticle->GetNDaughters();
2575 if(ndaughters != 2){
2576 ok=kFALSE;
2577 continue;
2578 }
2579
2580 TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2581 TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2582 if(!d1 || !d2) {
2583 ok=kFALSE;
2584 continue;
2585 }
2586
2587 //iparticle->Print();
2588 //d1->Print();
2589 //d2->Print();
2590
2591 Int_t pdgD1 = d1->GetPdgCode();
2592 Int_t pdgD2 = d2->GetPdgCode();
2593 //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2594 //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2595
2596 if(pdgD1 != 22 || pdgD2 != 22){
2597 ok = kFALSE;
2598 continue;
2599 }
2600
2601 //printf("\t accept decay\n");
2602
2603 //Trigger on the meson, not on the daughter
2604 if(!fDecayPhotonInCalo){
2605
2606 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2607 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2608
2609 if(CheckDetectorAcceptance(phi,eta,i))
2610 {
2611 //printf("\t Accept meson pdg %d\n",pdg);
2612 ok =kTRUE;
2613 AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2614 break;
2615 } else {
2616 ok=kFALSE;
2617 continue;
2618 }
2619 }
2620
2621 //printf("Check daughters acceptance\n");
2622
2623 //Trigger on the meson daughters
2624 //Photon 1
2625 Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2626 Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
2627 if(d1->Pt() > fTriggerParticleMinPt)
2628 {
2629 //printf("\t Check acceptance photon 1! \n");
2630 if(CheckDetectorAcceptance(phi,eta,i))
2631 {
2632 //printf("\t Accept Photon 1\n");
2633 ok =kTRUE;
2634 AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2635 break;
2636 }
2637 else ok = kFALSE;
2638 } // pt cut
2639 else ok = kFALSE;
2640
2641 //Photon 2
2642 phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2643 eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2644
2645 if(d2->Pt() > fTriggerParticleMinPt)
2646 {
2647 //printf("\t Check acceptance photon 2! \n");
2648 if(CheckDetectorAcceptance(phi,eta,i))
2649 {
2650 //printf("\t Accept Photon 2\n");
2651 ok =kTRUE;
2652 AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2653 break;
2654 }
2655 else ok = kFALSE;
2656 } // pt cut
2657 else ok = kFALSE;
2658 } // force 2 photon daughters in pi0/eta decays
2659 else ok = kFALSE;
2660 } else ok = kFALSE; // check acceptance
2661 } // primary loop
2662
2663 //
2664 // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2665 // A particle passing all trigger conditions except phi position in PHOS, is used as reference
2666 //
2667 if(fCheckPHOSeta)
2668 {
2669 RotatePhi(ok);
2670 }
2671
2672 return ok;
2673}
2674
1f6c58c1 2675void AliGenPythia::SetSeed(UInt_t seed)
2676{
2677 GetRandom()->SetSeed(seed);
2678}