Include full Z0/gamma* interference in Pythia
[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;
df607629 527 case kPyWPWHG:
589380c6 528 case kPyW:
0f6ee828 529 case kPyZ:
df607629 530 case kPyZgamma:
9a8774a1 531 case kPyMBRSingleDiffraction:
532 case kPyMBRDoubleDiffraction:
533 case kPyMBRCentralDiffraction:
589380c6 534 break;
8d2cd130 535 }
536//
592f8307 537//
538// JetFinder for Trigger
539//
540// Configure detector (EMCAL like)
541//
d7de4a67 542 fPythia->SetPARU(51, fPycellEtaMax);
543 fPythia->SetMSTU(51, fPycellNEta);
544 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 545//
546// Configure Jet Finder
547//
d7de4a67 548 fPythia->SetPARU(58, fPycellThreshold);
549 fPythia->SetPARU(52, fPycellEtSeed);
550 fPythia->SetPARU(53, fPycellMinEtJet);
551 fPythia->SetPARU(54, fPycellMaxRadius);
552 fPythia->SetMSTU(54, 2);
592f8307 553//
8d2cd130 554// This counts the total number of calls to Pyevnt() per run.
555 fTrialsRun = 0;
556 fQ = 0.;
557 fX1 = 0.;
558 fX2 = 0.;
559 fNev = 0 ;
560//
1d568bc2 561//
562//
8d2cd130 563 AliGenMC::Init();
fb355e51 564
565 // Reset Lorentz boost if demanded
566 if(!fUseLorentzBoost) {
567 fDyBoost = 0;
568 Warning("Init","Demand to discard Lorentz boost.\n");
569 }
1d568bc2 570//
571//
572//
573 if (fSetNuclei) {
fb355e51 574 fDyBoost = 0;
575 Warning("Init","Deprecated function SetNuclei() used (nPDFs + no boost). Use SetProjectile + SetTarget + SetUseNuclearPDF + SetUseLorentzBoost instead.\n");
1d568bc2 576 }
cd07c39b 577 fPythia->SetPARJ(200, 0.0);
578 fPythia->SetPARJ(199, 0.0);
579 fPythia->SetPARJ(198, 0.0);
580 fPythia->SetPARJ(197, 0.0);
581
582 if (fQuench == 1) {
5fa4b20b 583 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 584 }
3a709cfa 585
66b8652c 586 if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
587
b976f7d8 588 if (fQuench == 3) {
589 // Nestor's change of the splittings
590 fPythia->SetPARJ(200, 0.8);
591 fPythia->SetMSTJ(41, 1); // QCD radiation only
592 fPythia->SetMSTJ(42, 2); // angular ordering
593 fPythia->SetMSTJ(44, 2); // option to run alpha_s
594 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
595 fPythia->SetMSTJ(50, 0); // No coherence in first branching
596 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
cd07c39b 597 } else if (fQuench == 4) {
598 // Armesto-Cunqueiro-Salgado change of the splittings.
599 AliFastGlauber* glauber = AliFastGlauber::Instance();
600 glauber->Init(2);
601 //read and store transverse almonds corresponding to differnt
602 //impact parameters.
cd07c39b 603 glauber->SetCentralityClass(0.,0.1);
604 fPythia->SetPARJ(200, 1.);
605 fPythia->SetPARJ(198, fQhat);
606 fPythia->SetPARJ(199, fLength);
cd07c39b 607 fPythia->SetMSTJ(42, 2); // angular ordering
608 fPythia->SetMSTJ(44, 2); // option to run alpha_s
cd07c39b 609 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
b976f7d8 610 }
df607629 611
612 StdoutToAliDebug(1, fPythia->Pystat(4));
613 StdoutToAliDebug(1, fPythia->Pystat(5));
8d2cd130 614}
615
616void AliGenPythia::Generate()
617{
618// Generate one event
19ef8cf0 619 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
8d2cd130 620 fDecayer->ForceDecay();
621
13cca24a 622 Double_t polar[3] = {0,0,0};
623 Double_t origin[3] = {0,0,0};
624 Double_t p[4];
8d2cd130 625// converts from mm/c to s
626 const Float_t kconv=0.001/2.999792458e8;
627//
628 Int_t nt=0;
629 Int_t jev=0;
630 Int_t j, kf;
631 fTrials=0;
f913ec4f 632 fEventTime = 0.;
633
2590ccf9 634
8d2cd130 635
636 // Set collision vertex position
2590ccf9 637 if (fVertexSmear == kPerEvent) Vertex();
638
8d2cd130 639// event loop
640 while(1)
641 {
32d6ef7d 642//
5fa4b20b 643// Produce event
32d6ef7d 644//
32d6ef7d 645//
5fa4b20b 646// Switch hadronisation off
647//
648 fPythia->SetMSTJ(1, 0);
cd07c39b 649
650 if (fQuench ==4){
651 Double_t bimp;
652 // Quenching comes through medium-modified splitting functions.
653 AliFastGlauber::Instance()->GetRandomBHard(bimp);
e6fe9b82 654 fPythia->SetPARJ(197, bimp);
655 fImpact = bimp;
6c43eccb 656 fPythia->Qpygin0();
cd07c39b 657 }
32d6ef7d 658//
5fa4b20b 659// Either produce new event or read partons from file
660//
661 if (!fReadFromFile) {
beac474c 662 if (!fNewMIS) {
663 fPythia->Pyevnt();
664 } else {
665 fPythia->Pyevnw();
666 }
5fa4b20b 667 fNpartons = fPythia->GetN();
668 } else {
33c3c91a 669 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
670 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
5fa4b20b 671 fPythia->SetN(0);
672 LoadEvent(fRL->Stack(), 0 , 1);
673 fPythia->Pyedit(21);
674 }
675
32d6ef7d 676//
677// Run quenching routine
678//
5fa4b20b 679 if (fQuench == 1) {
680 fPythia->Quench();
681 } else if (fQuench == 2){
682 fPythia->Pyquen(208., 0, 0.);
b976f7d8 683 } else if (fQuench == 3) {
684 // Quenching is via multiplicative correction of the splittings
5fa4b20b 685 }
b976f7d8 686
32d6ef7d 687//
5fa4b20b 688// Switch hadronisation on
32d6ef7d 689//
20e47f08 690 if (fHadronisation) {
691 fPythia->SetMSTJ(1, 1);
5fa4b20b 692//
693// .. and perform hadronisation
aea21c57 694// printf("Calling hadronisation %d\n", fPythia->GetN());
03358a32 695
696 if (fPatchOmegaDalitz) {
697 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
698 fPythia->Pyexec();
699 fPythia->DalitzDecays();
700 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
701 }
702 fPythia->Pyexec();
20e47f08 703 }
8d2cd130 704 fTrials++;
8507138f 705 fPythia->ImportParticles(&fParticles,"All");
03358a32 706
df275cfa 707 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
e1cf8975 708 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
f64640b5 709
8d2cd130 710//
711//
712//
713 Int_t i;
714
dbd64db6 715 fNprimaries = 0;
8507138f 716 Int_t np = fParticles.GetEntriesFast();
5fa4b20b 717
7c21f297 718 if (np == 0) continue;
8d2cd130 719//
2590ccf9 720
8d2cd130 721//
722 Int_t* pParent = new Int_t[np];
723 Int_t* pSelected = new Int_t[np];
724 Int_t* trackIt = new Int_t[np];
5fa4b20b 725 for (i = 0; i < np; i++) {
8d2cd130 726 pParent[i] = -1;
727 pSelected[i] = 0;
728 trackIt[i] = 0;
729 }
730
731 Int_t nc = 0; // Total n. of selected particles
732 Int_t nParents = 0; // Selected parents
733 Int_t nTkbles = 0; // Trackable particles
f529e69b 734 if (fProcess != kPyMbDefault &&
735 fProcess != kPyMb &&
6d2ec66d 736 fProcess != kPyMbAtlasTuneMC09 &&
04081a91 737 fProcess != kPyMbWithDirectPhoton &&
f529e69b 738 fProcess != kPyJets &&
8d2cd130 739 fProcess != kPyDirectGamma &&
589380c6 740 fProcess != kPyMbNonDiffr &&
d7de4a67 741 fProcess != kPyMbMSEL1 &&
f529e69b 742 fProcess != kPyW &&
743 fProcess != kPyZ &&
df607629 744 fProcess != kPyZgamma &&
f529e69b 745 fProcess != kPyCharmppMNRwmi &&
9dfe63b3 746 fProcess != kPyBeautyppMNRwmi &&
64da86aa 747 fProcess != kPyBeautyJets &&
df607629 748 fProcess != kPyWPWHG &&
64da86aa 749 fProcess != kPyJetsPWHG &&
750 fProcess != kPyCharmPWHG &&
751 fProcess != kPyBeautyPWHG) {
8d2cd130 752
5fa4b20b 753 for (i = 0; i < np; i++) {
8507138f 754 TParticle* iparticle = (TParticle *) fParticles.At(i);
8d2cd130 755 Int_t ks = iparticle->GetStatusCode();
756 kf = CheckPDGCode(iparticle->GetPdgCode());
757// No initial state partons
758 if (ks==21) continue;
759//
760// Heavy Flavor Selection
761//
762 // quark ?
763 kf = TMath::Abs(kf);
764 Int_t kfl = kf;
9ff6c04c 765 // Resonance
f913ec4f 766
9ff6c04c 767 if (kfl > 100000) kfl %= 100000;
183a5ca9 768 if (kfl > 10000) kfl %= 10000;
8d2cd130 769 // meson ?
770 if (kfl > 10) kfl/=100;
771 // baryon
772 if (kfl > 10) kfl/=10;
8d2cd130 773 Int_t ipa = iparticle->GetFirstMother()-1;
774 Int_t kfMo = 0;
f913ec4f 775//
776// Establish mother daughter relation between heavy quarks and mesons
777//
778 if (kf >= fFlavorSelect && kf <= 6) {
779 Int_t idau = iparticle->GetFirstDaughter() - 1;
780 if (idau > -1) {
8507138f 781 TParticle* daughter = (TParticle *) fParticles.At(idau);
f913ec4f 782 Int_t pdgD = daughter->GetPdgCode();
783 if (pdgD == 91 || pdgD == 92) {
784 Int_t jmin = daughter->GetFirstDaughter() - 1;
785 Int_t jmax = daughter->GetLastDaughter() - 1;
2ab330c9 786 for (Int_t jp = jmin; jp <= jmax; jp++)
787 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
f913ec4f 788 } // is string or cluster
789 } // has daughter
790 } // heavy quark
8d2cd130 791
f913ec4f 792
8d2cd130 793 if (ipa > -1) {
8507138f 794 TParticle * mother = (TParticle *) fParticles.At(ipa);
8d2cd130 795 kfMo = TMath::Abs(mother->GetPdgCode());
796 }
f913ec4f 797
8d2cd130 798 // What to keep in Stack?
799 Bool_t flavorOK = kFALSE;
800 Bool_t selectOK = kFALSE;
801 if (fFeedDownOpt) {
32d6ef7d 802 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 803 } else {
32d6ef7d 804 if (kfl > fFlavorSelect) {
805 nc = -1;
806 break;
807 }
808 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 809 }
810 switch (fStackFillOpt) {
06956108 811 case kHeavyFlavor:
8d2cd130 812 case kFlavorSelection:
32d6ef7d 813 selectOK = kTRUE;
814 break;
8d2cd130 815 case kParentSelection:
32d6ef7d 816 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
817 break;
8d2cd130 818 }
819 if (flavorOK && selectOK) {
820//
821// Heavy flavor hadron or quark
822//
823// Kinematic seletion on final state heavy flavor mesons
824 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
825 {
9ff6c04c 826 continue;
8d2cd130 827 }
828 pSelected[i] = 1;
829 if (ParentSelected(kf)) ++nParents; // Update parent count
830// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
831 } else {
832// Kinematic seletion on decay products
833 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 834 && !KinematicSelection(iparticle, 1))
8d2cd130 835 {
9ff6c04c 836 continue;
8d2cd130 837 }
838//
839// Decay products
840// Select if mother was selected and is not tracked
841
842 if (pSelected[ipa] &&
843 !trackIt[ipa] && // mother will be tracked ?
844 kfMo != 5 && // mother is b-quark, don't store fragments
845 kfMo != 4 && // mother is c-quark, don't store fragments
846 kf != 92) // don't store string
847 {
848//
849// Semi-stable or de-selected: diselect decay products:
850//
851//
852 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
853 {
854 Int_t ipF = iparticle->GetFirstDaughter();
855 Int_t ipL = iparticle->GetLastDaughter();
856 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
857 }
858// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
859 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
860 }
861 }
862 if (pSelected[i] == -1) pSelected[i] = 0;
863 if (!pSelected[i]) continue;
864 // Count quarks only if you did not include fragmentation
865 if (fFragmentation && kf <= 10) continue;
9ff6c04c 866
8d2cd130 867 nc++;
868// Decision on tracking
869 trackIt[i] = 0;
870//
871// Track final state particle
872 if (ks == 1) trackIt[i] = 1;
873// Track semi-stable particles
d25cfd65 874 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 875// Track particles selected by process if undecayed.
876 if (fForceDecay == kNoDecay) {
877 if (ParentSelected(kf)) trackIt[i] = 1;
878 } else {
879 if (ParentSelected(kf)) trackIt[i] = 0;
880 }
881 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
882//
883//
884
885 } // particle selection loop
886 if (nc > 0) {
887 for (i = 0; i<np; i++) {
888 if (!pSelected[i]) continue;
8507138f 889 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 890 kf = CheckPDGCode(iparticle->GetPdgCode());
891 Int_t ks = iparticle->GetStatusCode();
892 p[0] = iparticle->Px();
893 p[1] = iparticle->Py();
894 p[2] = iparticle->Pz();
a920faf9 895 p[3] = iparticle->Energy();
896
2590ccf9 897 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
898 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
899 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
900
21391258 901 Float_t tof = fTime + kconv*iparticle->T();
8d2cd130 902 Int_t ipa = iparticle->GetFirstMother()-1;
903 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 904
905 PushTrack(fTrackIt*trackIt[i], iparent, kf,
906 p[0], p[1], p[2], p[3],
907 origin[0], origin[1], origin[2], tof,
908 polar[0], polar[1], polar[2],
909 kPPrimary, nt, 1., ks);
8d2cd130 910 pParent[i] = nt;
dbd64db6 911 KeepTrack(nt);
912 fNprimaries++;
642f15cf 913 } // PushTrack loop
8d2cd130 914 }
915 } else {
916 nc = GenerateMB();
917 } // mb ?
f913ec4f 918
919 GetSubEventTime();
8d2cd130 920
234f6d32 921 delete[] pParent;
922 delete[] pSelected;
923 delete[] trackIt;
8d2cd130 924
925 if (nc > 0) {
926 switch (fCountMode) {
927 case kCountAll:
928 // printf(" Count all \n");
929 jev += nc;
930 break;
931 case kCountParents:
932 // printf(" Count parents \n");
933 jev += nParents;
934 break;
935 case kCountTrackables:
936 // printf(" Count trackable \n");
937 jev += nTkbles;
938 break;
939 }
940 if (jev >= fNpart || fNpart == -1) {
941 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 942
8d2cd130 943 fQ += fPythia->GetVINT(51);
944 fX1 += fPythia->GetVINT(41);
945 fX2 += fPythia->GetVINT(42);
946 fTrialsRun += fTrials;
947 fNev++;
948 MakeHeader();
949 break;
950 }
951 }
952 } // event loop
953 SetHighWaterMark(nt);
954// adjust weight due to kinematic selection
b88f5cea 955// AdjustWeights();
8d2cd130 956// get cross-section
957 fXsection=fPythia->GetPARI(1);
958}
959
960Int_t AliGenPythia::GenerateMB()
961{
962//
963// Min Bias selection and other global selections
964//
965 Int_t i, kf, nt, iparent;
966 Int_t nc = 0;
13cca24a 967 Double_t p[4];
968 Double_t polar[3] = {0,0,0};
969 Double_t origin[3] = {0,0,0};
8d2cd130 970// converts from mm/c to s
971 const Float_t kconv=0.001/2.999792458e8;
972
e0e89f40 973
8507138f 974 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
aea21c57 975
5fa4b20b 976
e0e89f40 977
8d2cd130 978 Int_t* pParent = new Int_t[np];
979 for (i=0; i< np; i++) pParent[i] = -1;
64da86aa 980 if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
0e9e66f0 981 && fEtMinJet > 0.) {
bed3df2c 982 TParticle* jet1 = (TParticle *) fParticles.At(6);
8507138f 983 TParticle* jet2 = (TParticle *) fParticles.At(7);
bed3df2c 984
985 if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
234f6d32 986 delete [] pParent;
987 return 0;
988 }
8d2cd130 989 }
e0e89f40 990
40fe59d4 991
992 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
ac01c7c6 993 // implemented primaryly for kPyJets, but extended to any kind of process.
40fe59d4 994 if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
995 (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
996 Bool_t ok = TriggerOnSelectedParticles(np);
43e3b80a 997
998 if(!ok) {
999 delete[] pParent;
1000 return 0;
ec2c406e 1001 }
43e3b80a 1002 }
9dfe63b3 1003
800be8b7 1004 // Check for diffraction
1005 if(fkTuneForDiff) {
c5dfa3e4 1006 if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
800be8b7 1007 if(!CheckDiffraction()) {
1008 delete [] pParent;
1009 return 0;
1010 }
1011 }
1012 }
1013
700b9416 1014 // Check for minimum multiplicity
1015 if (fTriggerMultiplicity > 0) {
1016 Int_t multiplicity = 0;
1017 for (i = 0; i < np; i++) {
8507138f 1018 TParticle * iparticle = (TParticle *) fParticles.At(i);
700b9416 1019
1020 Int_t statusCode = iparticle->GetStatusCode();
1021
1022 // Initial state particle
e296848e 1023 if (statusCode != 1)
700b9416 1024 continue;
38112f3f 1025 // eta cut
700b9416 1026 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1027 continue;
38112f3f 1028 // pt cut
1029 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1030 continue;
1031
700b9416 1032 TParticlePDG* pdgPart = iparticle->GetPDG();
1033 if (pdgPart && pdgPart->Charge() == 0)
1034 continue;
1035
1036 ++multiplicity;
1037 }
1038
1039 if (multiplicity < fTriggerMultiplicity) {
1040 delete [] pParent;
1041 return 0;
1042 }
38112f3f 1043 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
700b9416 1044 }
90a236ce 1045
90a236ce 1046
7ce1321b 1047 if (fTriggerParticle) {
1048 Bool_t triggered = kFALSE;
1049 for (i = 0; i < np; i++) {
8507138f 1050 TParticle * iparticle = (TParticle *) fParticles.At(i);
7ce1321b 1051 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 1052 if (kf != fTriggerParticle) continue;
7ce1321b 1053 if (iparticle->Pt() == 0.) continue;
1054 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
2591bd0e 1055 if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
7ce1321b 1056 triggered = kTRUE;
1057 break;
1058 }
234f6d32 1059 if (!triggered) {
1060 delete [] pParent;
1061 return 0;
1062 }
7ce1321b 1063 }
e0e89f40 1064
1065
1066 // Check if there is a ccbar or bbbar pair with at least one of the two
1067 // in fYMin < y < fYMax
2f405d65 1068
9dfe63b3 1069 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
9ccc3d0e 1070 TParticle *partCheck;
1071 TParticle *mother;
e0e89f40 1072 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 1073 Bool_t theChild=kFALSE;
5d76923e 1074 Bool_t triggered=kFALSE;
9ccc3d0e 1075 Float_t y;
1076 Int_t pdg,mpdg,mpdgUpperFamily;
e0e89f40 1077 for(i=0; i<np; i++) {
9ccc3d0e 1078 partCheck = (TParticle*)fParticles.At(i);
1079 pdg = partCheck->GetPdgCode();
1080 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1081 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1082 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1083 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1084 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1085 }
5d76923e 1086 if(fTriggerParticle) {
1087 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1088 }
9ccc3d0e 1089 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1090 Int_t mi = partCheck->GetFirstMother() - 1;
1091 if(mi<0) continue;
1092 mother = (TParticle*)fParticles.At(mi);
1093 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 1094 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 1095 if ( ParentSelected(mpdg) ||
1096 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1097 if (KinematicSelection(partCheck,1)) {
1098 theChild=kTRUE;
1099 }
1100 }
1101 }
e0e89f40 1102 }
9ccc3d0e 1103 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
234f6d32 1104 delete[] pParent;
e0e89f40 1105 return 0;
1106 }
9ccc3d0e 1107 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1108 delete[] pParent;
1109 return 0;
1110 }
5d76923e 1111 if(fTriggerParticle && !triggered) { // particle requested is not produced
1112 delete[] pParent;
1113 return 0;
1114 }
9ccc3d0e 1115
e0e89f40 1116 }
aea21c57 1117
0f6ee828 1118 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
e136a92a 1119 if ( (
df607629 1120 fProcess == kPyWPWHG ||
e136a92a 1121 fProcess == kPyW ||
f529e69b 1122 fProcess == kPyZ ||
df607629 1123 fProcess == kPyZgamma ||
f529e69b 1124 fProcess == kPyMbDefault ||
1125 fProcess == kPyMb ||
6d2ec66d 1126 fProcess == kPyMbAtlasTuneMC09 ||
04081a91 1127 fProcess == kPyMbWithDirectPhoton ||
f529e69b 1128 fProcess == kPyMbNonDiffr)
0f6ee828 1129 && (fCutOnChild == 1) ) {
1130 if ( !CheckKinematicsOnChild() ) {
234f6d32 1131 delete[] pParent;
0f6ee828 1132 return 0;
1133 }
aea21c57 1134 }
1135
1136
f913ec4f 1137 for (i = 0; i < np; i++) {
8d2cd130 1138 Int_t trackIt = 0;
8507138f 1139 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 1140 kf = CheckPDGCode(iparticle->GetPdgCode());
1141 Int_t ks = iparticle->GetStatusCode();
1142 Int_t km = iparticle->GetFirstMother();
06956108 1143 if (
1144 (((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
64da86aa 1145 ((fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) && ks == 21 && km == 0 && i>1)
06956108 1146 )
1147 {
1148 nc++;
8d2cd130 1149 if (ks == 1) trackIt = 1;
1150 Int_t ipa = iparticle->GetFirstMother()-1;
1151
1152 iparent = (ipa > -1) ? pParent[ipa] : -1;
1153
1154//
1155// store track information
1156 p[0] = iparticle->Px();
1157 p[1] = iparticle->Py();
1158 p[2] = iparticle->Pz();
a920faf9 1159 p[3] = iparticle->Energy();
1406f599 1160
a920faf9 1161
2590ccf9 1162 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1163 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1164 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1165
21391258 1166 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
a920faf9 1167
1168 PushTrack(fTrackIt*trackIt, iparent, kf,
1169 p[0], p[1], p[2], p[3],
1170 origin[0], origin[1], origin[2], tof,
1171 polar[0], polar[1], polar[2],
1172 kPPrimary, nt, 1., ks);
dbd64db6 1173 fNprimaries++;
8d2cd130 1174 KeepTrack(nt);
1175 pParent[i] = nt;
f913ec4f 1176 SetHighWaterMark(nt);
1177
8d2cd130 1178 } // select particle
1179 } // particle loop
1180
234f6d32 1181 delete[] pParent;
e0e89f40 1182
f913ec4f 1183 return 1;
8d2cd130 1184}
1185
1186
1187void AliGenPythia::FinishRun()
1188{
1189// Print x-section summary
1190 fPythia->Pystat(1);
2779fc64 1191
1192 if (fNev > 0.) {
1193 fQ /= fNev;
1194 fX1 /= fNev;
1195 fX2 /= fNev;
1196 }
1197
8d2cd130 1198 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1199 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1200}
1201
7184e472 1202void AliGenPythia::AdjustWeights() const
8d2cd130 1203{
1204// Adjust the weights after generation of all events
1205//
e2bddf81 1206 if (gAlice) {
1207 TParticle *part;
1208 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1209 for (Int_t i=0; i<ntrack; i++) {
1210 part= gAlice->GetMCApp()->Particle(i);
1211 part->SetWeight(part->GetWeight()*fKineBias);
1212 }
8d2cd130 1213 }
1214}
1215
20e47f08 1216void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
8d2cd130 1217{
1218// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1219
fb355e51 1220 fAProjectile = a1;
1221 fATarget = a2;
1222 fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
1223 fUseNuclearPDF = kTRUE;
1224 fSetNuclei = kTRUE;
8d2cd130 1225}
1226
1227
1228void AliGenPythia::MakeHeader()
1229{
7184e472 1230//
1231// Make header for the simulated event
1232//
183a5ca9 1233 if (gAlice) {
1234 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1235 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1236 }
1237
8d2cd130 1238// Builds the event header, to be called after each event
e5c87a3d 1239 if (fHeader) delete fHeader;
1240 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1241//
1242// Event type
800be8b7 1243 if(fProcDiff>0){
1244 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1245 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1246 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1247 }
1248 else
e5c87a3d 1249 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1250//
1251// Number of trials
e5c87a3d 1252 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1253//
1254// Event Vertex
d25cfd65 1255 fHeader->SetPrimaryVertex(fVertex);
21391258 1256 fHeader->SetInteractionTime(fTime+fEventTime);
dbd64db6 1257//
1258// Number of primaries
1259 fHeader->SetNProduced(fNprimaries);
8d2cd130 1260//
a0027228 1261// Event weight
1262 fHeader->SetEventWeight(fPythia->GetVINT(97));
1263//
8d2cd130 1264// Jets that have triggered
f913ec4f 1265
9dfe63b3 1266 //Need to store jets for b-jet studies too!
64da86aa 1267 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
8d2cd130 1268 {
1269 Int_t ntrig, njet;
1270 Float_t jets[4][10];
1271 GetJets(njet, ntrig, jets);
9ff6c04c 1272
8d2cd130 1273
1274 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1275 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1276 jets[3][i]);
1277 }
1278 }
5fa4b20b 1279//
1280// Copy relevant information from external header, if present.
1281//
1282 Float_t uqJet[4];
1283
1284 if (fRL) {
1285 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1286 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1287 {
1288 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1289
1290
1291 exHeader->TriggerJet(i, uqJet);
1292 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1293 }
1294 }
1295//
1296// Store quenching parameters
1297//
1298 if (fQuench){
b6262a45 1299 Double_t z[4] = {0.};
1300 Double_t xp = 0.;
1301 Double_t yp = 0.;
1302
7c21f297 1303 if (fQuench == 1) {
1304 // Pythia::Quench()
1305 fPythia->GetQuenchingParameters(xp, yp, z);
1bab4b79 1306 } else if (fQuench == 2){
7c21f297 1307 // Pyquen
1308 Double_t r1 = PARIMP.rb1;
1309 Double_t r2 = PARIMP.rb2;
1310 Double_t b = PARIMP.b1;
1311 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1312 Double_t phi = PARIMP.psib1;
1313 xp = r * TMath::Cos(phi);
1314 yp = r * TMath::Sin(phi);
1315
1bab4b79 1316 } else if (fQuench == 4) {
1317 // QPythia
5831e75f 1318 Double_t xy[2];
e9719084 1319 Double_t i0i1[2];
1320 AliFastGlauber::Instance()->GetSavedXY(xy);
1321 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
5831e75f 1322 xp = xy[0];
1323 yp = xy[1];
e6fe9b82 1324 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
7c21f297 1325 }
1bab4b79 1326
7c21f297 1327 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1328 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1bab4b79 1329 }
beac474c 1330//
39ea8b7f 1331// Store pt^hard and cross-section
beac474c 1332 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
39ea8b7f 1333 ((AliGenPythiaEventHeader*) fHeader)->SetXsection(fPythia->GetPARI(1));
64da86aa 1334
1335//
1336// Store Event Weight
39ea8b7f 1337 ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7));
64da86aa 1338
5fa4b20b 1339//
cf57b268 1340// Pass header
5fa4b20b 1341//
cf57b268 1342 AddHeader(fHeader);
4c4eac97 1343 fHeader = 0x0;
8d2cd130 1344}
cf57b268 1345
c2fc9d31 1346Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
8d2cd130 1347{
1348// Check the kinematic trigger condition
1349//
1350 Double_t eta[2];
1351 eta[0] = jet1->Eta();
1352 eta[1] = jet2->Eta();
1353 Double_t phi[2];
1354 phi[0] = jet1->Phi();
1355 phi[1] = jet2->Phi();
1356 Int_t pdg[2];
1357 pdg[0] = jet1->GetPdgCode();
1358 pdg[1] = jet2->GetPdgCode();
1359 Bool_t triggered = kFALSE;
1360
64da86aa 1361 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) {
8d2cd130 1362 Int_t njets = 0;
1363 Int_t ntrig = 0;
1364 Float_t jets[4][10];
1365//
1366// Use Pythia clustering on parton level to determine jet axis
1367//
1368 GetJets(njets, ntrig, jets);
1369
76d6ba9a 1370 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1371//
1372 } else {
1373 Int_t ij = 0;
1374 Int_t ig = 1;
1375 if (pdg[0] == kGamma) {
1376 ij = 1;
1377 ig = 0;
1378 }
1379 //Check eta range first...
1380 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1381 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1382 {
1383 //Eta is okay, now check phi range
1384 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1385 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1386 {
1387 triggered = kTRUE;
1388 }
1389 }
1390 }
1391 return triggered;
1392}
aea21c57 1393
1394
aea21c57 1395
7184e472 1396Bool_t AliGenPythia::CheckKinematicsOnChild(){
1397//
1398//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1399//
aea21c57 1400 Bool_t checking = kFALSE;
1401 Int_t j, kcode, ks, km;
1402 Int_t nPartAcc = 0; //number of particles in the acceptance range
1403 Int_t numberOfAcceptedParticles = 1;
1404 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1405 Int_t npart = fParticles.GetEntriesFast();
aea21c57 1406
0f6ee828 1407 for (j = 0; j<npart; j++) {
8507138f 1408 TParticle * jparticle = (TParticle *) fParticles.At(j);
aea21c57 1409 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1410 ks = jparticle->GetStatusCode();
1411 km = jparticle->GetFirstMother();
1412
1413 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1414 nPartAcc++;
1415 }
0f6ee828 1416 if( numberOfAcceptedParticles <= nPartAcc){
1417 checking = kTRUE;
1418 break;
1419 }
aea21c57 1420 }
0f6ee828 1421
aea21c57 1422 return checking;
aea21c57 1423}
1424
5fa4b20b 1425void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1426{
1058d9df 1427 //
1428 // Load event into Pythia Common Block
1429 //
1430
1431 Int_t npart = stack -> GetNprimary();
1432 Int_t n0 = 0;
1433
1434 if (!flag) {
1435 (fPythia->GetPyjets())->N = npart;
1436 } else {
1437 n0 = (fPythia->GetPyjets())->N;
1438 (fPythia->GetPyjets())->N = n0 + npart;
1439 }
1440
1441
1442 for (Int_t part = 0; part < npart; part++) {
1443 TParticle *mPart = stack->Particle(part);
32d6ef7d 1444
1058d9df 1445 Int_t kf = mPart->GetPdgCode();
1446 Int_t ks = mPart->GetStatusCode();
1447 Int_t idf = mPart->GetFirstDaughter();
1448 Int_t idl = mPart->GetLastDaughter();
1449
1450 if (reHadr) {
1451 if (ks == 11 || ks == 12) {
1452 ks -= 10;
1453 idf = -1;
1454 idl = -1;
1455 }
32d6ef7d 1456 }
1457
1058d9df 1458 Float_t px = mPart->Px();
1459 Float_t py = mPart->Py();
1460 Float_t pz = mPart->Pz();
1461 Float_t e = mPart->Energy();
1462 Float_t m = mPart->GetCalcMass();
32d6ef7d 1463
1058d9df 1464
1465 (fPythia->GetPyjets())->P[0][part+n0] = px;
1466 (fPythia->GetPyjets())->P[1][part+n0] = py;
1467 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1468 (fPythia->GetPyjets())->P[3][part+n0] = e;
1469 (fPythia->GetPyjets())->P[4][part+n0] = m;
1470
1471 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1472 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1473 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1474 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1475 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1476 }
1477}
1478
c2fc9d31 1479void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1058d9df 1480{
1481 //
1482 // Load event into Pythia Common Block
1483 //
1484
1485 Int_t npart = stack -> GetEntries();
1486 Int_t n0 = 0;
1487
1488 if (!flag) {
1489 (fPythia->GetPyjets())->N = npart;
1490 } else {
1491 n0 = (fPythia->GetPyjets())->N;
1492 (fPythia->GetPyjets())->N = n0 + npart;
1493 }
1494
1495
1496 for (Int_t part = 0; part < npart; part++) {
1497 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
92847124 1498 if (!mPart) continue;
1499
1058d9df 1500 Int_t kf = mPart->GetPdgCode();
1501 Int_t ks = mPart->GetStatusCode();
1502 Int_t idf = mPart->GetFirstDaughter();
1503 Int_t idl = mPart->GetLastDaughter();
1504
1505 if (reHadr) {
92847124 1506 if (ks == 11 || ks == 12) {
1507 ks -= 10;
1508 idf = -1;
1509 idl = -1;
1510 }
8d2cd130 1511 }
1058d9df 1512
1513 Float_t px = mPart->Px();
1514 Float_t py = mPart->Py();
1515 Float_t pz = mPart->Pz();
1516 Float_t e = mPart->Energy();
1517 Float_t m = mPart->GetCalcMass();
1518
1519
1520 (fPythia->GetPyjets())->P[0][part+n0] = px;
1521 (fPythia->GetPyjets())->P[1][part+n0] = py;
1522 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1523 (fPythia->GetPyjets())->P[3][part+n0] = e;
1524 (fPythia->GetPyjets())->P[4][part+n0] = m;
1525
1526 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1527 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1528 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1529 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1530 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1531 }
8d2cd130 1532}
1533
5fa4b20b 1534
014a9521 1535void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1536{
1537//
1538// Calls the Pythia jet finding algorithm to find jets in the current event
1539//
1540//
8d2cd130 1541//
1542// Save jets
1543 Int_t n = fPythia->GetN();
1544
1545//
1546// Run Jet Finder
1547 fPythia->Pycell(njets);
1548 Int_t i;
1549 for (i = 0; i < njets; i++) {
1550 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1551 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1552 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1553 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1554
1555 jets[0][i] = px;
1556 jets[1][i] = py;
1557 jets[2][i] = pz;
1558 jets[3][i] = e;
1559 }
1560}
1561
1562
1563
1564void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1565{
1566//
1567// Calls the Pythia clustering algorithm to find jets in the current event
1568//
1569 Int_t n = fPythia->GetN();
1570 nJets = 0;
1571 nJetsTrig = 0;
1572 if (fJetReconstruction == kCluster) {
1573//
1574// Configure cluster algorithm
1575//
1576 fPythia->SetPARU(43, 2.);
1577 fPythia->SetMSTU(41, 1);
1578//
1579// Call cluster algorithm
1580//
1581 fPythia->Pyclus(nJets);
1582//
1583// Loading jets from common block
1584//
1585 } else {
592f8307 1586
8d2cd130 1587//
1588// Run Jet Finder
1589 fPythia->Pycell(nJets);
1590 }
1591
1592 Int_t i;
1593 for (i = 0; i < nJets; i++) {
1594 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1595 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1596 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1597 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1598 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1599 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1600 Float_t theta = TMath::ATan2(pt,pz);
1601 Float_t et = e * TMath::Sin(theta);
1602 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1603 if (
1604 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1605 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1606 et > fEtMinJet && et < fEtMaxJet
1607 )
1608 {
1609 jets[0][nJetsTrig] = px;
1610 jets[1][nJetsTrig] = py;
1611 jets[2][nJetsTrig] = pz;
1612 jets[3][nJetsTrig] = e;
1613 nJetsTrig++;
5fa4b20b 1614// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1615 } else {
1616// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1617 }
1618 }
1619}
1620
f913ec4f 1621void AliGenPythia::GetSubEventTime()
1622{
1623 // Calculates time of the next subevent
9d7108a7 1624 fEventTime = 0.;
1625 if (fEventsTime) {
1626 TArrayF &array = *fEventsTime;
1627 fEventTime = array[fCurSubEvent++];
1628 }
1629 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1630 return;
f913ec4f 1631}
8d2cd130 1632
e81f2bac 1633Bool_t AliGenPythia::IsInBarrel(Float_t eta) const
40fe59d4 1634{
1635 // Is particle in Central Barrel acceptance?
1636 // etamin=-etamax
1637 if( eta < fTriggerEta )
1638 return kTRUE;
1639 else
1640 return kFALSE;
1641}
1642
e81f2bac 1643Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
ec2c406e 1644{
1645 // Is particle in EMCAL acceptance?
ec2c406e 1646 // phi in degrees, etamin=-etamax
9fd8e520 1647 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1648 eta < fEMCALEta )
ec2c406e 1649 return kTRUE;
1650 else
1651 return kFALSE;
1652}
1653
e81f2bac 1654Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
ec2c406e 1655{
1656 // Is particle in PHOS acceptance?
1657 // Acceptance slightly larger considered.
1658 // phi in degrees, etamin=-etamax
40fe59d4 1659 // iparticle is the index of the particle to be checked, for PHOS rotation case
1660
9fd8e520 1661 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1662 eta < fPHOSEta )
ec2c406e 1663 return kTRUE;
1664 else
40fe59d4 1665 {
1666 if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1667
ec2c406e 1668 return kFALSE;
40fe59d4 1669 }
ec2c406e 1670}
1671
40fe59d4 1672void AliGenPythia::RotatePhi(Bool_t& okdd)
90a236ce 1673{
40fe59d4 1674 //Rotate event in phi to enhance events in PHOS acceptance
1675
1676 if(fPHOSRotateCandidate < 0) return ;
1677
90a236ce 1678 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1679 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1680 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
800be8b7 1681 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
90a236ce 1682
1683 //calculate deltaphi
40fe59d4 1684 TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
90a236ce 1685 Double_t phphi = ph->Phi();
1686 Double_t deltaphi = phiPHOS - phphi;
40fe59d4 1687
90a236ce 1688
1689
1690 //loop for all particles and produce the phi rotation
8507138f 1691 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
90a236ce 1692 Double_t oldphi, newphi;
777e69b0 1693 Double_t newVx, newVy, r, vZ, time;
1694 Double_t newPx, newPy, pt, pz, e;
90a236ce 1695 for(Int_t i=0; i< np; i++) {
40fe59d4 1696 TParticle* iparticle = (TParticle *) fParticles.At(i);
1697 oldphi = iparticle->Phi();
1698 newphi = oldphi + deltaphi;
1699 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1700 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1701
1702 r = iparticle->R();
1703 newVx = r * TMath::Cos(newphi);
1704 newVy = r * TMath::Sin(newphi);
1705 vZ = iparticle->Vz(); // don't transform
1706 time = iparticle->T(); // don't transform
1707
1708 pt = iparticle->Pt();
1709 newPx = pt * TMath::Cos(newphi);
1710 newPy = pt * TMath::Sin(newphi);
1711 pz = iparticle->Pz(); // don't transform
1712 e = iparticle->Energy(); // don't transform
1713
1714 // apply rotation
1715 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1716 iparticle->SetMomentum(newPx, newPy, pz, e);
1717
90a236ce 1718 } //end particle loop
1719
40fe59d4 1720 // now let's check that we put correctly the candidate photon in PHOS
1721 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1722 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1723 if(IsInPHOS(phi,eta,-1))
1724 okdd = kTRUE;
1725
1726 // reset the value for next event
1727 fPHOSRotateCandidate = -1;
1728
90a236ce 1729}
ec2c406e 1730
1731
800be8b7 1732Bool_t AliGenPythia::CheckDiffraction()
1733{
1734 // use this method only with Perugia-0 tune!
1735
1736 // printf("AAA\n");
1737
1738 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1739
1740 Int_t iPart1=-1;
1741 Int_t iPart2=-1;
1742
1743 Double_t y1 = 1e10;
1744 Double_t y2 = -1e10;
1745
1746 const Int_t kNstable=20;
1747 const Int_t pdgStable[20] = {
1748 22, // Photon
1749 11, // Electron
1750 12, // Electron Neutrino
1751 13, // Muon
1752 14, // Muon Neutrino
1753 15, // Tau
1754 16, // Tau Neutrino
1755 211, // Pion
1756 321, // Kaon
1757 311, // K0
1758 130, // K0s
1759 310, // K0l
1760 2212, // Proton
1761 2112, // Neutron
1762 3122, // Lambda_0
1763 3112, // Sigma Minus
1764 3222, // Sigma Plus
1765 3312, // Xsi Minus
1766 3322, // Xsi0
1767 3334 // Omega
1768 };
1769
1770 for (Int_t i = 0; i < np; i++) {
1771 TParticle * part = (TParticle *) fParticles.At(i);
1772
1773 Int_t statusCode = part->GetStatusCode();
1774
1775 // Initial state particle
1776 if (statusCode != 1)
1777 continue;
1778
1779 Int_t pdg = TMath::Abs(part->GetPdgCode());
1780 Bool_t isStable = kFALSE;
1781 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1782 if (pdg == pdgStable[i1]) {
1783 isStable = kTRUE;
1784 break;
1785 }
1786 }
1787 if(!isStable)
1788 continue;
1789
1790 Double_t y = part->Y();
1791
1792 if (y < y1)
1793 {
1794 y1 = y;
1795 iPart1 = i;
1796 }
1797 if (y > y2)
1798 {
1799 y2 = y;
1800 iPart2 = i;
1801 }
1802 }
1803
1804 if(iPart1<0 || iPart2<0) return kFALSE;
1805
1806 y1=TMath::Abs(y1);
1807 y2=TMath::Abs(y2);
1808
1809 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1810 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1811
1812 Int_t pdg1 = part1->GetPdgCode();
1813 Int_t pdg2 = part2->GetPdgCode();
1814
1815
1816 Int_t iPart = -1;
1817 if (pdg1 == 2212 && pdg2 == 2212)
1818 {
1819 if(y1 > y2)
1820 iPart = iPart1;
1821 else if(y1 < y2)
1822 iPart = iPart2;
1823 else {
1824 iPart = iPart1;
1825 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1826 }
1827 }
1828 else if (pdg1 == 2212)
1829 iPart = iPart1;
1830 else if (pdg2 == 2212)
1831 iPart = iPart2;
1832
1833
1834
1835
1836
1837 Double_t M=-1.;
1838 if(iPart>0) {
1839 TParticle * part = (TParticle *) fParticles.At(iPart);
1840 Double_t E= part->Energy();
1841 Double_t P= part->P();
4f813e90 1842 Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
1843 if(M2<0) return kFALSE;
1844 M= TMath::Sqrt(M2);
800be8b7 1845 }
1846
c5dfa3e4 1847 Double_t Mmin, Mmax, wSD, wDD, wND;
1848 if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1849
1850 if(M>-1 && M<Mmin) return kFALSE;
1851 if(M>Mmax) M=-1;
1852
1853 Int_t procType=fPythia->GetMSTI(1);
1854 Int_t proc0=2;
1855 if(procType== 94) proc0=1;
1856 if(procType== 92 || procType== 93) proc0=0;
1857
1858 Int_t proc=2;
1859 if(M>0) proc=0;
1860 else if(proc0==1) proc=1;
1861
1862 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1863 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
800be8b7 1864
1865
1866 // if(proc==1 || proc==2) return kFALSE;
1867
c5dfa3e4 1868 if(proc!=0) {
1869 if(proc0!=0) fProcDiff = procType;
1870 else fProcDiff = 95;
1871 return kTRUE;
1872 }
1873
1874 if(wSD<0) AliError("wSD<0 ! \n");
1875
1876 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1877
1878 // printf("iPart = %d\n", iPart);
1879
1880 if(iPart==iPart1) fProcDiff=93;
1881 else if(iPart==iPart2) fProcDiff=92;
1882 else {
1883 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1884
800be8b7 1885 }
1886
c5dfa3e4 1887 return kTRUE;
1888}
1889
1890
1891
1892Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1893 Double_t &wSD, Double_t &wDD, Double_t &wND)
1894{
1895
1896 // 900 GeV
1897 if(TMath::Abs(fEnergyCMS-900)<1 ){
1898
1899const Int_t nbin=400;
1900Double_t bin[]={
19011.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
19024.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
19037.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
190410.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
190513.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
190615.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
190718.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
190821.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
190924.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
191027.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
191130.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
191233.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
191336.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
191439.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
191542.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
191645.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
191748.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
191851.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
191954.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
192057.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
192160.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
192263.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
192366.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
192469.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
192572.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
192675.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
192778.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
192881.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
192984.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
193087.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
193190.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
193293.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
193396.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
193499.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1935102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1936105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1937108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1938111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1939114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1940117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1941120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1942123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1943126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
1944129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
1945132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
1946135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
1947138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
1948141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
1949144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
1950147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
1951150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
1952153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
1953156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
1954159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
1955162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
1956165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
1957168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
1958171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
1959174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
1960177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
1961180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
1962183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
1963186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
1964189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
1965192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
1966195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
1967198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1968Double_t w[]={
19691.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
19700.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
19710.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
19720.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
19730.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
19740.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
19750.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
19760.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
19770.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
19780.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
19790.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
19800.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
19810.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
19820.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
19830.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
19840.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
19850.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
19860.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
19870.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
19880.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
19890.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
19900.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
19910.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
19920.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
19930.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
19940.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
19950.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
19960.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
19970.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
19980.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
19990.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
20000.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
20010.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
20020.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
20030.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
20040.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
20050.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
20060.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
20070.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
20080.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
20090.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
20100.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
20110.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
20120.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
20130.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
20140.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
20150.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
20160.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
20170.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
20180.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
20190.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
20200.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
20210.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
20220.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
20230.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
20240.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
20250.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
20260.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
20270.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
20280.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
20290.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
20300.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
20310.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
20320.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
20330.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
20340.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
20350.386112, 0.364314, 0.434375, 0.334629};
2036wDD = 0.379611;
2037wND = 0.496961;
2038wSD = -1;
2039
2040 Mmin = bin[0];
2041 Mmax = bin[nbin];
2042 if(M<Mmin || M>Mmax) return kTRUE;
2043
7873275c 2044 Int_t ibin=nbin-1;
93a60586 2045 for(Int_t i=1; i<=nbin; i++)
2ff57420 2046 if(M<=bin[i]) {
2047 ibin=i-1;
2048 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
800be8b7 2049 break;
2050 }
c5dfa3e4 2051 wSD=w[ibin];
2052 return kTRUE;
2053 }
2054 else if(TMath::Abs(fEnergyCMS-2760)<1 ){
800be8b7 2055
c5dfa3e4 2056const Int_t nbin=400;
2057Double_t bin[]={
20581.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
20594.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
20607.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
206110.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
206213.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
206315.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
206418.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
206521.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
206624.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
206727.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
206830.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
206933.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
207036.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
207139.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
207242.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
207345.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
207448.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
207551.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
207654.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
207757.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
207860.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
207963.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
208066.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
208169.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
208272.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
208375.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
208478.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
208581.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
208684.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
208787.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
208890.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
208993.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
209096.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
209199.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2092102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2093105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2094108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2095111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2096114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2097117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2098120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2099123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2100126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2101129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2102132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2103135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2104138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2105141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2106144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2107147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2108150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2109153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2110156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2111159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2112162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2113165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2114168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2115171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2116174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2117177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2118180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2119183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2120186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2121189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2122192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2123195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2124198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2125Double_t w[]={
21261.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
21270.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
21280.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
21290.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
21300.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
21310.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
21320.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
21330.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
21340.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
21350.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
21360.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
21370.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
21380.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
21390.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
21400.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
21410.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
21420.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
21430.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
21440.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
21450.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
21460.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
21470.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
21480.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
21490.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
21500.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
21510.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
21520.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
21530.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
21540.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
21550.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
21560.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
21570.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
21580.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
21590.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
21600.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
21610.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
21620.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
21630.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
21640.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
21650.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
21660.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
21670.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
21680.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
21690.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
21700.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
21710.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
21720.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
21730.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
21740.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
21750.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
21760.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
21770.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
21780.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
21790.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
21800.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
21810.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
21820.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
21830.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
21840.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
21850.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
21860.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
21870.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
21880.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
21890.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
21900.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
21910.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
21920.201712, 0.242225, 0.255565, 0.258738};
2193wDD = 0.512813;
2194wND = 0.518820;
2195wSD = -1;
2196
2197 Mmin = bin[0];
2198 Mmax = bin[nbin];
2199 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2200
c5dfa3e4 2201 Int_t ibin=nbin-1;
2202 for(Int_t i=1; i<=nbin; i++)
2203 if(M<=bin[i]) {
2204 ibin=i-1;
2205 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2206 break;
2207 }
2208 wSD=w[ibin];
2209 return kTRUE;
2210 }
800be8b7 2211
800be8b7 2212
c5dfa3e4 2213 else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2214const Int_t nbin=400;
2215Double_t bin[]={
22161.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
22174.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
22187.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
221910.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
222013.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
222115.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
222218.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
222321.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
222424.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
222527.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
222630.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
222733.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
222836.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
222939.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
223042.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
223145.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
223248.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
223351.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
223454.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
223557.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
223660.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
223763.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
223866.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
223969.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
224072.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
224175.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
224278.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
224381.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
224484.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
224587.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
224690.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
224793.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
224896.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
224999.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2250102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2251105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2252108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2253111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2254114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2255117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2256120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2257123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2258126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2259129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2260132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2261135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2262138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2263141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2264144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2265147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2266150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2267153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2268156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2269159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2270162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2271165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2272168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2273171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2274174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2275177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2276180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2277183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2278186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2279189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2280192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2281195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2282198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2283Double_t w[]={
22841.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
22850.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
22860.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
22870.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
22880.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
22890.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
22900.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
22910.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
22920.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
22930.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
22940.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
22950.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
22960.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
22970.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
22980.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
22990.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
23000.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
23010.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
23020.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
23030.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
23040.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
23050.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
23060.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
23070.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
23080.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
23090.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
23100.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
23110.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
23120.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
23130.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
23140.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
23150.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
23160.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
23170.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
23180.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
23190.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
23200.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
23210.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
23220.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
23230.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
23240.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
23250.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
23260.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
23270.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
23280.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
23290.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
23300.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
23310.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
23320.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
23330.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
23340.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
23350.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
23360.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
23370.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
23380.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
23390.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
23400.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
23410.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
23420.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
23430.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
23440.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
23450.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
23460.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
23470.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
23480.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
23490.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
23500.175006, 0.223482, 0.202706, 0.218108};
2351wDD = 0.207705;
2352wND = 0.289628;
2353wSD = -1;
2354
2355 Mmin = bin[0];
2356 Mmax = bin[nbin];
2357
2358 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2359
c5dfa3e4 2360 Int_t ibin=nbin-1;
2361 for(Int_t i=1; i<=nbin; i++)
2362 if(M<=bin[i]) {
2363 ibin=i-1;
2364 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2365 break;
2366 }
2367 wSD=w[ibin];
800be8b7 2368 return kTRUE;
c5dfa3e4 2369 }
2370
2371 return kFALSE;
800be8b7 2372}
06956108 2373
2374Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2375{
2376// Check if this is a heavy flavor decay product
2377 TParticle * part = (TParticle *) fParticles.At(ipart);
2378 Int_t mpdg = TMath::Abs(part->GetPdgCode());
2379 Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2380 //
2381 // Light hadron
2382 if (mfl >= 4 && mfl < 6) return kTRUE;
2383
2384 Int_t imo = part->GetFirstMother()-1;
2385 TParticle* pm = part;
2386 //
2387 // Heavy flavor hadron produced by generator
2388 while (imo > -1) {
2389 pm = (TParticle*)fParticles.At(imo);
2390 mpdg = TMath::Abs(pm->GetPdgCode());
2391 mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2392 if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2393 imo = pm->GetFirstMother()-1;
2394 }
2395 return kFALSE;
2396}
40fe59d4 2397
e81f2bac 2398Bool_t AliGenPythia::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
40fe59d4 2399{
2400 // check the eta/phi correspond to the detectors acceptance
2401 // iparticle is the index of the particle to be checked, for PHOS rotation case
2402 if (fCheckPHOS && IsInPHOS (phi,eta,iparticle)) return kTRUE;
2403 else if(fCheckEMCAL && IsInEMCAL (phi,eta)) return kTRUE;
2404 else if(fCheckBarrel && IsInBarrel( eta)) return kTRUE;
2405 else return kFALSE;
2406}
2407
e81f2bac 2408Bool_t AliGenPythia::TriggerOnSelectedParticles(Int_t np)
40fe59d4 2409{
2410 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2411 // implemented primaryly for kPyJets, but extended to any kind of process.
2412
2413 //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2414 // fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
2415
2416 Bool_t ok = kFALSE;
2417 for (Int_t i=0; i< np; i++) {
2418
2419 TParticle* iparticle = (TParticle *) fParticles.At(i);
2420
2421 Int_t pdg = iparticle->GetPdgCode();
2422 Int_t status = iparticle->GetStatusCode();
2423 Int_t imother = iparticle->GetFirstMother() - 1;
2424
2425 TParticle* pmother = 0x0;
2426 Int_t momStatus = -1;
2427 Int_t momPdg = -1;
2428 if(imother > 0 ){
2429 pmother = (TParticle *) fParticles.At(imother);
2430 momStatus = pmother->GetStatusCode();
2431 momPdg = pmother->GetPdgCode();
2432 }
2433
2434 ok = kFALSE;
2435
2436 //
2437 // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2438 //
2439 // Hadron
2440 if (fHadronInCalo && status == 1)
2441 {
2442 if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
2443 // (in case neutral mesons were declared stable)
2444 ok = kTRUE;
2445 }
2446 //Electron
2447 else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2448 {
2449 ok = kTRUE;
2450 }
2451 //Fragmentation photon
2452 else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2453 {
2454 if(momStatus != 11) ok = kTRUE ; // No photon from hadron decay
2455 }
2456 // Decay photon
2457 else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2458 {
2459 if( momStatus == 11)
2460 {
2461 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2462 // pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2463 ok = kTRUE ; // photon from hadron decay
2464
2465 //In case only decays from pi0 or eta requested
2466 if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2467 if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2468 }
2469
2470 }
2471 // Pi0 or Eta particle
2472 else if ((fPi0InCalo || fEtaInCalo))
2473 {
2474 if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2475
2476 if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2477 {
2478 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2479 ok = kTRUE;
2480 }
2481 else if (fEtaInCalo && pdg == 221)
2482 {
2483 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2484 ok = kTRUE;
2485 }
2486
2487 }// pi0 or eta
2488
2489 //
2490 // Check that the selected particle is in the calorimeter acceptance
2491 //
2492 if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2493 {
2494 //Just check if the selected particle falls in the acceptance
2495 if(!fForceNeutralMeson2PhotonDecay )
2496 {
2497 //printf("\t Check acceptance! \n");
2498 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2499 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2500
2501 if(CheckDetectorAcceptance(phi,eta,i))
2502 {
2503 ok =kTRUE;
2504 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));
2505 //printf("\t Accept \n");
2506 break;
2507 }
2508 else ok = kFALSE;
2509 }
2510 //Mesons have several decay modes, select only those decaying into 2 photons
2511 else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2512 {
2513 // In case we want the pi0/eta trigger,
2514 // check the decay mode (2 photons)
2515
2516 //printf("\t Force decay 2 gamma\n");
2517
2518 Int_t ndaughters = iparticle->GetNDaughters();
2519 if(ndaughters != 2){
2520 ok=kFALSE;
2521 continue;
2522 }
2523
2524 TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2525 TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2526 if(!d1 || !d2) {
2527 ok=kFALSE;
2528 continue;
2529 }
2530
2531 //iparticle->Print();
2532 //d1->Print();
2533 //d2->Print();
2534
2535 Int_t pdgD1 = d1->GetPdgCode();
2536 Int_t pdgD2 = d2->GetPdgCode();
2537 //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2538 //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2539
2540 if(pdgD1 != 22 || pdgD2 != 22){
2541 ok = kFALSE;
2542 continue;
2543 }
2544
2545 //printf("\t accept decay\n");
2546
2547 //Trigger on the meson, not on the daughter
2548 if(!fDecayPhotonInCalo){
2549
2550 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2551 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2552
2553 if(CheckDetectorAcceptance(phi,eta,i))
2554 {
2555 //printf("\t Accept meson pdg %d\n",pdg);
2556 ok =kTRUE;
2557 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));
2558 break;
2559 } else {
2560 ok=kFALSE;
2561 continue;
2562 }
2563 }
2564
2565 //printf("Check daughters acceptance\n");
2566
2567 //Trigger on the meson daughters
2568 //Photon 1
2569 Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2570 Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
2571 if(d1->Pt() > fTriggerParticleMinPt)
2572 {
2573 //printf("\t Check acceptance photon 1! \n");
2574 if(CheckDetectorAcceptance(phi,eta,i))
2575 {
2576 //printf("\t Accept Photon 1\n");
2577 ok =kTRUE;
2578 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));
2579 break;
2580 }
2581 else ok = kFALSE;
2582 } // pt cut
2583 else ok = kFALSE;
2584
2585 //Photon 2
2586 phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2587 eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2588
2589 if(d2->Pt() > fTriggerParticleMinPt)
2590 {
2591 //printf("\t Check acceptance photon 2! \n");
2592 if(CheckDetectorAcceptance(phi,eta,i))
2593 {
2594 //printf("\t Accept Photon 2\n");
2595 ok =kTRUE;
2596 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));
2597 break;
2598 }
2599 else ok = kFALSE;
2600 } // pt cut
2601 else ok = kFALSE;
2602 } // force 2 photon daughters in pi0/eta decays
2603 else ok = kFALSE;
2604 } else ok = kFALSE; // check acceptance
2605 } // primary loop
2606
2607 //
2608 // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2609 // A particle passing all trigger conditions except phi position in PHOS, is used as reference
2610 //
2611 if(fCheckPHOSeta)
2612 {
2613 RotatePhi(ok);
2614 }
2615
2616 return ok;
2617}
2618