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