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