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