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