]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA6/AliGenPythia.cxx
Renamed local method rnorml to rnormlEP to avoid conflict with
[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
e28ccdaf 382 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
383 if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
384
8d2cd130 385 // Fragmentation?
386 if (fFragmentation) {
387 fPythia->SetMSTP(111,1);
388 } else {
389 fPythia->SetMSTP(111,0);
390 }
391
392
393// initial state radiation
394 fPythia->SetMSTP(61,fGinit);
395// final state radiation
396 fPythia->SetMSTP(71,fGfinal);
397// pt - kick
398 if (fPtKick > 0.) {
399 fPythia->SetMSTP(91,1);
688950ef 400 fPythia->SetPARP(91,fPtKick);
401 fPythia->SetPARP(93, 4. * fPtKick);
8d2cd130 402 } else {
403 fPythia->SetMSTP(91,0);
404 }
405
5fa4b20b 406
407 if (fReadFromFile) {
777e69b0 408 fRL = AliRunLoader::Open(fkFileName, "Partons");
5fa4b20b 409 fRL->LoadKinematics();
410 fRL->LoadHeader();
411 } else {
412 fRL = 0x0;
413 }
fdea4387 414 //
efe3b1cd 415 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
fdea4387 416 // Forward Paramters to the AliPythia object
417 fDecayer->SetForceDecay(fForceDecay);
beac474c 418// Switch off Heavy Flavors on request
419 if (fHFoff) {
51387949 420 // Maximum number of quark flavours used in pdf
beac474c 421 fPythia->SetMSTP(58, 3);
51387949 422 // Maximum number of flavors that can be used in showers
beac474c 423 fPythia->SetMSTJ(45, 3);
51387949 424 // Switch off g->QQbar splitting in decay table
425 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
beac474c 426 }
fdea4387 427
51387949 428 fDecayer->Init();
429
8d2cd130 430
431// Parent and Children Selection
432 switch (fProcess)
433 {
65f2626c 434 case kPyOldUEQ2ordered:
435 case kPyOldUEQ2ordered2:
436 case kPyOldPopcorn:
437 break;
8d2cd130 438 case kPyCharm:
439 case kPyCharmUnforced:
adf4d898 440 case kPyCharmPbPbMNR:
aabc7187 441 case kPyCharmpPbMNR:
e0e89f40 442 case kPyCharmppMNR:
443 case kPyCharmppMNRwmi:
8d2cd130 444 fParentSelect[0] = 411;
445 fParentSelect[1] = 421;
446 fParentSelect[2] = 431;
447 fParentSelect[3] = 4122;
9ccc3d0e 448 fParentSelect[4] = 4232;
449 fParentSelect[5] = 4132;
450 fParentSelect[6] = 4332;
8d2cd130 451 fFlavorSelect = 4;
452 break;
adf4d898 453 case kPyD0PbPbMNR:
454 case kPyD0pPbMNR:
455 case kPyD0ppMNR:
8d2cd130 456 fParentSelect[0] = 421;
457 fFlavorSelect = 4;
458 break;
90d7b703 459 case kPyDPlusPbPbMNR:
460 case kPyDPluspPbMNR:
461 case kPyDPlusppMNR:
462 fParentSelect[0] = 411;
463 fFlavorSelect = 4;
464 break;
e0e89f40 465 case kPyDPlusStrangePbPbMNR:
466 case kPyDPlusStrangepPbMNR:
467 case kPyDPlusStrangeppMNR:
468 fParentSelect[0] = 431;
469 fFlavorSelect = 4;
470 break;
68504d42 471 case kPyLambdacppMNR:
472 fParentSelect[0] = 4122;
473 fFlavorSelect = 4;
474 break;
8d2cd130 475 case kPyBeauty:
9dfe63b3 476 case kPyBeautyJets:
adf4d898 477 case kPyBeautyPbPbMNR:
478 case kPyBeautypPbMNR:
479 case kPyBeautyppMNR:
e0e89f40 480 case kPyBeautyppMNRwmi:
8d2cd130 481 fParentSelect[0]= 511;
482 fParentSelect[1]= 521;
483 fParentSelect[2]= 531;
484 fParentSelect[3]= 5122;
485 fParentSelect[4]= 5132;
486 fParentSelect[5]= 5232;
487 fParentSelect[6]= 5332;
488 fFlavorSelect = 5;
489 break;
490 case kPyBeautyUnforced:
491 fParentSelect[0] = 511;
492 fParentSelect[1] = 521;
493 fParentSelect[2] = 531;
494 fParentSelect[3] = 5122;
495 fParentSelect[4] = 5132;
496 fParentSelect[5] = 5232;
497 fParentSelect[6] = 5332;
498 fFlavorSelect = 5;
499 break;
500 case kPyJpsiChi:
501 case kPyJpsi:
502 fParentSelect[0] = 443;
503 break;
f529e69b 504 case kPyMbDefault:
0bd3d7c5 505 case kPyMbAtlasTuneMC09:
8d2cd130 506 case kPyMb:
04081a91 507 case kPyMbWithDirectPhoton:
8d2cd130 508 case kPyMbNonDiffr:
d7de4a67 509 case kPyMbMSEL1:
8d2cd130 510 case kPyJets:
511 case kPyDirectGamma:
e33acb67 512 case kPyLhwgMb:
8d2cd130 513 break;
589380c6 514 case kPyW:
0f6ee828 515 case kPyZ:
9a8774a1 516 case kPyMBRSingleDiffraction:
517 case kPyMBRDoubleDiffraction:
518 case kPyMBRCentralDiffraction:
589380c6 519 break;
8d2cd130 520 }
521//
592f8307 522//
523// JetFinder for Trigger
524//
525// Configure detector (EMCAL like)
526//
d7de4a67 527 fPythia->SetPARU(51, fPycellEtaMax);
528 fPythia->SetMSTU(51, fPycellNEta);
529 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 530//
531// Configure Jet Finder
532//
d7de4a67 533 fPythia->SetPARU(58, fPycellThreshold);
534 fPythia->SetPARU(52, fPycellEtSeed);
535 fPythia->SetPARU(53, fPycellMinEtJet);
536 fPythia->SetPARU(54, fPycellMaxRadius);
537 fPythia->SetMSTU(54, 2);
592f8307 538//
8d2cd130 539// This counts the total number of calls to Pyevnt() per run.
540 fTrialsRun = 0;
541 fQ = 0.;
542 fX1 = 0.;
543 fX2 = 0.;
544 fNev = 0 ;
545//
1d568bc2 546//
547//
8d2cd130 548 AliGenMC::Init();
1d568bc2 549//
550//
551//
552 if (fSetNuclei) {
553 fDyBoost = 0;
554 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
555 }
d7de4a67 556
cd07c39b 557 fPythia->SetPARJ(200, 0.0);
558 fPythia->SetPARJ(199, 0.0);
559 fPythia->SetPARJ(198, 0.0);
560 fPythia->SetPARJ(197, 0.0);
561
562 if (fQuench == 1) {
5fa4b20b 563 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 564 }
3a709cfa 565
66b8652c 566 if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
567
b976f7d8 568 if (fQuench == 3) {
569 // Nestor's change of the splittings
570 fPythia->SetPARJ(200, 0.8);
571 fPythia->SetMSTJ(41, 1); // QCD radiation only
572 fPythia->SetMSTJ(42, 2); // angular ordering
573 fPythia->SetMSTJ(44, 2); // option to run alpha_s
574 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
575 fPythia->SetMSTJ(50, 0); // No coherence in first branching
576 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
cd07c39b 577 } else if (fQuench == 4) {
578 // Armesto-Cunqueiro-Salgado change of the splittings.
579 AliFastGlauber* glauber = AliFastGlauber::Instance();
580 glauber->Init(2);
581 //read and store transverse almonds corresponding to differnt
582 //impact parameters.
cd07c39b 583 glauber->SetCentralityClass(0.,0.1);
584 fPythia->SetPARJ(200, 1.);
585 fPythia->SetPARJ(198, fQhat);
586 fPythia->SetPARJ(199, fLength);
cd07c39b 587 fPythia->SetMSTJ(42, 2); // angular ordering
588 fPythia->SetMSTJ(44, 2); // option to run alpha_s
cd07c39b 589 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
b976f7d8 590 }
8d2cd130 591}
592
593void AliGenPythia::Generate()
594{
595// Generate one event
19ef8cf0 596 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
8d2cd130 597 fDecayer->ForceDecay();
598
13cca24a 599 Double_t polar[3] = {0,0,0};
600 Double_t origin[3] = {0,0,0};
601 Double_t p[4];
8d2cd130 602// converts from mm/c to s
603 const Float_t kconv=0.001/2.999792458e8;
604//
605 Int_t nt=0;
606 Int_t jev=0;
607 Int_t j, kf;
608 fTrials=0;
f913ec4f 609 fEventTime = 0.;
610
2590ccf9 611
8d2cd130 612
613 // Set collision vertex position
2590ccf9 614 if (fVertexSmear == kPerEvent) Vertex();
615
8d2cd130 616// event loop
617 while(1)
618 {
32d6ef7d 619//
5fa4b20b 620// Produce event
32d6ef7d 621//
32d6ef7d 622//
5fa4b20b 623// Switch hadronisation off
624//
625 fPythia->SetMSTJ(1, 0);
cd07c39b 626
627 if (fQuench ==4){
628 Double_t bimp;
629 // Quenching comes through medium-modified splitting functions.
630 AliFastGlauber::Instance()->GetRandomBHard(bimp);
e6fe9b82 631 fPythia->SetPARJ(197, bimp);
632 fImpact = bimp;
6c43eccb 633 fPythia->Qpygin0();
cd07c39b 634 }
32d6ef7d 635//
5fa4b20b 636// Either produce new event or read partons from file
637//
638 if (!fReadFromFile) {
beac474c 639 if (!fNewMIS) {
640 fPythia->Pyevnt();
641 } else {
642 fPythia->Pyevnw();
643 }
5fa4b20b 644 fNpartons = fPythia->GetN();
645 } else {
33c3c91a 646 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
647 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
5fa4b20b 648 fPythia->SetN(0);
649 LoadEvent(fRL->Stack(), 0 , 1);
650 fPythia->Pyedit(21);
651 }
652
32d6ef7d 653//
654// Run quenching routine
655//
5fa4b20b 656 if (fQuench == 1) {
657 fPythia->Quench();
658 } else if (fQuench == 2){
659 fPythia->Pyquen(208., 0, 0.);
b976f7d8 660 } else if (fQuench == 3) {
661 // Quenching is via multiplicative correction of the splittings
5fa4b20b 662 }
b976f7d8 663
32d6ef7d 664//
5fa4b20b 665// Switch hadronisation on
32d6ef7d 666//
20e47f08 667 if (fHadronisation) {
668 fPythia->SetMSTJ(1, 1);
5fa4b20b 669//
670// .. and perform hadronisation
aea21c57 671// printf("Calling hadronisation %d\n", fPythia->GetN());
03358a32 672
673 if (fPatchOmegaDalitz) {
674 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
675 fPythia->Pyexec();
676 fPythia->DalitzDecays();
677 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
678 }
679 fPythia->Pyexec();
20e47f08 680 }
8d2cd130 681 fTrials++;
8507138f 682 fPythia->ImportParticles(&fParticles,"All");
03358a32 683
df275cfa 684 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
e1cf8975 685 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
f64640b5 686
8d2cd130 687//
688//
689//
690 Int_t i;
691
dbd64db6 692 fNprimaries = 0;
8507138f 693 Int_t np = fParticles.GetEntriesFast();
5fa4b20b 694
7c21f297 695 if (np == 0) continue;
8d2cd130 696//
2590ccf9 697
8d2cd130 698//
699 Int_t* pParent = new Int_t[np];
700 Int_t* pSelected = new Int_t[np];
701 Int_t* trackIt = new Int_t[np];
5fa4b20b 702 for (i = 0; i < np; i++) {
8d2cd130 703 pParent[i] = -1;
704 pSelected[i] = 0;
705 trackIt[i] = 0;
706 }
707
708 Int_t nc = 0; // Total n. of selected particles
709 Int_t nParents = 0; // Selected parents
710 Int_t nTkbles = 0; // Trackable particles
f529e69b 711 if (fProcess != kPyMbDefault &&
712 fProcess != kPyMb &&
6d2ec66d 713 fProcess != kPyMbAtlasTuneMC09 &&
04081a91 714 fProcess != kPyMbWithDirectPhoton &&
f529e69b 715 fProcess != kPyJets &&
8d2cd130 716 fProcess != kPyDirectGamma &&
589380c6 717 fProcess != kPyMbNonDiffr &&
d7de4a67 718 fProcess != kPyMbMSEL1 &&
f529e69b 719 fProcess != kPyW &&
720 fProcess != kPyZ &&
721 fProcess != kPyCharmppMNRwmi &&
9dfe63b3 722 fProcess != kPyBeautyppMNRwmi &&
723 fProcess != kPyBeautyJets) {
8d2cd130 724
5fa4b20b 725 for (i = 0; i < np; i++) {
8507138f 726 TParticle* iparticle = (TParticle *) fParticles.At(i);
8d2cd130 727 Int_t ks = iparticle->GetStatusCode();
728 kf = CheckPDGCode(iparticle->GetPdgCode());
729// No initial state partons
730 if (ks==21) continue;
731//
732// Heavy Flavor Selection
733//
734 // quark ?
735 kf = TMath::Abs(kf);
736 Int_t kfl = kf;
9ff6c04c 737 // Resonance
f913ec4f 738
9ff6c04c 739 if (kfl > 100000) kfl %= 100000;
183a5ca9 740 if (kfl > 10000) kfl %= 10000;
8d2cd130 741 // meson ?
742 if (kfl > 10) kfl/=100;
743 // baryon
744 if (kfl > 10) kfl/=10;
8d2cd130 745 Int_t ipa = iparticle->GetFirstMother()-1;
746 Int_t kfMo = 0;
f913ec4f 747//
748// Establish mother daughter relation between heavy quarks and mesons
749//
750 if (kf >= fFlavorSelect && kf <= 6) {
751 Int_t idau = iparticle->GetFirstDaughter() - 1;
752 if (idau > -1) {
8507138f 753 TParticle* daughter = (TParticle *) fParticles.At(idau);
f913ec4f 754 Int_t pdgD = daughter->GetPdgCode();
755 if (pdgD == 91 || pdgD == 92) {
756 Int_t jmin = daughter->GetFirstDaughter() - 1;
757 Int_t jmax = daughter->GetLastDaughter() - 1;
2ab330c9 758 for (Int_t jp = jmin; jp <= jmax; jp++)
759 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
f913ec4f 760 } // is string or cluster
761 } // has daughter
762 } // heavy quark
8d2cd130 763
f913ec4f 764
8d2cd130 765 if (ipa > -1) {
8507138f 766 TParticle * mother = (TParticle *) fParticles.At(ipa);
8d2cd130 767 kfMo = TMath::Abs(mother->GetPdgCode());
768 }
f913ec4f 769
8d2cd130 770 // What to keep in Stack?
771 Bool_t flavorOK = kFALSE;
772 Bool_t selectOK = kFALSE;
773 if (fFeedDownOpt) {
32d6ef7d 774 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 775 } else {
32d6ef7d 776 if (kfl > fFlavorSelect) {
777 nc = -1;
778 break;
779 }
780 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 781 }
782 switch (fStackFillOpt) {
06956108 783 case kHeavyFlavor:
8d2cd130 784 case kFlavorSelection:
32d6ef7d 785 selectOK = kTRUE;
786 break;
8d2cd130 787 case kParentSelection:
32d6ef7d 788 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
789 break;
8d2cd130 790 }
791 if (flavorOK && selectOK) {
792//
793// Heavy flavor hadron or quark
794//
795// Kinematic seletion on final state heavy flavor mesons
796 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
797 {
9ff6c04c 798 continue;
8d2cd130 799 }
800 pSelected[i] = 1;
801 if (ParentSelected(kf)) ++nParents; // Update parent count
802// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
803 } else {
804// Kinematic seletion on decay products
805 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 806 && !KinematicSelection(iparticle, 1))
8d2cd130 807 {
9ff6c04c 808 continue;
8d2cd130 809 }
810//
811// Decay products
812// Select if mother was selected and is not tracked
813
814 if (pSelected[ipa] &&
815 !trackIt[ipa] && // mother will be tracked ?
816 kfMo != 5 && // mother is b-quark, don't store fragments
817 kfMo != 4 && // mother is c-quark, don't store fragments
818 kf != 92) // don't store string
819 {
820//
821// Semi-stable or de-selected: diselect decay products:
822//
823//
824 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
825 {
826 Int_t ipF = iparticle->GetFirstDaughter();
827 Int_t ipL = iparticle->GetLastDaughter();
828 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
829 }
830// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
831 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
832 }
833 }
834 if (pSelected[i] == -1) pSelected[i] = 0;
835 if (!pSelected[i]) continue;
836 // Count quarks only if you did not include fragmentation
837 if (fFragmentation && kf <= 10) continue;
9ff6c04c 838
8d2cd130 839 nc++;
840// Decision on tracking
841 trackIt[i] = 0;
842//
843// Track final state particle
844 if (ks == 1) trackIt[i] = 1;
845// Track semi-stable particles
d25cfd65 846 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 847// Track particles selected by process if undecayed.
848 if (fForceDecay == kNoDecay) {
849 if (ParentSelected(kf)) trackIt[i] = 1;
850 } else {
851 if (ParentSelected(kf)) trackIt[i] = 0;
852 }
853 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
854//
855//
856
857 } // particle selection loop
858 if (nc > 0) {
859 for (i = 0; i<np; i++) {
860 if (!pSelected[i]) continue;
8507138f 861 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 862 kf = CheckPDGCode(iparticle->GetPdgCode());
863 Int_t ks = iparticle->GetStatusCode();
864 p[0] = iparticle->Px();
865 p[1] = iparticle->Py();
866 p[2] = iparticle->Pz();
a920faf9 867 p[3] = iparticle->Energy();
868
2590ccf9 869 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
870 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
871 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
872
21391258 873 Float_t tof = fTime + kconv*iparticle->T();
8d2cd130 874 Int_t ipa = iparticle->GetFirstMother()-1;
875 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 876
877 PushTrack(fTrackIt*trackIt[i], iparent, kf,
878 p[0], p[1], p[2], p[3],
879 origin[0], origin[1], origin[2], tof,
880 polar[0], polar[1], polar[2],
881 kPPrimary, nt, 1., ks);
8d2cd130 882 pParent[i] = nt;
dbd64db6 883 KeepTrack(nt);
884 fNprimaries++;
642f15cf 885 } // PushTrack loop
8d2cd130 886 }
887 } else {
888 nc = GenerateMB();
889 } // mb ?
f913ec4f 890
891 GetSubEventTime();
8d2cd130 892
234f6d32 893 delete[] pParent;
894 delete[] pSelected;
895 delete[] trackIt;
8d2cd130 896
897 if (nc > 0) {
898 switch (fCountMode) {
899 case kCountAll:
900 // printf(" Count all \n");
901 jev += nc;
902 break;
903 case kCountParents:
904 // printf(" Count parents \n");
905 jev += nParents;
906 break;
907 case kCountTrackables:
908 // printf(" Count trackable \n");
909 jev += nTkbles;
910 break;
911 }
912 if (jev >= fNpart || fNpart == -1) {
913 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 914
8d2cd130 915 fQ += fPythia->GetVINT(51);
916 fX1 += fPythia->GetVINT(41);
917 fX2 += fPythia->GetVINT(42);
918 fTrialsRun += fTrials;
919 fNev++;
920 MakeHeader();
921 break;
922 }
923 }
924 } // event loop
925 SetHighWaterMark(nt);
926// adjust weight due to kinematic selection
b88f5cea 927// AdjustWeights();
8d2cd130 928// get cross-section
929 fXsection=fPythia->GetPARI(1);
930}
931
932Int_t AliGenPythia::GenerateMB()
933{
934//
935// Min Bias selection and other global selections
936//
937 Int_t i, kf, nt, iparent;
938 Int_t nc = 0;
13cca24a 939 Double_t p[4];
940 Double_t polar[3] = {0,0,0};
941 Double_t origin[3] = {0,0,0};
8d2cd130 942// converts from mm/c to s
943 const Float_t kconv=0.001/2.999792458e8;
944
e0e89f40 945
8507138f 946 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
aea21c57 947
5fa4b20b 948
e0e89f40 949
8d2cd130 950 Int_t* pParent = new Int_t[np];
951 for (i=0; i< np; i++) pParent[i] = -1;
0e9e66f0 952 if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
953 && fEtMinJet > 0.) {
bed3df2c 954 TParticle* jet1 = (TParticle *) fParticles.At(6);
8507138f 955 TParticle* jet2 = (TParticle *) fParticles.At(7);
bed3df2c 956
957 if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
234f6d32 958 delete [] pParent;
959 return 0;
960 }
8d2cd130 961 }
e0e89f40 962
40fe59d4 963
964 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
ac01c7c6 965 // implemented primaryly for kPyJets, but extended to any kind of process.
40fe59d4 966 if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
967 (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
968 Bool_t ok = TriggerOnSelectedParticles(np);
43e3b80a 969
970 if(!ok) {
971 delete[] pParent;
972 return 0;
ec2c406e 973 }
43e3b80a 974 }
9dfe63b3 975
800be8b7 976 // Check for diffraction
977 if(fkTuneForDiff) {
c5dfa3e4 978 if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
800be8b7 979 if(!CheckDiffraction()) {
980 delete [] pParent;
981 return 0;
982 }
983 }
984 }
985
700b9416 986 // Check for minimum multiplicity
987 if (fTriggerMultiplicity > 0) {
988 Int_t multiplicity = 0;
989 for (i = 0; i < np; i++) {
8507138f 990 TParticle * iparticle = (TParticle *) fParticles.At(i);
700b9416 991
992 Int_t statusCode = iparticle->GetStatusCode();
993
994 // Initial state particle
e296848e 995 if (statusCode != 1)
700b9416 996 continue;
38112f3f 997 // eta cut
700b9416 998 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
999 continue;
38112f3f 1000 // pt cut
1001 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1002 continue;
1003
700b9416 1004 TParticlePDG* pdgPart = iparticle->GetPDG();
1005 if (pdgPart && pdgPart->Charge() == 0)
1006 continue;
1007
1008 ++multiplicity;
1009 }
1010
1011 if (multiplicity < fTriggerMultiplicity) {
1012 delete [] pParent;
1013 return 0;
1014 }
38112f3f 1015 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
700b9416 1016 }
90a236ce 1017
90a236ce 1018
7ce1321b 1019 if (fTriggerParticle) {
1020 Bool_t triggered = kFALSE;
1021 for (i = 0; i < np; i++) {
8507138f 1022 TParticle * iparticle = (TParticle *) fParticles.At(i);
7ce1321b 1023 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 1024 if (kf != fTriggerParticle) continue;
7ce1321b 1025 if (iparticle->Pt() == 0.) continue;
1026 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
2591bd0e 1027 if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
7ce1321b 1028 triggered = kTRUE;
1029 break;
1030 }
234f6d32 1031 if (!triggered) {
1032 delete [] pParent;
1033 return 0;
1034 }
7ce1321b 1035 }
e0e89f40 1036
1037
1038 // Check if there is a ccbar or bbbar pair with at least one of the two
1039 // in fYMin < y < fYMax
2f405d65 1040
9dfe63b3 1041 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
9ccc3d0e 1042 TParticle *partCheck;
1043 TParticle *mother;
e0e89f40 1044 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 1045 Bool_t theChild=kFALSE;
5d76923e 1046 Bool_t triggered=kFALSE;
9ccc3d0e 1047 Float_t y;
1048 Int_t pdg,mpdg,mpdgUpperFamily;
e0e89f40 1049 for(i=0; i<np; i++) {
9ccc3d0e 1050 partCheck = (TParticle*)fParticles.At(i);
1051 pdg = partCheck->GetPdgCode();
1052 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1053 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1054 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1055 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1056 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1057 }
5d76923e 1058 if(fTriggerParticle) {
1059 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1060 }
9ccc3d0e 1061 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1062 Int_t mi = partCheck->GetFirstMother() - 1;
1063 if(mi<0) continue;
1064 mother = (TParticle*)fParticles.At(mi);
1065 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 1066 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 1067 if ( ParentSelected(mpdg) ||
1068 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1069 if (KinematicSelection(partCheck,1)) {
1070 theChild=kTRUE;
1071 }
1072 }
1073 }
e0e89f40 1074 }
9ccc3d0e 1075 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
234f6d32 1076 delete[] pParent;
e0e89f40 1077 return 0;
1078 }
9ccc3d0e 1079 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1080 delete[] pParent;
1081 return 0;
1082 }
5d76923e 1083 if(fTriggerParticle && !triggered) { // particle requested is not produced
1084 delete[] pParent;
1085 return 0;
1086 }
9ccc3d0e 1087
e0e89f40 1088 }
aea21c57 1089
0f6ee828 1090 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
f529e69b 1091 if ( (fProcess == kPyW ||
1092 fProcess == kPyZ ||
1093 fProcess == kPyMbDefault ||
1094 fProcess == kPyMb ||
6d2ec66d 1095 fProcess == kPyMbAtlasTuneMC09 ||
04081a91 1096 fProcess == kPyMbWithDirectPhoton ||
f529e69b 1097 fProcess == kPyMbNonDiffr)
0f6ee828 1098 && (fCutOnChild == 1) ) {
1099 if ( !CheckKinematicsOnChild() ) {
234f6d32 1100 delete[] pParent;
0f6ee828 1101 return 0;
1102 }
aea21c57 1103 }
1104
1105
f913ec4f 1106 for (i = 0; i < np; i++) {
8d2cd130 1107 Int_t trackIt = 0;
8507138f 1108 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 1109 kf = CheckPDGCode(iparticle->GetPdgCode());
1110 Int_t ks = iparticle->GetStatusCode();
1111 Int_t km = iparticle->GetFirstMother();
06956108 1112 if (
1113 (((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
1114 ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
1115 )
1116 {
1117 nc++;
8d2cd130 1118 if (ks == 1) trackIt = 1;
1119 Int_t ipa = iparticle->GetFirstMother()-1;
1120
1121 iparent = (ipa > -1) ? pParent[ipa] : -1;
1122
1123//
1124// store track information
1125 p[0] = iparticle->Px();
1126 p[1] = iparticle->Py();
1127 p[2] = iparticle->Pz();
a920faf9 1128 p[3] = iparticle->Energy();
1406f599 1129
a920faf9 1130
2590ccf9 1131 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1132 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1133 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1134
21391258 1135 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
a920faf9 1136
1137 PushTrack(fTrackIt*trackIt, iparent, kf,
1138 p[0], p[1], p[2], p[3],
1139 origin[0], origin[1], origin[2], tof,
1140 polar[0], polar[1], polar[2],
1141 kPPrimary, nt, 1., ks);
dbd64db6 1142 fNprimaries++;
8d2cd130 1143 KeepTrack(nt);
1144 pParent[i] = nt;
f913ec4f 1145 SetHighWaterMark(nt);
1146
8d2cd130 1147 } // select particle
1148 } // particle loop
1149
234f6d32 1150 delete[] pParent;
e0e89f40 1151
f913ec4f 1152 return 1;
8d2cd130 1153}
1154
1155
1156void AliGenPythia::FinishRun()
1157{
1158// Print x-section summary
1159 fPythia->Pystat(1);
2779fc64 1160
1161 if (fNev > 0.) {
1162 fQ /= fNev;
1163 fX1 /= fNev;
1164 fX2 /= fNev;
1165 }
1166
8d2cd130 1167 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1168 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1169}
1170
7184e472 1171void AliGenPythia::AdjustWeights() const
8d2cd130 1172{
1173// Adjust the weights after generation of all events
1174//
e2bddf81 1175 if (gAlice) {
1176 TParticle *part;
1177 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1178 for (Int_t i=0; i<ntrack; i++) {
1179 part= gAlice->GetMCApp()->Particle(i);
1180 part->SetWeight(part->GetWeight()*fKineBias);
1181 }
8d2cd130 1182 }
1183}
1184
20e47f08 1185void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
8d2cd130 1186{
1187// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1188
1a626d4e 1189 fAProjectile = a1;
1190 fATarget = a2;
66f02a7f 1191 fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
1d568bc2 1192 fSetNuclei = kTRUE;
8d2cd130 1193}
1194
1195
1196void AliGenPythia::MakeHeader()
1197{
7184e472 1198//
1199// Make header for the simulated event
1200//
183a5ca9 1201 if (gAlice) {
1202 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1203 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1204 }
1205
8d2cd130 1206// Builds the event header, to be called after each event
e5c87a3d 1207 if (fHeader) delete fHeader;
1208 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1209//
1210// Event type
800be8b7 1211 if(fProcDiff>0){
1212 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1213 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1214 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1215 }
1216 else
e5c87a3d 1217 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1218//
1219// Number of trials
e5c87a3d 1220 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1221//
1222// Event Vertex
d25cfd65 1223 fHeader->SetPrimaryVertex(fVertex);
21391258 1224 fHeader->SetInteractionTime(fTime+fEventTime);
dbd64db6 1225//
1226// Number of primaries
1227 fHeader->SetNProduced(fNprimaries);
8d2cd130 1228//
1229// Jets that have triggered
f913ec4f 1230
9dfe63b3 1231 //Need to store jets for b-jet studies too!
2f405d65 1232 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
8d2cd130 1233 {
1234 Int_t ntrig, njet;
1235 Float_t jets[4][10];
1236 GetJets(njet, ntrig, jets);
9ff6c04c 1237
8d2cd130 1238
1239 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1240 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1241 jets[3][i]);
1242 }
1243 }
5fa4b20b 1244//
1245// Copy relevant information from external header, if present.
1246//
1247 Float_t uqJet[4];
1248
1249 if (fRL) {
1250 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1251 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1252 {
1253 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1254
1255
1256 exHeader->TriggerJet(i, uqJet);
1257 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1258 }
1259 }
1260//
1261// Store quenching parameters
1262//
1263 if (fQuench){
b6262a45 1264 Double_t z[4] = {0.};
1265 Double_t xp = 0.;
1266 Double_t yp = 0.;
1267
7c21f297 1268 if (fQuench == 1) {
1269 // Pythia::Quench()
1270 fPythia->GetQuenchingParameters(xp, yp, z);
1bab4b79 1271 } else if (fQuench == 2){
7c21f297 1272 // Pyquen
1273 Double_t r1 = PARIMP.rb1;
1274 Double_t r2 = PARIMP.rb2;
1275 Double_t b = PARIMP.b1;
1276 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1277 Double_t phi = PARIMP.psib1;
1278 xp = r * TMath::Cos(phi);
1279 yp = r * TMath::Sin(phi);
1280
1bab4b79 1281 } else if (fQuench == 4) {
1282 // QPythia
5831e75f 1283 Double_t xy[2];
e9719084 1284 Double_t i0i1[2];
1285 AliFastGlauber::Instance()->GetSavedXY(xy);
1286 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
5831e75f 1287 xp = xy[0];
1288 yp = xy[1];
e6fe9b82 1289 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
7c21f297 1290 }
1bab4b79 1291
7c21f297 1292 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1293 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1bab4b79 1294 }
beac474c 1295//
1296// Store pt^hard
1297 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 1298//
cf57b268 1299// Pass header
5fa4b20b 1300//
cf57b268 1301 AddHeader(fHeader);
4c4eac97 1302 fHeader = 0x0;
8d2cd130 1303}
cf57b268 1304
c2fc9d31 1305Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
8d2cd130 1306{
1307// Check the kinematic trigger condition
1308//
1309 Double_t eta[2];
1310 eta[0] = jet1->Eta();
1311 eta[1] = jet2->Eta();
1312 Double_t phi[2];
1313 phi[0] = jet1->Phi();
1314 phi[1] = jet2->Phi();
1315 Int_t pdg[2];
1316 pdg[0] = jet1->GetPdgCode();
1317 pdg[1] = jet2->GetPdgCode();
1318 Bool_t triggered = kFALSE;
1319
2f405d65 1320 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8d2cd130 1321 Int_t njets = 0;
1322 Int_t ntrig = 0;
1323 Float_t jets[4][10];
1324//
1325// Use Pythia clustering on parton level to determine jet axis
1326//
1327 GetJets(njets, ntrig, jets);
1328
76d6ba9a 1329 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1330//
1331 } else {
1332 Int_t ij = 0;
1333 Int_t ig = 1;
1334 if (pdg[0] == kGamma) {
1335 ij = 1;
1336 ig = 0;
1337 }
1338 //Check eta range first...
1339 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1340 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1341 {
1342 //Eta is okay, now check phi range
1343 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1344 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1345 {
1346 triggered = kTRUE;
1347 }
1348 }
1349 }
1350 return triggered;
1351}
aea21c57 1352
1353
aea21c57 1354
7184e472 1355Bool_t AliGenPythia::CheckKinematicsOnChild(){
1356//
1357//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1358//
aea21c57 1359 Bool_t checking = kFALSE;
1360 Int_t j, kcode, ks, km;
1361 Int_t nPartAcc = 0; //number of particles in the acceptance range
1362 Int_t numberOfAcceptedParticles = 1;
1363 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1364 Int_t npart = fParticles.GetEntriesFast();
aea21c57 1365
0f6ee828 1366 for (j = 0; j<npart; j++) {
8507138f 1367 TParticle * jparticle = (TParticle *) fParticles.At(j);
aea21c57 1368 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1369 ks = jparticle->GetStatusCode();
1370 km = jparticle->GetFirstMother();
1371
1372 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1373 nPartAcc++;
1374 }
0f6ee828 1375 if( numberOfAcceptedParticles <= nPartAcc){
1376 checking = kTRUE;
1377 break;
1378 }
aea21c57 1379 }
0f6ee828 1380
aea21c57 1381 return checking;
aea21c57 1382}
1383
5fa4b20b 1384void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1385{
1058d9df 1386 //
1387 // Load event into Pythia Common Block
1388 //
1389
1390 Int_t npart = stack -> GetNprimary();
1391 Int_t n0 = 0;
1392
1393 if (!flag) {
1394 (fPythia->GetPyjets())->N = npart;
1395 } else {
1396 n0 = (fPythia->GetPyjets())->N;
1397 (fPythia->GetPyjets())->N = n0 + npart;
1398 }
1399
1400
1401 for (Int_t part = 0; part < npart; part++) {
1402 TParticle *mPart = stack->Particle(part);
32d6ef7d 1403
1058d9df 1404 Int_t kf = mPart->GetPdgCode();
1405 Int_t ks = mPart->GetStatusCode();
1406 Int_t idf = mPart->GetFirstDaughter();
1407 Int_t idl = mPart->GetLastDaughter();
1408
1409 if (reHadr) {
1410 if (ks == 11 || ks == 12) {
1411 ks -= 10;
1412 idf = -1;
1413 idl = -1;
1414 }
32d6ef7d 1415 }
1416
1058d9df 1417 Float_t px = mPart->Px();
1418 Float_t py = mPart->Py();
1419 Float_t pz = mPart->Pz();
1420 Float_t e = mPart->Energy();
1421 Float_t m = mPart->GetCalcMass();
32d6ef7d 1422
1058d9df 1423
1424 (fPythia->GetPyjets())->P[0][part+n0] = px;
1425 (fPythia->GetPyjets())->P[1][part+n0] = py;
1426 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1427 (fPythia->GetPyjets())->P[3][part+n0] = e;
1428 (fPythia->GetPyjets())->P[4][part+n0] = m;
1429
1430 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1431 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1432 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1433 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1434 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1435 }
1436}
1437
c2fc9d31 1438void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1058d9df 1439{
1440 //
1441 // Load event into Pythia Common Block
1442 //
1443
1444 Int_t npart = stack -> GetEntries();
1445 Int_t n0 = 0;
1446
1447 if (!flag) {
1448 (fPythia->GetPyjets())->N = npart;
1449 } else {
1450 n0 = (fPythia->GetPyjets())->N;
1451 (fPythia->GetPyjets())->N = n0 + npart;
1452 }
1453
1454
1455 for (Int_t part = 0; part < npart; part++) {
1456 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
92847124 1457 if (!mPart) continue;
1458
1058d9df 1459 Int_t kf = mPart->GetPdgCode();
1460 Int_t ks = mPart->GetStatusCode();
1461 Int_t idf = mPart->GetFirstDaughter();
1462 Int_t idl = mPart->GetLastDaughter();
1463
1464 if (reHadr) {
92847124 1465 if (ks == 11 || ks == 12) {
1466 ks -= 10;
1467 idf = -1;
1468 idl = -1;
1469 }
8d2cd130 1470 }
1058d9df 1471
1472 Float_t px = mPart->Px();
1473 Float_t py = mPart->Py();
1474 Float_t pz = mPart->Pz();
1475 Float_t e = mPart->Energy();
1476 Float_t m = mPart->GetCalcMass();
1477
1478
1479 (fPythia->GetPyjets())->P[0][part+n0] = px;
1480 (fPythia->GetPyjets())->P[1][part+n0] = py;
1481 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1482 (fPythia->GetPyjets())->P[3][part+n0] = e;
1483 (fPythia->GetPyjets())->P[4][part+n0] = m;
1484
1485 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1486 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1487 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1488 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1489 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1490 }
8d2cd130 1491}
1492
5fa4b20b 1493
014a9521 1494void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1495{
1496//
1497// Calls the Pythia jet finding algorithm to find jets in the current event
1498//
1499//
8d2cd130 1500//
1501// Save jets
1502 Int_t n = fPythia->GetN();
1503
1504//
1505// Run Jet Finder
1506 fPythia->Pycell(njets);
1507 Int_t i;
1508 for (i = 0; i < njets; i++) {
1509 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1510 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1511 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1512 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1513
1514 jets[0][i] = px;
1515 jets[1][i] = py;
1516 jets[2][i] = pz;
1517 jets[3][i] = e;
1518 }
1519}
1520
1521
1522
1523void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1524{
1525//
1526// Calls the Pythia clustering algorithm to find jets in the current event
1527//
1528 Int_t n = fPythia->GetN();
1529 nJets = 0;
1530 nJetsTrig = 0;
1531 if (fJetReconstruction == kCluster) {
1532//
1533// Configure cluster algorithm
1534//
1535 fPythia->SetPARU(43, 2.);
1536 fPythia->SetMSTU(41, 1);
1537//
1538// Call cluster algorithm
1539//
1540 fPythia->Pyclus(nJets);
1541//
1542// Loading jets from common block
1543//
1544 } else {
592f8307 1545
8d2cd130 1546//
1547// Run Jet Finder
1548 fPythia->Pycell(nJets);
1549 }
1550
1551 Int_t i;
1552 for (i = 0; i < nJets; i++) {
1553 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1554 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1555 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1556 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1557 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1558 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1559 Float_t theta = TMath::ATan2(pt,pz);
1560 Float_t et = e * TMath::Sin(theta);
1561 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1562 if (
1563 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1564 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1565 et > fEtMinJet && et < fEtMaxJet
1566 )
1567 {
1568 jets[0][nJetsTrig] = px;
1569 jets[1][nJetsTrig] = py;
1570 jets[2][nJetsTrig] = pz;
1571 jets[3][nJetsTrig] = e;
1572 nJetsTrig++;
5fa4b20b 1573// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1574 } else {
1575// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1576 }
1577 }
1578}
1579
f913ec4f 1580void AliGenPythia::GetSubEventTime()
1581{
1582 // Calculates time of the next subevent
9d7108a7 1583 fEventTime = 0.;
1584 if (fEventsTime) {
1585 TArrayF &array = *fEventsTime;
1586 fEventTime = array[fCurSubEvent++];
1587 }
1588 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1589 return;
f913ec4f 1590}
8d2cd130 1591
e81f2bac 1592Bool_t AliGenPythia::IsInBarrel(Float_t eta) const
40fe59d4 1593{
1594 // Is particle in Central Barrel acceptance?
1595 // etamin=-etamax
1596 if( eta < fTriggerEta )
1597 return kTRUE;
1598 else
1599 return kFALSE;
1600}
1601
e81f2bac 1602Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
ec2c406e 1603{
1604 // Is particle in EMCAL acceptance?
ec2c406e 1605 // phi in degrees, etamin=-etamax
9fd8e520 1606 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1607 eta < fEMCALEta )
ec2c406e 1608 return kTRUE;
1609 else
1610 return kFALSE;
1611}
1612
e81f2bac 1613Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
ec2c406e 1614{
1615 // Is particle in PHOS acceptance?
1616 // Acceptance slightly larger considered.
1617 // phi in degrees, etamin=-etamax
40fe59d4 1618 // iparticle is the index of the particle to be checked, for PHOS rotation case
1619
9fd8e520 1620 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1621 eta < fPHOSEta )
ec2c406e 1622 return kTRUE;
1623 else
40fe59d4 1624 {
1625 if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1626
ec2c406e 1627 return kFALSE;
40fe59d4 1628 }
ec2c406e 1629}
1630
40fe59d4 1631void AliGenPythia::RotatePhi(Bool_t& okdd)
90a236ce 1632{
40fe59d4 1633 //Rotate event in phi to enhance events in PHOS acceptance
1634
1635 if(fPHOSRotateCandidate < 0) return ;
1636
90a236ce 1637 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1638 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1639 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
800be8b7 1640 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
90a236ce 1641
1642 //calculate deltaphi
40fe59d4 1643 TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
90a236ce 1644 Double_t phphi = ph->Phi();
1645 Double_t deltaphi = phiPHOS - phphi;
40fe59d4 1646
90a236ce 1647
1648
1649 //loop for all particles and produce the phi rotation
8507138f 1650 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
90a236ce 1651 Double_t oldphi, newphi;
777e69b0 1652 Double_t newVx, newVy, r, vZ, time;
1653 Double_t newPx, newPy, pt, pz, e;
90a236ce 1654 for(Int_t i=0; i< np; i++) {
40fe59d4 1655 TParticle* iparticle = (TParticle *) fParticles.At(i);
1656 oldphi = iparticle->Phi();
1657 newphi = oldphi + deltaphi;
1658 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1659 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1660
1661 r = iparticle->R();
1662 newVx = r * TMath::Cos(newphi);
1663 newVy = r * TMath::Sin(newphi);
1664 vZ = iparticle->Vz(); // don't transform
1665 time = iparticle->T(); // don't transform
1666
1667 pt = iparticle->Pt();
1668 newPx = pt * TMath::Cos(newphi);
1669 newPy = pt * TMath::Sin(newphi);
1670 pz = iparticle->Pz(); // don't transform
1671 e = iparticle->Energy(); // don't transform
1672
1673 // apply rotation
1674 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1675 iparticle->SetMomentum(newPx, newPy, pz, e);
1676
90a236ce 1677 } //end particle loop
1678
40fe59d4 1679 // now let's check that we put correctly the candidate photon in PHOS
1680 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1681 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1682 if(IsInPHOS(phi,eta,-1))
1683 okdd = kTRUE;
1684
1685 // reset the value for next event
1686 fPHOSRotateCandidate = -1;
1687
90a236ce 1688}
ec2c406e 1689
1690
800be8b7 1691Bool_t AliGenPythia::CheckDiffraction()
1692{
1693 // use this method only with Perugia-0 tune!
1694
1695 // printf("AAA\n");
1696
1697 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1698
1699 Int_t iPart1=-1;
1700 Int_t iPart2=-1;
1701
1702 Double_t y1 = 1e10;
1703 Double_t y2 = -1e10;
1704
1705 const Int_t kNstable=20;
1706 const Int_t pdgStable[20] = {
1707 22, // Photon
1708 11, // Electron
1709 12, // Electron Neutrino
1710 13, // Muon
1711 14, // Muon Neutrino
1712 15, // Tau
1713 16, // Tau Neutrino
1714 211, // Pion
1715 321, // Kaon
1716 311, // K0
1717 130, // K0s
1718 310, // K0l
1719 2212, // Proton
1720 2112, // Neutron
1721 3122, // Lambda_0
1722 3112, // Sigma Minus
1723 3222, // Sigma Plus
1724 3312, // Xsi Minus
1725 3322, // Xsi0
1726 3334 // Omega
1727 };
1728
1729 for (Int_t i = 0; i < np; i++) {
1730 TParticle * part = (TParticle *) fParticles.At(i);
1731
1732 Int_t statusCode = part->GetStatusCode();
1733
1734 // Initial state particle
1735 if (statusCode != 1)
1736 continue;
1737
1738 Int_t pdg = TMath::Abs(part->GetPdgCode());
1739 Bool_t isStable = kFALSE;
1740 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1741 if (pdg == pdgStable[i1]) {
1742 isStable = kTRUE;
1743 break;
1744 }
1745 }
1746 if(!isStable)
1747 continue;
1748
1749 Double_t y = part->Y();
1750
1751 if (y < y1)
1752 {
1753 y1 = y;
1754 iPart1 = i;
1755 }
1756 if (y > y2)
1757 {
1758 y2 = y;
1759 iPart2 = i;
1760 }
1761 }
1762
1763 if(iPart1<0 || iPart2<0) return kFALSE;
1764
1765 y1=TMath::Abs(y1);
1766 y2=TMath::Abs(y2);
1767
1768 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1769 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1770
1771 Int_t pdg1 = part1->GetPdgCode();
1772 Int_t pdg2 = part2->GetPdgCode();
1773
1774
1775 Int_t iPart = -1;
1776 if (pdg1 == 2212 && pdg2 == 2212)
1777 {
1778 if(y1 > y2)
1779 iPart = iPart1;
1780 else if(y1 < y2)
1781 iPart = iPart2;
1782 else {
1783 iPart = iPart1;
1784 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1785 }
1786 }
1787 else if (pdg1 == 2212)
1788 iPart = iPart1;
1789 else if (pdg2 == 2212)
1790 iPart = iPart2;
1791
1792
1793
1794
1795
1796 Double_t M=-1.;
1797 if(iPart>0) {
1798 TParticle * part = (TParticle *) fParticles.At(iPart);
1799 Double_t E= part->Energy();
1800 Double_t P= part->P();
4f813e90 1801 Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
1802 if(M2<0) return kFALSE;
1803 M= TMath::Sqrt(M2);
800be8b7 1804 }
1805
c5dfa3e4 1806 Double_t Mmin, Mmax, wSD, wDD, wND;
1807 if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1808
1809 if(M>-1 && M<Mmin) return kFALSE;
1810 if(M>Mmax) M=-1;
1811
1812 Int_t procType=fPythia->GetMSTI(1);
1813 Int_t proc0=2;
1814 if(procType== 94) proc0=1;
1815 if(procType== 92 || procType== 93) proc0=0;
1816
1817 Int_t proc=2;
1818 if(M>0) proc=0;
1819 else if(proc0==1) proc=1;
1820
1821 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1822 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
800be8b7 1823
1824
1825 // if(proc==1 || proc==2) return kFALSE;
1826
c5dfa3e4 1827 if(proc!=0) {
1828 if(proc0!=0) fProcDiff = procType;
1829 else fProcDiff = 95;
1830 return kTRUE;
1831 }
1832
1833 if(wSD<0) AliError("wSD<0 ! \n");
1834
1835 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1836
1837 // printf("iPart = %d\n", iPart);
1838
1839 if(iPart==iPart1) fProcDiff=93;
1840 else if(iPart==iPart2) fProcDiff=92;
1841 else {
1842 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1843
800be8b7 1844 }
1845
c5dfa3e4 1846 return kTRUE;
1847}
1848
1849
1850
1851Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1852 Double_t &wSD, Double_t &wDD, Double_t &wND)
1853{
1854
1855 // 900 GeV
1856 if(TMath::Abs(fEnergyCMS-900)<1 ){
1857
1858const Int_t nbin=400;
1859Double_t bin[]={
18601.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
18614.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
18627.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
186310.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
186413.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
186515.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
186618.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
186721.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
186824.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
186927.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
187030.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
187133.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
187236.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
187339.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
187442.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
187545.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
187648.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
187751.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
187854.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
187957.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
188060.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
188163.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
188266.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
188369.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
188472.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
188575.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
188678.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
188781.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
188884.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
188987.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
189090.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
189193.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
189296.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
189399.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1894102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1895105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1896108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1897111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1898114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1899117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1900120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1901123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1902126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
1903129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
1904132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
1905135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
1906138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
1907141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
1908144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
1909147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
1910150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
1911153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
1912156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
1913159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
1914162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
1915165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
1916168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
1917171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
1918174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
1919177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
1920180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
1921183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
1922186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
1923189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
1924192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
1925195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
1926198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1927Double_t w[]={
19281.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
19290.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
19300.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
19310.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
19320.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
19330.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
19340.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
19350.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
19360.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
19370.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
19380.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
19390.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
19400.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
19410.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
19420.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
19430.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
19440.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
19450.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
19460.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
19470.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
19480.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
19490.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
19500.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
19510.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
19520.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
19530.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
19540.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
19550.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
19560.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
19570.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
19580.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
19590.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
19600.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
19610.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
19620.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
19630.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
19640.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
19650.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
19660.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
19670.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
19680.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
19690.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
19700.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
19710.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
19720.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
19730.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
19740.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
19750.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
19760.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
19770.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
19780.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
19790.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
19800.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
19810.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
19820.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
19830.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
19840.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
19850.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
19860.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
19870.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
19880.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
19890.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
19900.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
19910.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
19920.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
19930.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
19940.386112, 0.364314, 0.434375, 0.334629};
1995wDD = 0.379611;
1996wND = 0.496961;
1997wSD = -1;
1998
1999 Mmin = bin[0];
2000 Mmax = bin[nbin];
2001 if(M<Mmin || M>Mmax) return kTRUE;
2002
7873275c 2003 Int_t ibin=nbin-1;
93a60586 2004 for(Int_t i=1; i<=nbin; i++)
2ff57420 2005 if(M<=bin[i]) {
2006 ibin=i-1;
2007 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
800be8b7 2008 break;
2009 }
c5dfa3e4 2010 wSD=w[ibin];
2011 return kTRUE;
2012 }
2013 else if(TMath::Abs(fEnergyCMS-2760)<1 ){
800be8b7 2014
c5dfa3e4 2015const Int_t nbin=400;
2016Double_t bin[]={
20171.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
20184.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
20197.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
202010.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
202113.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
202215.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
202318.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
202421.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
202524.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
202627.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
202730.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
202833.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
202936.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
203039.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
203142.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
203245.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
203348.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
203451.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
203554.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
203657.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
203760.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
203863.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
203966.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
204069.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
204172.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
204275.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
204378.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
204481.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
204584.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
204687.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
204790.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
204893.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
204996.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
205099.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2051102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2052105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2053108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2054111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2055114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2056117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2057120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2058123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2059126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2060129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2061132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2062135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2063138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2064141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2065144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2066147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2067150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2068153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2069156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2070159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2071162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2072165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2073168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2074171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2075174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2076177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2077180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2078183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2079186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2080189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2081192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2082195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2083198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2084Double_t w[]={
20851.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
20860.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
20870.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
20880.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
20890.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
20900.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
20910.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
20920.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
20930.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
20940.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
20950.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
20960.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
20970.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
20980.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
20990.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
21000.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
21010.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
21020.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
21030.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
21040.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
21050.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
21060.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
21070.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
21080.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
21090.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
21100.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
21110.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
21120.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
21130.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
21140.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
21150.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
21160.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
21170.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
21180.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
21190.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
21200.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
21210.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
21220.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
21230.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
21240.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
21250.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
21260.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
21270.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
21280.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
21290.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
21300.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
21310.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
21320.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
21330.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
21340.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
21350.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
21360.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
21370.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
21380.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
21390.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
21400.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
21410.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
21420.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
21430.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
21440.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
21450.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
21460.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
21470.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
21480.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
21490.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
21500.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
21510.201712, 0.242225, 0.255565, 0.258738};
2152wDD = 0.512813;
2153wND = 0.518820;
2154wSD = -1;
2155
2156 Mmin = bin[0];
2157 Mmax = bin[nbin];
2158 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2159
c5dfa3e4 2160 Int_t ibin=nbin-1;
2161 for(Int_t i=1; i<=nbin; i++)
2162 if(M<=bin[i]) {
2163 ibin=i-1;
2164 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2165 break;
2166 }
2167 wSD=w[ibin];
2168 return kTRUE;
2169 }
800be8b7 2170
800be8b7 2171
c5dfa3e4 2172 else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2173const Int_t nbin=400;
2174Double_t bin[]={
21751.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
21764.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
21777.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
217810.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
217913.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
218015.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
218118.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
218221.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
218324.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
218427.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
218530.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
218633.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
218736.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
218839.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
218942.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
219045.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
219148.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
219251.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
219354.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
219457.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
219560.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
219663.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
219766.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
219869.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
219972.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
220075.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
220178.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
220281.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
220384.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
220487.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
220590.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
220693.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
220796.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
220899.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2209102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2210105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2211108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2212111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2213114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2214117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2215120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2216123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2217126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2218129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2219132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2220135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2221138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2222141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2223144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2224147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2225150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2226153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2227156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2228159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2229162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2230165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2231168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2232171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2233174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2234177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2235180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2236183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2237186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2238189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2239192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2240195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2241198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2242Double_t w[]={
22431.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
22440.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
22450.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
22460.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
22470.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
22480.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
22490.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
22500.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
22510.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
22520.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
22530.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
22540.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
22550.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
22560.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
22570.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
22580.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
22590.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
22600.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
22610.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
22620.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
22630.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
22640.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
22650.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
22660.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
22670.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
22680.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
22690.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
22700.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
22710.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
22720.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
22730.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
22740.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
22750.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
22760.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
22770.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
22780.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
22790.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
22800.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
22810.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
22820.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
22830.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
22840.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
22850.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
22860.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
22870.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
22880.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
22890.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
22900.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
22910.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
22920.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
22930.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
22940.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
22950.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
22960.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
22970.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
22980.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
22990.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
23000.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
23010.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
23020.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
23030.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
23040.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
23050.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
23060.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
23070.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
23080.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
23090.175006, 0.223482, 0.202706, 0.218108};
2310wDD = 0.207705;
2311wND = 0.289628;
2312wSD = -1;
2313
2314 Mmin = bin[0];
2315 Mmax = bin[nbin];
2316
2317 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2318
c5dfa3e4 2319 Int_t ibin=nbin-1;
2320 for(Int_t i=1; i<=nbin; i++)
2321 if(M<=bin[i]) {
2322 ibin=i-1;
2323 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2324 break;
2325 }
2326 wSD=w[ibin];
800be8b7 2327 return kTRUE;
c5dfa3e4 2328 }
2329
2330 return kFALSE;
800be8b7 2331}
06956108 2332
2333Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2334{
2335// Check if this is a heavy flavor decay product
2336 TParticle * part = (TParticle *) fParticles.At(ipart);
2337 Int_t mpdg = TMath::Abs(part->GetPdgCode());
2338 Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2339 //
2340 // Light hadron
2341 if (mfl >= 4 && mfl < 6) return kTRUE;
2342
2343 Int_t imo = part->GetFirstMother()-1;
2344 TParticle* pm = part;
2345 //
2346 // Heavy flavor hadron produced by generator
2347 while (imo > -1) {
2348 pm = (TParticle*)fParticles.At(imo);
2349 mpdg = TMath::Abs(pm->GetPdgCode());
2350 mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2351 if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2352 imo = pm->GetFirstMother()-1;
2353 }
2354 return kFALSE;
2355}
40fe59d4 2356
e81f2bac 2357Bool_t AliGenPythia::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
40fe59d4 2358{
2359 // check the eta/phi correspond to the detectors acceptance
2360 // iparticle is the index of the particle to be checked, for PHOS rotation case
2361 if (fCheckPHOS && IsInPHOS (phi,eta,iparticle)) return kTRUE;
2362 else if(fCheckEMCAL && IsInEMCAL (phi,eta)) return kTRUE;
2363 else if(fCheckBarrel && IsInBarrel( eta)) return kTRUE;
2364 else return kFALSE;
2365}
2366
e81f2bac 2367Bool_t AliGenPythia::TriggerOnSelectedParticles(Int_t np)
40fe59d4 2368{
2369 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2370 // implemented primaryly for kPyJets, but extended to any kind of process.
2371
2372 //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2373 // fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
2374
2375 Bool_t ok = kFALSE;
2376 for (Int_t i=0; i< np; i++) {
2377
2378 TParticle* iparticle = (TParticle *) fParticles.At(i);
2379
2380 Int_t pdg = iparticle->GetPdgCode();
2381 Int_t status = iparticle->GetStatusCode();
2382 Int_t imother = iparticle->GetFirstMother() - 1;
2383
2384 TParticle* pmother = 0x0;
2385 Int_t momStatus = -1;
2386 Int_t momPdg = -1;
2387 if(imother > 0 ){
2388 pmother = (TParticle *) fParticles.At(imother);
2389 momStatus = pmother->GetStatusCode();
2390 momPdg = pmother->GetPdgCode();
2391 }
2392
2393 ok = kFALSE;
2394
2395 //
2396 // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2397 //
2398 // Hadron
2399 if (fHadronInCalo && status == 1)
2400 {
2401 if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
2402 // (in case neutral mesons were declared stable)
2403 ok = kTRUE;
2404 }
2405 //Electron
2406 else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2407 {
2408 ok = kTRUE;
2409 }
2410 //Fragmentation photon
2411 else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2412 {
2413 if(momStatus != 11) ok = kTRUE ; // No photon from hadron decay
2414 }
2415 // Decay photon
2416 else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2417 {
2418 if( momStatus == 11)
2419 {
2420 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2421 // pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2422 ok = kTRUE ; // photon from hadron decay
2423
2424 //In case only decays from pi0 or eta requested
2425 if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2426 if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2427 }
2428
2429 }
2430 // Pi0 or Eta particle
2431 else if ((fPi0InCalo || fEtaInCalo))
2432 {
2433 if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2434
2435 if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2436 {
2437 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2438 ok = kTRUE;
2439 }
2440 else if (fEtaInCalo && pdg == 221)
2441 {
2442 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2443 ok = kTRUE;
2444 }
2445
2446 }// pi0 or eta
2447
2448 //
2449 // Check that the selected particle is in the calorimeter acceptance
2450 //
2451 if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2452 {
2453 //Just check if the selected particle falls in the acceptance
2454 if(!fForceNeutralMeson2PhotonDecay )
2455 {
2456 //printf("\t Check acceptance! \n");
2457 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2458 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2459
2460 if(CheckDetectorAcceptance(phi,eta,i))
2461 {
2462 ok =kTRUE;
2463 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));
2464 //printf("\t Accept \n");
2465 break;
2466 }
2467 else ok = kFALSE;
2468 }
2469 //Mesons have several decay modes, select only those decaying into 2 photons
2470 else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2471 {
2472 // In case we want the pi0/eta trigger,
2473 // check the decay mode (2 photons)
2474
2475 //printf("\t Force decay 2 gamma\n");
2476
2477 Int_t ndaughters = iparticle->GetNDaughters();
2478 if(ndaughters != 2){
2479 ok=kFALSE;
2480 continue;
2481 }
2482
2483 TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2484 TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2485 if(!d1 || !d2) {
2486 ok=kFALSE;
2487 continue;
2488 }
2489
2490 //iparticle->Print();
2491 //d1->Print();
2492 //d2->Print();
2493
2494 Int_t pdgD1 = d1->GetPdgCode();
2495 Int_t pdgD2 = d2->GetPdgCode();
2496 //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2497 //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2498
2499 if(pdgD1 != 22 || pdgD2 != 22){
2500 ok = kFALSE;
2501 continue;
2502 }
2503
2504 //printf("\t accept decay\n");
2505
2506 //Trigger on the meson, not on the daughter
2507 if(!fDecayPhotonInCalo){
2508
2509 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2510 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2511
2512 if(CheckDetectorAcceptance(phi,eta,i))
2513 {
2514 //printf("\t Accept meson pdg %d\n",pdg);
2515 ok =kTRUE;
2516 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));
2517 break;
2518 } else {
2519 ok=kFALSE;
2520 continue;
2521 }
2522 }
2523
2524 //printf("Check daughters acceptance\n");
2525
2526 //Trigger on the meson daughters
2527 //Photon 1
2528 Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2529 Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
2530 if(d1->Pt() > fTriggerParticleMinPt)
2531 {
2532 //printf("\t Check acceptance photon 1! \n");
2533 if(CheckDetectorAcceptance(phi,eta,i))
2534 {
2535 //printf("\t Accept Photon 1\n");
2536 ok =kTRUE;
2537 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));
2538 break;
2539 }
2540 else ok = kFALSE;
2541 } // pt cut
2542 else ok = kFALSE;
2543
2544 //Photon 2
2545 phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2546 eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2547
2548 if(d2->Pt() > fTriggerParticleMinPt)
2549 {
2550 //printf("\t Check acceptance photon 2! \n");
2551 if(CheckDetectorAcceptance(phi,eta,i))
2552 {
2553 //printf("\t Accept Photon 2\n");
2554 ok =kTRUE;
2555 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));
2556 break;
2557 }
2558 else ok = kFALSE;
2559 } // pt cut
2560 else ok = kFALSE;
2561 } // force 2 photon daughters in pi0/eta decays
2562 else ok = kFALSE;
2563 } else ok = kFALSE; // check acceptance
2564 } // primary loop
2565
2566 //
2567 // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2568 // A particle passing all trigger conditions except phi position in PHOS, is used as reference
2569 //
2570 if(fCheckPHOSeta)
2571 {
2572 RotatePhi(ok);
2573 }
2574
2575 return ok;
2576}
2577