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