]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA6/AliGenPythia.cxx
Adding analysis task
[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:
589380c6 514 break;
8d2cd130 515 }
516//
592f8307 517//
518// JetFinder for Trigger
519//
520// Configure detector (EMCAL like)
521//
d7de4a67 522 fPythia->SetPARU(51, fPycellEtaMax);
523 fPythia->SetMSTU(51, fPycellNEta);
524 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 525//
526// Configure Jet Finder
527//
d7de4a67 528 fPythia->SetPARU(58, fPycellThreshold);
529 fPythia->SetPARU(52, fPycellEtSeed);
530 fPythia->SetPARU(53, fPycellMinEtJet);
531 fPythia->SetPARU(54, fPycellMaxRadius);
532 fPythia->SetMSTU(54, 2);
592f8307 533//
8d2cd130 534// This counts the total number of calls to Pyevnt() per run.
535 fTrialsRun = 0;
536 fQ = 0.;
537 fX1 = 0.;
538 fX2 = 0.;
539 fNev = 0 ;
540//
1d568bc2 541//
542//
8d2cd130 543 AliGenMC::Init();
1d568bc2 544//
545//
546//
547 if (fSetNuclei) {
548 fDyBoost = 0;
549 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
550 }
d7de4a67 551
cd07c39b 552 fPythia->SetPARJ(200, 0.0);
553 fPythia->SetPARJ(199, 0.0);
554 fPythia->SetPARJ(198, 0.0);
555 fPythia->SetPARJ(197, 0.0);
556
557 if (fQuench == 1) {
5fa4b20b 558 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 559 }
3a709cfa 560
66b8652c 561 if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
562
b976f7d8 563 if (fQuench == 3) {
564 // Nestor's change of the splittings
565 fPythia->SetPARJ(200, 0.8);
566 fPythia->SetMSTJ(41, 1); // QCD radiation only
567 fPythia->SetMSTJ(42, 2); // angular ordering
568 fPythia->SetMSTJ(44, 2); // option to run alpha_s
569 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
570 fPythia->SetMSTJ(50, 0); // No coherence in first branching
571 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
cd07c39b 572 } else if (fQuench == 4) {
573 // Armesto-Cunqueiro-Salgado change of the splittings.
574 AliFastGlauber* glauber = AliFastGlauber::Instance();
575 glauber->Init(2);
576 //read and store transverse almonds corresponding to differnt
577 //impact parameters.
cd07c39b 578 glauber->SetCentralityClass(0.,0.1);
579 fPythia->SetPARJ(200, 1.);
580 fPythia->SetPARJ(198, fQhat);
581 fPythia->SetPARJ(199, fLength);
cd07c39b 582 fPythia->SetMSTJ(42, 2); // angular ordering
583 fPythia->SetMSTJ(44, 2); // option to run alpha_s
cd07c39b 584 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
b976f7d8 585 }
8d2cd130 586}
587
588void AliGenPythia::Generate()
589{
590// Generate one event
19ef8cf0 591 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
8d2cd130 592 fDecayer->ForceDecay();
593
13cca24a 594 Double_t polar[3] = {0,0,0};
595 Double_t origin[3] = {0,0,0};
596 Double_t p[4];
8d2cd130 597// converts from mm/c to s
598 const Float_t kconv=0.001/2.999792458e8;
599//
600 Int_t nt=0;
601 Int_t jev=0;
602 Int_t j, kf;
603 fTrials=0;
f913ec4f 604 fEventTime = 0.;
605
2590ccf9 606
8d2cd130 607
608 // Set collision vertex position
2590ccf9 609 if (fVertexSmear == kPerEvent) Vertex();
610
8d2cd130 611// event loop
612 while(1)
613 {
32d6ef7d 614//
5fa4b20b 615// Produce event
32d6ef7d 616//
32d6ef7d 617//
5fa4b20b 618// Switch hadronisation off
619//
620 fPythia->SetMSTJ(1, 0);
cd07c39b 621
622 if (fQuench ==4){
623 Double_t bimp;
624 // Quenching comes through medium-modified splitting functions.
625 AliFastGlauber::Instance()->GetRandomBHard(bimp);
e6fe9b82 626 fPythia->SetPARJ(197, bimp);
627 fImpact = bimp;
6c43eccb 628 fPythia->Qpygin0();
cd07c39b 629 }
32d6ef7d 630//
5fa4b20b 631// Either produce new event or read partons from file
632//
633 if (!fReadFromFile) {
beac474c 634 if (!fNewMIS) {
635 fPythia->Pyevnt();
636 } else {
637 fPythia->Pyevnw();
638 }
5fa4b20b 639 fNpartons = fPythia->GetN();
640 } else {
33c3c91a 641 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
642 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
5fa4b20b 643 fPythia->SetN(0);
644 LoadEvent(fRL->Stack(), 0 , 1);
645 fPythia->Pyedit(21);
646 }
647
32d6ef7d 648//
649// Run quenching routine
650//
5fa4b20b 651 if (fQuench == 1) {
652 fPythia->Quench();
653 } else if (fQuench == 2){
654 fPythia->Pyquen(208., 0, 0.);
b976f7d8 655 } else if (fQuench == 3) {
656 // Quenching is via multiplicative correction of the splittings
5fa4b20b 657 }
b976f7d8 658
32d6ef7d 659//
5fa4b20b 660// Switch hadronisation on
32d6ef7d 661//
20e47f08 662 if (fHadronisation) {
663 fPythia->SetMSTJ(1, 1);
5fa4b20b 664//
665// .. and perform hadronisation
aea21c57 666// printf("Calling hadronisation %d\n", fPythia->GetN());
03358a32 667
668 if (fPatchOmegaDalitz) {
669 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
670 fPythia->Pyexec();
671 fPythia->DalitzDecays();
672 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
673 }
674 fPythia->Pyexec();
20e47f08 675 }
8d2cd130 676 fTrials++;
8507138f 677 fPythia->ImportParticles(&fParticles,"All");
03358a32 678
df275cfa 679 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
e1cf8975 680 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
f64640b5 681
8d2cd130 682//
683//
684//
685 Int_t i;
686
dbd64db6 687 fNprimaries = 0;
8507138f 688 Int_t np = fParticles.GetEntriesFast();
5fa4b20b 689
7c21f297 690 if (np == 0) continue;
8d2cd130 691//
2590ccf9 692
8d2cd130 693//
694 Int_t* pParent = new Int_t[np];
695 Int_t* pSelected = new Int_t[np];
696 Int_t* trackIt = new Int_t[np];
5fa4b20b 697 for (i = 0; i < np; i++) {
8d2cd130 698 pParent[i] = -1;
699 pSelected[i] = 0;
700 trackIt[i] = 0;
701 }
702
703 Int_t nc = 0; // Total n. of selected particles
704 Int_t nParents = 0; // Selected parents
705 Int_t nTkbles = 0; // Trackable particles
f529e69b 706 if (fProcess != kPyMbDefault &&
707 fProcess != kPyMb &&
6d2ec66d 708 fProcess != kPyMbAtlasTuneMC09 &&
04081a91 709 fProcess != kPyMbWithDirectPhoton &&
f529e69b 710 fProcess != kPyJets &&
8d2cd130 711 fProcess != kPyDirectGamma &&
589380c6 712 fProcess != kPyMbNonDiffr &&
d7de4a67 713 fProcess != kPyMbMSEL1 &&
f529e69b 714 fProcess != kPyW &&
715 fProcess != kPyZ &&
716 fProcess != kPyCharmppMNRwmi &&
9dfe63b3 717 fProcess != kPyBeautyppMNRwmi &&
718 fProcess != kPyBeautyJets) {
8d2cd130 719
5fa4b20b 720 for (i = 0; i < np; i++) {
8507138f 721 TParticle* iparticle = (TParticle *) fParticles.At(i);
8d2cd130 722 Int_t ks = iparticle->GetStatusCode();
723 kf = CheckPDGCode(iparticle->GetPdgCode());
724// No initial state partons
725 if (ks==21) continue;
726//
727// Heavy Flavor Selection
728//
729 // quark ?
730 kf = TMath::Abs(kf);
731 Int_t kfl = kf;
9ff6c04c 732 // Resonance
f913ec4f 733
9ff6c04c 734 if (kfl > 100000) kfl %= 100000;
183a5ca9 735 if (kfl > 10000) kfl %= 10000;
8d2cd130 736 // meson ?
737 if (kfl > 10) kfl/=100;
738 // baryon
739 if (kfl > 10) kfl/=10;
8d2cd130 740 Int_t ipa = iparticle->GetFirstMother()-1;
741 Int_t kfMo = 0;
f913ec4f 742//
743// Establish mother daughter relation between heavy quarks and mesons
744//
745 if (kf >= fFlavorSelect && kf <= 6) {
746 Int_t idau = iparticle->GetFirstDaughter() - 1;
747 if (idau > -1) {
8507138f 748 TParticle* daughter = (TParticle *) fParticles.At(idau);
f913ec4f 749 Int_t pdgD = daughter->GetPdgCode();
750 if (pdgD == 91 || pdgD == 92) {
751 Int_t jmin = daughter->GetFirstDaughter() - 1;
752 Int_t jmax = daughter->GetLastDaughter() - 1;
2ab330c9 753 for (Int_t jp = jmin; jp <= jmax; jp++)
754 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
f913ec4f 755 } // is string or cluster
756 } // has daughter
757 } // heavy quark
8d2cd130 758
f913ec4f 759
8d2cd130 760 if (ipa > -1) {
8507138f 761 TParticle * mother = (TParticle *) fParticles.At(ipa);
8d2cd130 762 kfMo = TMath::Abs(mother->GetPdgCode());
763 }
f913ec4f 764
8d2cd130 765 // What to keep in Stack?
766 Bool_t flavorOK = kFALSE;
767 Bool_t selectOK = kFALSE;
768 if (fFeedDownOpt) {
32d6ef7d 769 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 770 } else {
32d6ef7d 771 if (kfl > fFlavorSelect) {
772 nc = -1;
773 break;
774 }
775 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 776 }
777 switch (fStackFillOpt) {
06956108 778 case kHeavyFlavor:
8d2cd130 779 case kFlavorSelection:
32d6ef7d 780 selectOK = kTRUE;
781 break;
8d2cd130 782 case kParentSelection:
32d6ef7d 783 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
784 break;
8d2cd130 785 }
786 if (flavorOK && selectOK) {
787//
788// Heavy flavor hadron or quark
789//
790// Kinematic seletion on final state heavy flavor mesons
791 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
792 {
9ff6c04c 793 continue;
8d2cd130 794 }
795 pSelected[i] = 1;
796 if (ParentSelected(kf)) ++nParents; // Update parent count
797// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
798 } else {
799// Kinematic seletion on decay products
800 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 801 && !KinematicSelection(iparticle, 1))
8d2cd130 802 {
9ff6c04c 803 continue;
8d2cd130 804 }
805//
806// Decay products
807// Select if mother was selected and is not tracked
808
809 if (pSelected[ipa] &&
810 !trackIt[ipa] && // mother will be tracked ?
811 kfMo != 5 && // mother is b-quark, don't store fragments
812 kfMo != 4 && // mother is c-quark, don't store fragments
813 kf != 92) // don't store string
814 {
815//
816// Semi-stable or de-selected: diselect decay products:
817//
818//
819 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
820 {
821 Int_t ipF = iparticle->GetFirstDaughter();
822 Int_t ipL = iparticle->GetLastDaughter();
823 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
824 }
825// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
826 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
827 }
828 }
829 if (pSelected[i] == -1) pSelected[i] = 0;
830 if (!pSelected[i]) continue;
831 // Count quarks only if you did not include fragmentation
832 if (fFragmentation && kf <= 10) continue;
9ff6c04c 833
8d2cd130 834 nc++;
835// Decision on tracking
836 trackIt[i] = 0;
837//
838// Track final state particle
839 if (ks == 1) trackIt[i] = 1;
840// Track semi-stable particles
d25cfd65 841 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 842// Track particles selected by process if undecayed.
843 if (fForceDecay == kNoDecay) {
844 if (ParentSelected(kf)) trackIt[i] = 1;
845 } else {
846 if (ParentSelected(kf)) trackIt[i] = 0;
847 }
848 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
849//
850//
851
852 } // particle selection loop
853 if (nc > 0) {
854 for (i = 0; i<np; i++) {
855 if (!pSelected[i]) continue;
8507138f 856 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 857 kf = CheckPDGCode(iparticle->GetPdgCode());
858 Int_t ks = iparticle->GetStatusCode();
859 p[0] = iparticle->Px();
860 p[1] = iparticle->Py();
861 p[2] = iparticle->Pz();
a920faf9 862 p[3] = iparticle->Energy();
863
2590ccf9 864 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
865 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
866 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
867
21391258 868 Float_t tof = fTime + kconv*iparticle->T();
8d2cd130 869 Int_t ipa = iparticle->GetFirstMother()-1;
870 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 871
872 PushTrack(fTrackIt*trackIt[i], iparent, kf,
873 p[0], p[1], p[2], p[3],
874 origin[0], origin[1], origin[2], tof,
875 polar[0], polar[1], polar[2],
876 kPPrimary, nt, 1., ks);
8d2cd130 877 pParent[i] = nt;
dbd64db6 878 KeepTrack(nt);
879 fNprimaries++;
642f15cf 880 } // PushTrack loop
8d2cd130 881 }
882 } else {
883 nc = GenerateMB();
884 } // mb ?
f913ec4f 885
886 GetSubEventTime();
8d2cd130 887
234f6d32 888 delete[] pParent;
889 delete[] pSelected;
890 delete[] trackIt;
8d2cd130 891
892 if (nc > 0) {
893 switch (fCountMode) {
894 case kCountAll:
895 // printf(" Count all \n");
896 jev += nc;
897 break;
898 case kCountParents:
899 // printf(" Count parents \n");
900 jev += nParents;
901 break;
902 case kCountTrackables:
903 // printf(" Count trackable \n");
904 jev += nTkbles;
905 break;
906 }
907 if (jev >= fNpart || fNpart == -1) {
908 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 909
8d2cd130 910 fQ += fPythia->GetVINT(51);
911 fX1 += fPythia->GetVINT(41);
912 fX2 += fPythia->GetVINT(42);
913 fTrialsRun += fTrials;
914 fNev++;
915 MakeHeader();
916 break;
917 }
918 }
919 } // event loop
920 SetHighWaterMark(nt);
921// adjust weight due to kinematic selection
b88f5cea 922// AdjustWeights();
8d2cd130 923// get cross-section
924 fXsection=fPythia->GetPARI(1);
925}
926
927Int_t AliGenPythia::GenerateMB()
928{
929//
930// Min Bias selection and other global selections
931//
932 Int_t i, kf, nt, iparent;
933 Int_t nc = 0;
13cca24a 934 Double_t p[4];
935 Double_t polar[3] = {0,0,0};
936 Double_t origin[3] = {0,0,0};
8d2cd130 937// converts from mm/c to s
938 const Float_t kconv=0.001/2.999792458e8;
939
e0e89f40 940
8507138f 941 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
aea21c57 942
5fa4b20b 943
e0e89f40 944
8d2cd130 945 Int_t* pParent = new Int_t[np];
946 for (i=0; i< np; i++) pParent[i] = -1;
2f405d65 947 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8507138f 948 TParticle* jet1 = (TParticle *) fParticles.At(6);
949 TParticle* jet2 = (TParticle *) fParticles.At(7);
234f6d32 950 if (!CheckTrigger(jet1, jet2)) {
951 delete [] pParent;
952 return 0;
953 }
8d2cd130 954 }
e0e89f40 955
40fe59d4 956
957 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
ac01c7c6 958 // implemented primaryly for kPyJets, but extended to any kind of process.
40fe59d4 959 if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
960 (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
961 Bool_t ok = TriggerOnSelectedParticles(np);
43e3b80a 962
963 if(!ok) {
964 delete[] pParent;
965 return 0;
ec2c406e 966 }
43e3b80a 967 }
9dfe63b3 968
800be8b7 969 // Check for diffraction
970 if(fkTuneForDiff) {
c5dfa3e4 971 if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
800be8b7 972 if(!CheckDiffraction()) {
973 delete [] pParent;
974 return 0;
975 }
976 }
977 }
978
700b9416 979 // Check for minimum multiplicity
980 if (fTriggerMultiplicity > 0) {
981 Int_t multiplicity = 0;
982 for (i = 0; i < np; i++) {
8507138f 983 TParticle * iparticle = (TParticle *) fParticles.At(i);
700b9416 984
985 Int_t statusCode = iparticle->GetStatusCode();
986
987 // Initial state particle
e296848e 988 if (statusCode != 1)
700b9416 989 continue;
38112f3f 990 // eta cut
700b9416 991 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
992 continue;
38112f3f 993 // pt cut
994 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
995 continue;
996
700b9416 997 TParticlePDG* pdgPart = iparticle->GetPDG();
998 if (pdgPart && pdgPart->Charge() == 0)
999 continue;
1000
1001 ++multiplicity;
1002 }
1003
1004 if (multiplicity < fTriggerMultiplicity) {
1005 delete [] pParent;
1006 return 0;
1007 }
38112f3f 1008 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
700b9416 1009 }
90a236ce 1010
90a236ce 1011
7ce1321b 1012 if (fTriggerParticle) {
1013 Bool_t triggered = kFALSE;
1014 for (i = 0; i < np; i++) {
8507138f 1015 TParticle * iparticle = (TParticle *) fParticles.At(i);
7ce1321b 1016 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 1017 if (kf != fTriggerParticle) continue;
7ce1321b 1018 if (iparticle->Pt() == 0.) continue;
1019 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
2591bd0e 1020 if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
7ce1321b 1021 triggered = kTRUE;
1022 break;
1023 }
234f6d32 1024 if (!triggered) {
1025 delete [] pParent;
1026 return 0;
1027 }
7ce1321b 1028 }
e0e89f40 1029
1030
1031 // Check if there is a ccbar or bbbar pair with at least one of the two
1032 // in fYMin < y < fYMax
2f405d65 1033
9dfe63b3 1034 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
9ccc3d0e 1035 TParticle *partCheck;
1036 TParticle *mother;
e0e89f40 1037 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 1038 Bool_t theChild=kFALSE;
5d76923e 1039 Bool_t triggered=kFALSE;
9ccc3d0e 1040 Float_t y;
1041 Int_t pdg,mpdg,mpdgUpperFamily;
e0e89f40 1042 for(i=0; i<np; i++) {
9ccc3d0e 1043 partCheck = (TParticle*)fParticles.At(i);
1044 pdg = partCheck->GetPdgCode();
1045 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1046 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1047 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1048 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1049 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1050 }
5d76923e 1051 if(fTriggerParticle) {
1052 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1053 }
9ccc3d0e 1054 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1055 Int_t mi = partCheck->GetFirstMother() - 1;
1056 if(mi<0) continue;
1057 mother = (TParticle*)fParticles.At(mi);
1058 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 1059 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 1060 if ( ParentSelected(mpdg) ||
1061 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1062 if (KinematicSelection(partCheck,1)) {
1063 theChild=kTRUE;
1064 }
1065 }
1066 }
e0e89f40 1067 }
9ccc3d0e 1068 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
234f6d32 1069 delete[] pParent;
e0e89f40 1070 return 0;
1071 }
9ccc3d0e 1072 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1073 delete[] pParent;
1074 return 0;
1075 }
5d76923e 1076 if(fTriggerParticle && !triggered) { // particle requested is not produced
1077 delete[] pParent;
1078 return 0;
1079 }
9ccc3d0e 1080
e0e89f40 1081 }
aea21c57 1082
0f6ee828 1083 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
f529e69b 1084 if ( (fProcess == kPyW ||
1085 fProcess == kPyZ ||
1086 fProcess == kPyMbDefault ||
1087 fProcess == kPyMb ||
6d2ec66d 1088 fProcess == kPyMbAtlasTuneMC09 ||
04081a91 1089 fProcess == kPyMbWithDirectPhoton ||
f529e69b 1090 fProcess == kPyMbNonDiffr)
0f6ee828 1091 && (fCutOnChild == 1) ) {
1092 if ( !CheckKinematicsOnChild() ) {
234f6d32 1093 delete[] pParent;
0f6ee828 1094 return 0;
1095 }
aea21c57 1096 }
1097
1098
f913ec4f 1099 for (i = 0; i < np; i++) {
8d2cd130 1100 Int_t trackIt = 0;
8507138f 1101 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 1102 kf = CheckPDGCode(iparticle->GetPdgCode());
1103 Int_t ks = iparticle->GetStatusCode();
1104 Int_t km = iparticle->GetFirstMother();
06956108 1105 if (
1106 (((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
1107 ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
1108 )
1109 {
1110 nc++;
8d2cd130 1111 if (ks == 1) trackIt = 1;
1112 Int_t ipa = iparticle->GetFirstMother()-1;
1113
1114 iparent = (ipa > -1) ? pParent[ipa] : -1;
1115
1116//
1117// store track information
1118 p[0] = iparticle->Px();
1119 p[1] = iparticle->Py();
1120 p[2] = iparticle->Pz();
a920faf9 1121 p[3] = iparticle->Energy();
1406f599 1122
a920faf9 1123
2590ccf9 1124 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1125 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1126 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1127
21391258 1128 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
a920faf9 1129
1130 PushTrack(fTrackIt*trackIt, iparent, kf,
1131 p[0], p[1], p[2], p[3],
1132 origin[0], origin[1], origin[2], tof,
1133 polar[0], polar[1], polar[2],
1134 kPPrimary, nt, 1., ks);
dbd64db6 1135 fNprimaries++;
8d2cd130 1136 KeepTrack(nt);
1137 pParent[i] = nt;
f913ec4f 1138 SetHighWaterMark(nt);
1139
8d2cd130 1140 } // select particle
1141 } // particle loop
1142
234f6d32 1143 delete[] pParent;
e0e89f40 1144
f913ec4f 1145 return 1;
8d2cd130 1146}
1147
1148
1149void AliGenPythia::FinishRun()
1150{
1151// Print x-section summary
1152 fPythia->Pystat(1);
2779fc64 1153
1154 if (fNev > 0.) {
1155 fQ /= fNev;
1156 fX1 /= fNev;
1157 fX2 /= fNev;
1158 }
1159
8d2cd130 1160 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1161 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1162}
1163
7184e472 1164void AliGenPythia::AdjustWeights() const
8d2cd130 1165{
1166// Adjust the weights after generation of all events
1167//
e2bddf81 1168 if (gAlice) {
1169 TParticle *part;
1170 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1171 for (Int_t i=0; i<ntrack; i++) {
1172 part= gAlice->GetMCApp()->Particle(i);
1173 part->SetWeight(part->GetWeight()*fKineBias);
1174 }
8d2cd130 1175 }
1176}
1177
20e47f08 1178void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
8d2cd130 1179{
1180// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1181
1a626d4e 1182 fAProjectile = a1;
1183 fATarget = a2;
66f02a7f 1184 fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
1d568bc2 1185 fSetNuclei = kTRUE;
8d2cd130 1186}
1187
1188
1189void AliGenPythia::MakeHeader()
1190{
7184e472 1191//
1192// Make header for the simulated event
1193//
183a5ca9 1194 if (gAlice) {
1195 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1196 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1197 }
1198
8d2cd130 1199// Builds the event header, to be called after each event
e5c87a3d 1200 if (fHeader) delete fHeader;
1201 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1202//
1203// Event type
800be8b7 1204 if(fProcDiff>0){
1205 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1206 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1207 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1208 }
1209 else
e5c87a3d 1210 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1211//
1212// Number of trials
e5c87a3d 1213 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1214//
1215// Event Vertex
d25cfd65 1216 fHeader->SetPrimaryVertex(fVertex);
21391258 1217 fHeader->SetInteractionTime(fTime+fEventTime);
dbd64db6 1218//
1219// Number of primaries
1220 fHeader->SetNProduced(fNprimaries);
8d2cd130 1221//
1222// Jets that have triggered
f913ec4f 1223
9dfe63b3 1224 //Need to store jets for b-jet studies too!
2f405d65 1225 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
8d2cd130 1226 {
1227 Int_t ntrig, njet;
1228 Float_t jets[4][10];
1229 GetJets(njet, ntrig, jets);
9ff6c04c 1230
8d2cd130 1231
1232 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1233 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1234 jets[3][i]);
1235 }
1236 }
5fa4b20b 1237//
1238// Copy relevant information from external header, if present.
1239//
1240 Float_t uqJet[4];
1241
1242 if (fRL) {
1243 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1244 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1245 {
1246 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1247
1248
1249 exHeader->TriggerJet(i, uqJet);
1250 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1251 }
1252 }
1253//
1254// Store quenching parameters
1255//
1256 if (fQuench){
b6262a45 1257 Double_t z[4] = {0.};
1258 Double_t xp = 0.;
1259 Double_t yp = 0.;
1260
7c21f297 1261 if (fQuench == 1) {
1262 // Pythia::Quench()
1263 fPythia->GetQuenchingParameters(xp, yp, z);
1bab4b79 1264 } else if (fQuench == 2){
7c21f297 1265 // Pyquen
1266 Double_t r1 = PARIMP.rb1;
1267 Double_t r2 = PARIMP.rb2;
1268 Double_t b = PARIMP.b1;
1269 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1270 Double_t phi = PARIMP.psib1;
1271 xp = r * TMath::Cos(phi);
1272 yp = r * TMath::Sin(phi);
1273
1bab4b79 1274 } else if (fQuench == 4) {
1275 // QPythia
5831e75f 1276 Double_t xy[2];
e9719084 1277 Double_t i0i1[2];
1278 AliFastGlauber::Instance()->GetSavedXY(xy);
1279 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
5831e75f 1280 xp = xy[0];
1281 yp = xy[1];
e6fe9b82 1282 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
7c21f297 1283 }
1bab4b79 1284
7c21f297 1285 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1286 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1bab4b79 1287 }
beac474c 1288//
1289// Store pt^hard
1290 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 1291//
cf57b268 1292// Pass header
5fa4b20b 1293//
cf57b268 1294 AddHeader(fHeader);
4c4eac97 1295 fHeader = 0x0;
8d2cd130 1296}
cf57b268 1297
c2fc9d31 1298Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
8d2cd130 1299{
1300// Check the kinematic trigger condition
1301//
1302 Double_t eta[2];
1303 eta[0] = jet1->Eta();
1304 eta[1] = jet2->Eta();
1305 Double_t phi[2];
1306 phi[0] = jet1->Phi();
1307 phi[1] = jet2->Phi();
1308 Int_t pdg[2];
1309 pdg[0] = jet1->GetPdgCode();
1310 pdg[1] = jet2->GetPdgCode();
1311 Bool_t triggered = kFALSE;
1312
2f405d65 1313 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8d2cd130 1314 Int_t njets = 0;
1315 Int_t ntrig = 0;
1316 Float_t jets[4][10];
1317//
1318// Use Pythia clustering on parton level to determine jet axis
1319//
1320 GetJets(njets, ntrig, jets);
1321
76d6ba9a 1322 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1323//
1324 } else {
1325 Int_t ij = 0;
1326 Int_t ig = 1;
1327 if (pdg[0] == kGamma) {
1328 ij = 1;
1329 ig = 0;
1330 }
1331 //Check eta range first...
1332 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1333 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1334 {
1335 //Eta is okay, now check phi range
1336 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1337 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1338 {
1339 triggered = kTRUE;
1340 }
1341 }
1342 }
1343 return triggered;
1344}
aea21c57 1345
1346
aea21c57 1347
7184e472 1348Bool_t AliGenPythia::CheckKinematicsOnChild(){
1349//
1350//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1351//
aea21c57 1352 Bool_t checking = kFALSE;
1353 Int_t j, kcode, ks, km;
1354 Int_t nPartAcc = 0; //number of particles in the acceptance range
1355 Int_t numberOfAcceptedParticles = 1;
1356 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1357 Int_t npart = fParticles.GetEntriesFast();
aea21c57 1358
0f6ee828 1359 for (j = 0; j<npart; j++) {
8507138f 1360 TParticle * jparticle = (TParticle *) fParticles.At(j);
aea21c57 1361 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1362 ks = jparticle->GetStatusCode();
1363 km = jparticle->GetFirstMother();
1364
1365 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1366 nPartAcc++;
1367 }
0f6ee828 1368 if( numberOfAcceptedParticles <= nPartAcc){
1369 checking = kTRUE;
1370 break;
1371 }
aea21c57 1372 }
0f6ee828 1373
aea21c57 1374 return checking;
aea21c57 1375}
1376
5fa4b20b 1377void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1378{
1058d9df 1379 //
1380 // Load event into Pythia Common Block
1381 //
1382
1383 Int_t npart = stack -> GetNprimary();
1384 Int_t n0 = 0;
1385
1386 if (!flag) {
1387 (fPythia->GetPyjets())->N = npart;
1388 } else {
1389 n0 = (fPythia->GetPyjets())->N;
1390 (fPythia->GetPyjets())->N = n0 + npart;
1391 }
1392
1393
1394 for (Int_t part = 0; part < npart; part++) {
1395 TParticle *mPart = stack->Particle(part);
32d6ef7d 1396
1058d9df 1397 Int_t kf = mPart->GetPdgCode();
1398 Int_t ks = mPart->GetStatusCode();
1399 Int_t idf = mPart->GetFirstDaughter();
1400 Int_t idl = mPart->GetLastDaughter();
1401
1402 if (reHadr) {
1403 if (ks == 11 || ks == 12) {
1404 ks -= 10;
1405 idf = -1;
1406 idl = -1;
1407 }
32d6ef7d 1408 }
1409
1058d9df 1410 Float_t px = mPart->Px();
1411 Float_t py = mPart->Py();
1412 Float_t pz = mPart->Pz();
1413 Float_t e = mPart->Energy();
1414 Float_t m = mPart->GetCalcMass();
32d6ef7d 1415
1058d9df 1416
1417 (fPythia->GetPyjets())->P[0][part+n0] = px;
1418 (fPythia->GetPyjets())->P[1][part+n0] = py;
1419 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1420 (fPythia->GetPyjets())->P[3][part+n0] = e;
1421 (fPythia->GetPyjets())->P[4][part+n0] = m;
1422
1423 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1424 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1425 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1426 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1427 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1428 }
1429}
1430
c2fc9d31 1431void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1058d9df 1432{
1433 //
1434 // Load event into Pythia Common Block
1435 //
1436
1437 Int_t npart = stack -> GetEntries();
1438 Int_t n0 = 0;
1439
1440 if (!flag) {
1441 (fPythia->GetPyjets())->N = npart;
1442 } else {
1443 n0 = (fPythia->GetPyjets())->N;
1444 (fPythia->GetPyjets())->N = n0 + npart;
1445 }
1446
1447
1448 for (Int_t part = 0; part < npart; part++) {
1449 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
92847124 1450 if (!mPart) continue;
1451
1058d9df 1452 Int_t kf = mPart->GetPdgCode();
1453 Int_t ks = mPart->GetStatusCode();
1454 Int_t idf = mPart->GetFirstDaughter();
1455 Int_t idl = mPart->GetLastDaughter();
1456
1457 if (reHadr) {
92847124 1458 if (ks == 11 || ks == 12) {
1459 ks -= 10;
1460 idf = -1;
1461 idl = -1;
1462 }
8d2cd130 1463 }
1058d9df 1464
1465 Float_t px = mPart->Px();
1466 Float_t py = mPart->Py();
1467 Float_t pz = mPart->Pz();
1468 Float_t e = mPart->Energy();
1469 Float_t m = mPart->GetCalcMass();
1470
1471
1472 (fPythia->GetPyjets())->P[0][part+n0] = px;
1473 (fPythia->GetPyjets())->P[1][part+n0] = py;
1474 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1475 (fPythia->GetPyjets())->P[3][part+n0] = e;
1476 (fPythia->GetPyjets())->P[4][part+n0] = m;
1477
1478 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1479 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1480 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1481 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1482 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1483 }
8d2cd130 1484}
1485
5fa4b20b 1486
014a9521 1487void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1488{
1489//
1490// Calls the Pythia jet finding algorithm to find jets in the current event
1491//
1492//
8d2cd130 1493//
1494// Save jets
1495 Int_t n = fPythia->GetN();
1496
1497//
1498// Run Jet Finder
1499 fPythia->Pycell(njets);
1500 Int_t i;
1501 for (i = 0; i < njets; i++) {
1502 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1503 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1504 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1505 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1506
1507 jets[0][i] = px;
1508 jets[1][i] = py;
1509 jets[2][i] = pz;
1510 jets[3][i] = e;
1511 }
1512}
1513
1514
1515
1516void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1517{
1518//
1519// Calls the Pythia clustering algorithm to find jets in the current event
1520//
1521 Int_t n = fPythia->GetN();
1522 nJets = 0;
1523 nJetsTrig = 0;
1524 if (fJetReconstruction == kCluster) {
1525//
1526// Configure cluster algorithm
1527//
1528 fPythia->SetPARU(43, 2.);
1529 fPythia->SetMSTU(41, 1);
1530//
1531// Call cluster algorithm
1532//
1533 fPythia->Pyclus(nJets);
1534//
1535// Loading jets from common block
1536//
1537 } else {
592f8307 1538
8d2cd130 1539//
1540// Run Jet Finder
1541 fPythia->Pycell(nJets);
1542 }
1543
1544 Int_t i;
1545 for (i = 0; i < nJets; i++) {
1546 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1547 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1548 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1549 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1550 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1551 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1552 Float_t theta = TMath::ATan2(pt,pz);
1553 Float_t et = e * TMath::Sin(theta);
1554 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1555 if (
1556 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1557 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1558 et > fEtMinJet && et < fEtMaxJet
1559 )
1560 {
1561 jets[0][nJetsTrig] = px;
1562 jets[1][nJetsTrig] = py;
1563 jets[2][nJetsTrig] = pz;
1564 jets[3][nJetsTrig] = e;
1565 nJetsTrig++;
5fa4b20b 1566// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1567 } else {
1568// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1569 }
1570 }
1571}
1572
f913ec4f 1573void AliGenPythia::GetSubEventTime()
1574{
1575 // Calculates time of the next subevent
9d7108a7 1576 fEventTime = 0.;
1577 if (fEventsTime) {
1578 TArrayF &array = *fEventsTime;
1579 fEventTime = array[fCurSubEvent++];
1580 }
1581 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1582 return;
f913ec4f 1583}
8d2cd130 1584
40fe59d4 1585Bool_t AliGenPythia::IsInBarrel(const Float_t eta) const
1586{
1587 // Is particle in Central Barrel acceptance?
1588 // etamin=-etamax
1589 if( eta < fTriggerEta )
1590 return kTRUE;
1591 else
1592 return kFALSE;
1593}
1594
1595Bool_t AliGenPythia::IsInEMCAL(const Float_t phi, const Float_t eta) const
ec2c406e 1596{
1597 // Is particle in EMCAL acceptance?
ec2c406e 1598 // phi in degrees, etamin=-etamax
9fd8e520 1599 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1600 eta < fEMCALEta )
ec2c406e 1601 return kTRUE;
1602 else
1603 return kFALSE;
1604}
1605
40fe59d4 1606Bool_t AliGenPythia::IsInPHOS(const Float_t phi, const Float_t eta, const Int_t iparticle)
ec2c406e 1607{
1608 // Is particle in PHOS acceptance?
1609 // Acceptance slightly larger considered.
1610 // phi in degrees, etamin=-etamax
40fe59d4 1611 // iparticle is the index of the particle to be checked, for PHOS rotation case
1612
9fd8e520 1613 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1614 eta < fPHOSEta )
ec2c406e 1615 return kTRUE;
1616 else
40fe59d4 1617 {
1618 if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1619
ec2c406e 1620 return kFALSE;
40fe59d4 1621 }
ec2c406e 1622}
1623
40fe59d4 1624void AliGenPythia::RotatePhi(Bool_t& okdd)
90a236ce 1625{
40fe59d4 1626 //Rotate event in phi to enhance events in PHOS acceptance
1627
1628 if(fPHOSRotateCandidate < 0) return ;
1629
90a236ce 1630 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1631 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1632 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
800be8b7 1633 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
90a236ce 1634
1635 //calculate deltaphi
40fe59d4 1636 TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
90a236ce 1637 Double_t phphi = ph->Phi();
1638 Double_t deltaphi = phiPHOS - phphi;
40fe59d4 1639
90a236ce 1640
1641
1642 //loop for all particles and produce the phi rotation
8507138f 1643 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
90a236ce 1644 Double_t oldphi, newphi;
777e69b0 1645 Double_t newVx, newVy, r, vZ, time;
1646 Double_t newPx, newPy, pt, pz, e;
90a236ce 1647 for(Int_t i=0; i< np; i++) {
40fe59d4 1648 TParticle* iparticle = (TParticle *) fParticles.At(i);
1649 oldphi = iparticle->Phi();
1650 newphi = oldphi + deltaphi;
1651 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1652 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1653
1654 r = iparticle->R();
1655 newVx = r * TMath::Cos(newphi);
1656 newVy = r * TMath::Sin(newphi);
1657 vZ = iparticle->Vz(); // don't transform
1658 time = iparticle->T(); // don't transform
1659
1660 pt = iparticle->Pt();
1661 newPx = pt * TMath::Cos(newphi);
1662 newPy = pt * TMath::Sin(newphi);
1663 pz = iparticle->Pz(); // don't transform
1664 e = iparticle->Energy(); // don't transform
1665
1666 // apply rotation
1667 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1668 iparticle->SetMomentum(newPx, newPy, pz, e);
1669
90a236ce 1670 } //end particle loop
1671
40fe59d4 1672 // now let's check that we put correctly the candidate photon in PHOS
1673 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1674 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1675 if(IsInPHOS(phi,eta,-1))
1676 okdd = kTRUE;
1677
1678 // reset the value for next event
1679 fPHOSRotateCandidate = -1;
1680
90a236ce 1681}
ec2c406e 1682
1683
800be8b7 1684Bool_t AliGenPythia::CheckDiffraction()
1685{
1686 // use this method only with Perugia-0 tune!
1687
1688 // printf("AAA\n");
1689
1690 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1691
1692 Int_t iPart1=-1;
1693 Int_t iPart2=-1;
1694
1695 Double_t y1 = 1e10;
1696 Double_t y2 = -1e10;
1697
1698 const Int_t kNstable=20;
1699 const Int_t pdgStable[20] = {
1700 22, // Photon
1701 11, // Electron
1702 12, // Electron Neutrino
1703 13, // Muon
1704 14, // Muon Neutrino
1705 15, // Tau
1706 16, // Tau Neutrino
1707 211, // Pion
1708 321, // Kaon
1709 311, // K0
1710 130, // K0s
1711 310, // K0l
1712 2212, // Proton
1713 2112, // Neutron
1714 3122, // Lambda_0
1715 3112, // Sigma Minus
1716 3222, // Sigma Plus
1717 3312, // Xsi Minus
1718 3322, // Xsi0
1719 3334 // Omega
1720 };
1721
1722 for (Int_t i = 0; i < np; i++) {
1723 TParticle * part = (TParticle *) fParticles.At(i);
1724
1725 Int_t statusCode = part->GetStatusCode();
1726
1727 // Initial state particle
1728 if (statusCode != 1)
1729 continue;
1730
1731 Int_t pdg = TMath::Abs(part->GetPdgCode());
1732 Bool_t isStable = kFALSE;
1733 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1734 if (pdg == pdgStable[i1]) {
1735 isStable = kTRUE;
1736 break;
1737 }
1738 }
1739 if(!isStable)
1740 continue;
1741
1742 Double_t y = part->Y();
1743
1744 if (y < y1)
1745 {
1746 y1 = y;
1747 iPart1 = i;
1748 }
1749 if (y > y2)
1750 {
1751 y2 = y;
1752 iPart2 = i;
1753 }
1754 }
1755
1756 if(iPart1<0 || iPart2<0) return kFALSE;
1757
1758 y1=TMath::Abs(y1);
1759 y2=TMath::Abs(y2);
1760
1761 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1762 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1763
1764 Int_t pdg1 = part1->GetPdgCode();
1765 Int_t pdg2 = part2->GetPdgCode();
1766
1767
1768 Int_t iPart = -1;
1769 if (pdg1 == 2212 && pdg2 == 2212)
1770 {
1771 if(y1 > y2)
1772 iPart = iPart1;
1773 else if(y1 < y2)
1774 iPart = iPart2;
1775 else {
1776 iPart = iPart1;
1777 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1778 }
1779 }
1780 else if (pdg1 == 2212)
1781 iPart = iPart1;
1782 else if (pdg2 == 2212)
1783 iPart = iPart2;
1784
1785
1786
1787
1788
1789 Double_t M=-1.;
1790 if(iPart>0) {
1791 TParticle * part = (TParticle *) fParticles.At(iPart);
1792 Double_t E= part->Energy();
1793 Double_t P= part->P();
1794 M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
1795 }
1796
c5dfa3e4 1797 Double_t Mmin, Mmax, wSD, wDD, wND;
1798 if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1799
1800 if(M>-1 && M<Mmin) return kFALSE;
1801 if(M>Mmax) M=-1;
1802
1803 Int_t procType=fPythia->GetMSTI(1);
1804 Int_t proc0=2;
1805 if(procType== 94) proc0=1;
1806 if(procType== 92 || procType== 93) proc0=0;
1807
1808 Int_t proc=2;
1809 if(M>0) proc=0;
1810 else if(proc0==1) proc=1;
1811
1812 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1813 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
800be8b7 1814
1815
1816 // if(proc==1 || proc==2) return kFALSE;
1817
c5dfa3e4 1818 if(proc!=0) {
1819 if(proc0!=0) fProcDiff = procType;
1820 else fProcDiff = 95;
1821 return kTRUE;
1822 }
1823
1824 if(wSD<0) AliError("wSD<0 ! \n");
1825
1826 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1827
1828 // printf("iPart = %d\n", iPart);
1829
1830 if(iPart==iPart1) fProcDiff=93;
1831 else if(iPart==iPart2) fProcDiff=92;
1832 else {
1833 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1834
800be8b7 1835 }
1836
c5dfa3e4 1837 return kTRUE;
1838}
1839
1840
1841
1842Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1843 Double_t &wSD, Double_t &wDD, Double_t &wND)
1844{
1845
1846 // 900 GeV
1847 if(TMath::Abs(fEnergyCMS-900)<1 ){
1848
1849const Int_t nbin=400;
1850Double_t bin[]={
18511.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
18524.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
18537.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
185410.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
185513.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
185615.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
185718.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
185821.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
185924.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
186027.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
186130.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
186233.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
186336.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
186439.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
186542.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
186645.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
186748.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
186851.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
186954.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
187057.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
187160.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
187263.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
187366.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
187469.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
187572.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
187675.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
187778.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
187881.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
187984.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
188087.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
188190.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
188293.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
188396.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
188499.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1885102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1886105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1887108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1888111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1889114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1890117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1891120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1892123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1893126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
1894129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
1895132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
1896135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
1897138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
1898141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
1899144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
1900147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
1901150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
1902153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
1903156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
1904159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
1905162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
1906165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
1907168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
1908171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
1909174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
1910177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
1911180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
1912183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
1913186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
1914189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
1915192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
1916195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
1917198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1918Double_t w[]={
19191.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
19200.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
19210.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
19220.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
19230.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
19240.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
19250.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
19260.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
19270.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
19280.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
19290.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
19300.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
19310.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
19320.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
19330.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
19340.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
19350.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
19360.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
19370.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
19380.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
19390.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
19400.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
19410.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
19420.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
19430.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
19440.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
19450.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
19460.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
19470.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
19480.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
19490.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
19500.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
19510.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
19520.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
19530.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
19540.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
19550.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
19560.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
19570.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
19580.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
19590.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
19600.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
19610.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
19620.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
19630.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
19640.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
19650.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
19660.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
19670.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
19680.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
19690.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
19700.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
19710.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
19720.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
19730.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
19740.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
19750.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
19760.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
19770.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
19780.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
19790.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
19800.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
19810.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
19820.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
19830.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
19840.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
19850.386112, 0.364314, 0.434375, 0.334629};
1986wDD = 0.379611;
1987wND = 0.496961;
1988wSD = -1;
1989
1990 Mmin = bin[0];
1991 Mmax = bin[nbin];
1992 if(M<Mmin || M>Mmax) return kTRUE;
1993
7873275c 1994 Int_t ibin=nbin-1;
93a60586 1995 for(Int_t i=1; i<=nbin; i++)
2ff57420 1996 if(M<=bin[i]) {
1997 ibin=i-1;
1998 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
800be8b7 1999 break;
2000 }
c5dfa3e4 2001 wSD=w[ibin];
2002 return kTRUE;
2003 }
2004 else if(TMath::Abs(fEnergyCMS-2760)<1 ){
800be8b7 2005
c5dfa3e4 2006const Int_t nbin=400;
2007Double_t bin[]={
20081.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
20094.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
20107.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
201110.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
201213.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
201315.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
201418.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
201521.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
201624.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
201727.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
201830.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
201933.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
202036.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
202139.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
202242.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
202345.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
202448.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
202551.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
202654.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
202757.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
202860.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
202963.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
203066.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
203169.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
203272.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
203375.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
203478.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
203581.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
203684.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
203787.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
203890.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
203993.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
204096.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
204199.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2042102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2043105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2044108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2045111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2046114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2047117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2048120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2049123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2050126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2051129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2052132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2053135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2054138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2055141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2056144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2057147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2058150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2059153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2060156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2061159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2062162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2063165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2064168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2065171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2066174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2067177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2068180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2069183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2070186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2071189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2072192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2073195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2074198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2075Double_t w[]={
20761.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
20770.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
20780.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
20790.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
20800.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
20810.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
20820.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
20830.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
20840.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
20850.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
20860.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
20870.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
20880.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
20890.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
20900.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
20910.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
20920.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
20930.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
20940.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
20950.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
20960.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
20970.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
20980.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
20990.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
21000.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
21010.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
21020.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
21030.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
21040.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
21050.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
21060.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
21070.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
21080.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
21090.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
21100.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
21110.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
21120.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
21130.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
21140.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
21150.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
21160.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
21170.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
21180.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
21190.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
21200.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
21210.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
21220.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
21230.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
21240.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
21250.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
21260.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
21270.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
21280.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
21290.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
21300.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
21310.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
21320.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
21330.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
21340.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
21350.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
21360.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
21370.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
21380.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
21390.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
21400.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
21410.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
21420.201712, 0.242225, 0.255565, 0.258738};
2143wDD = 0.512813;
2144wND = 0.518820;
2145wSD = -1;
2146
2147 Mmin = bin[0];
2148 Mmax = bin[nbin];
2149 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2150
c5dfa3e4 2151 Int_t ibin=nbin-1;
2152 for(Int_t i=1; i<=nbin; i++)
2153 if(M<=bin[i]) {
2154 ibin=i-1;
2155 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2156 break;
2157 }
2158 wSD=w[ibin];
2159 return kTRUE;
2160 }
800be8b7 2161
800be8b7 2162
c5dfa3e4 2163 else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2164const Int_t nbin=400;
2165Double_t bin[]={
21661.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
21674.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
21687.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
216910.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
217013.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
217115.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
217218.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
217321.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
217424.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
217527.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
217630.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
217733.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
217836.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
217939.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
218042.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
218145.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
218248.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
218351.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
218454.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
218557.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
218660.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
218763.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
218866.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
218969.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
219072.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
219175.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
219278.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
219381.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
219484.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
219587.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
219690.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
219793.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
219896.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
219999.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2200102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2201105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2202108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2203111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2204114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2205117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2206120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2207123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2208126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2209129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2210132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2211135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2212138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2213141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2214144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2215147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2216150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2217153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2218156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2219159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2220162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2221165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2222168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2223171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2224174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2225177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2226180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2227183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2228186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2229189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2230192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2231195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2232198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2233Double_t w[]={
22341.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
22350.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
22360.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
22370.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
22380.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
22390.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
22400.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
22410.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
22420.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
22430.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
22440.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
22450.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
22460.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
22470.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
22480.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
22490.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
22500.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
22510.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
22520.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
22530.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
22540.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
22550.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
22560.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
22570.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
22580.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
22590.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
22600.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
22610.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
22620.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
22630.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
22640.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
22650.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
22660.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
22670.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
22680.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
22690.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
22700.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
22710.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
22720.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
22730.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
22740.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
22750.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
22760.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
22770.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
22780.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
22790.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
22800.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
22810.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
22820.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
22830.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
22840.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
22850.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
22860.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
22870.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
22880.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
22890.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
22900.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
22910.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
22920.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
22930.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
22940.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
22950.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
22960.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
22970.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
22980.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
22990.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
23000.175006, 0.223482, 0.202706, 0.218108};
2301wDD = 0.207705;
2302wND = 0.289628;
2303wSD = -1;
2304
2305 Mmin = bin[0];
2306 Mmax = bin[nbin];
2307
2308 if(M<Mmin || M>Mmax) return kTRUE;
800be8b7 2309
c5dfa3e4 2310 Int_t ibin=nbin-1;
2311 for(Int_t i=1; i<=nbin; i++)
2312 if(M<=bin[i]) {
2313 ibin=i-1;
2314 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2315 break;
2316 }
2317 wSD=w[ibin];
800be8b7 2318 return kTRUE;
c5dfa3e4 2319 }
2320
2321 return kFALSE;
800be8b7 2322}
06956108 2323
2324Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2325{
2326// Check if this is a heavy flavor decay product
2327 TParticle * part = (TParticle *) fParticles.At(ipart);
2328 Int_t mpdg = TMath::Abs(part->GetPdgCode());
2329 Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2330 //
2331 // Light hadron
2332 if (mfl >= 4 && mfl < 6) return kTRUE;
2333
2334 Int_t imo = part->GetFirstMother()-1;
2335 TParticle* pm = part;
2336 //
2337 // Heavy flavor hadron produced by generator
2338 while (imo > -1) {
2339 pm = (TParticle*)fParticles.At(imo);
2340 mpdg = TMath::Abs(pm->GetPdgCode());
2341 mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2342 if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2343 imo = pm->GetFirstMother()-1;
2344 }
2345 return kFALSE;
2346}
40fe59d4 2347
2348Bool_t AliGenPythia::CheckDetectorAcceptance(const Float_t phi, const Float_t eta, const Int_t iparticle)
2349{
2350 // check the eta/phi correspond to the detectors acceptance
2351 // iparticle is the index of the particle to be checked, for PHOS rotation case
2352 if (fCheckPHOS && IsInPHOS (phi,eta,iparticle)) return kTRUE;
2353 else if(fCheckEMCAL && IsInEMCAL (phi,eta)) return kTRUE;
2354 else if(fCheckBarrel && IsInBarrel( eta)) return kTRUE;
2355 else return kFALSE;
2356}
2357
2358Bool_t AliGenPythia::TriggerOnSelectedParticles(const Int_t np)
2359{
2360 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2361 // implemented primaryly for kPyJets, but extended to any kind of process.
2362
2363 //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2364 // fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
2365
2366 Bool_t ok = kFALSE;
2367 for (Int_t i=0; i< np; i++) {
2368
2369 TParticle* iparticle = (TParticle *) fParticles.At(i);
2370
2371 Int_t pdg = iparticle->GetPdgCode();
2372 Int_t status = iparticle->GetStatusCode();
2373 Int_t imother = iparticle->GetFirstMother() - 1;
2374
2375 TParticle* pmother = 0x0;
2376 Int_t momStatus = -1;
2377 Int_t momPdg = -1;
2378 if(imother > 0 ){
2379 pmother = (TParticle *) fParticles.At(imother);
2380 momStatus = pmother->GetStatusCode();
2381 momPdg = pmother->GetPdgCode();
2382 }
2383
2384 ok = kFALSE;
2385
2386 //
2387 // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2388 //
2389 // Hadron
2390 if (fHadronInCalo && status == 1)
2391 {
2392 if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
2393 // (in case neutral mesons were declared stable)
2394 ok = kTRUE;
2395 }
2396 //Electron
2397 else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2398 {
2399 ok = kTRUE;
2400 }
2401 //Fragmentation photon
2402 else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2403 {
2404 if(momStatus != 11) ok = kTRUE ; // No photon from hadron decay
2405 }
2406 // Decay photon
2407 else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2408 {
2409 if( momStatus == 11)
2410 {
2411 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2412 // pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2413 ok = kTRUE ; // photon from hadron decay
2414
2415 //In case only decays from pi0 or eta requested
2416 if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2417 if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2418 }
2419
2420 }
2421 // Pi0 or Eta particle
2422 else if ((fPi0InCalo || fEtaInCalo))
2423 {
2424 if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2425
2426 if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2427 {
2428 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2429 ok = kTRUE;
2430 }
2431 else if (fEtaInCalo && pdg == 221)
2432 {
2433 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2434 ok = kTRUE;
2435 }
2436
2437 }// pi0 or eta
2438
2439 //
2440 // Check that the selected particle is in the calorimeter acceptance
2441 //
2442 if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2443 {
2444 //Just check if the selected particle falls in the acceptance
2445 if(!fForceNeutralMeson2PhotonDecay )
2446 {
2447 //printf("\t Check acceptance! \n");
2448 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2449 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2450
2451 if(CheckDetectorAcceptance(phi,eta,i))
2452 {
2453 ok =kTRUE;
2454 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));
2455 //printf("\t Accept \n");
2456 break;
2457 }
2458 else ok = kFALSE;
2459 }
2460 //Mesons have several decay modes, select only those decaying into 2 photons
2461 else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2462 {
2463 // In case we want the pi0/eta trigger,
2464 // check the decay mode (2 photons)
2465
2466 //printf("\t Force decay 2 gamma\n");
2467
2468 Int_t ndaughters = iparticle->GetNDaughters();
2469 if(ndaughters != 2){
2470 ok=kFALSE;
2471 continue;
2472 }
2473
2474 TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2475 TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2476 if(!d1 || !d2) {
2477 ok=kFALSE;
2478 continue;
2479 }
2480
2481 //iparticle->Print();
2482 //d1->Print();
2483 //d2->Print();
2484
2485 Int_t pdgD1 = d1->GetPdgCode();
2486 Int_t pdgD2 = d2->GetPdgCode();
2487 //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2488 //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2489
2490 if(pdgD1 != 22 || pdgD2 != 22){
2491 ok = kFALSE;
2492 continue;
2493 }
2494
2495 //printf("\t accept decay\n");
2496
2497 //Trigger on the meson, not on the daughter
2498 if(!fDecayPhotonInCalo){
2499
2500 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2501 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2502
2503 if(CheckDetectorAcceptance(phi,eta,i))
2504 {
2505 //printf("\t Accept meson pdg %d\n",pdg);
2506 ok =kTRUE;
2507 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));
2508 break;
2509 } else {
2510 ok=kFALSE;
2511 continue;
2512 }
2513 }
2514
2515 //printf("Check daughters acceptance\n");
2516
2517 //Trigger on the meson daughters
2518 //Photon 1
2519 Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2520 Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
2521 if(d1->Pt() > fTriggerParticleMinPt)
2522 {
2523 //printf("\t Check acceptance photon 1! \n");
2524 if(CheckDetectorAcceptance(phi,eta,i))
2525 {
2526 //printf("\t Accept Photon 1\n");
2527 ok =kTRUE;
2528 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));
2529 break;
2530 }
2531 else ok = kFALSE;
2532 } // pt cut
2533 else ok = kFALSE;
2534
2535 //Photon 2
2536 phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2537 eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2538
2539 if(d2->Pt() > fTriggerParticleMinPt)
2540 {
2541 //printf("\t Check acceptance photon 2! \n");
2542 if(CheckDetectorAcceptance(phi,eta,i))
2543 {
2544 //printf("\t Accept Photon 2\n");
2545 ok =kTRUE;
2546 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));
2547 break;
2548 }
2549 else ok = kFALSE;
2550 } // pt cut
2551 else ok = kFALSE;
2552 } // force 2 photon daughters in pi0/eta decays
2553 else ok = kFALSE;
2554 } else ok = kFALSE; // check acceptance
2555 } // primary loop
2556
2557 //
2558 // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2559 // A particle passing all trigger conditions except phi position in PHOS, is used as reference
2560 //
2561 if(fCheckPHOSeta)
2562 {
2563 RotatePhi(ok);
2564 }
2565
2566 return ok;
2567}
2568