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