git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Tuned crossections according to the ALICE measurements (Martin)
[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.),
e6fe9b82 87 fImpact(0.),
e8a8adcd 88 fPtKick(1.),
89 fFullEvent(kTRUE),
90 fDecayer(new AliDecayerPythia()),
91 fDebugEventFirst(-1),
92 fDebugEventLast(-1),
93 fEtMinJet(0.),
94 fEtMaxJet(1.e4),
95 fEtaMinJet(-20.),
96 fEtaMaxJet(20.),
97 fPhiMinJet(0.),
98 fPhiMaxJet(2.* TMath::Pi()),
99 fJetReconstruction(kCell),
100 fEtaMinGamma(-20.),
101 fEtaMaxGamma(20.),
102 fPhiMinGamma(0.),
103 fPhiMaxGamma(2. * TMath::Pi()),
104 fPycellEtaMax(2.),
105 fPycellNEta(274),
106 fPycellNPhi(432),
107 fPycellThreshold(0.),
108 fPycellEtSeed(4.),
109 fPycellMinEtJet(10.),
110 fPycellMaxRadius(1.),
111 fStackFillOpt(kFlavorSelection),
112 fFeedDownOpt(kTRUE),
113 fFragmentation(kTRUE),
114 fSetNuclei(kFALSE),
115 fNewMIS(kFALSE),
116 fHFoff(kFALSE),
20e47f08 117 fNucPdf(0),
e8a8adcd 118 fTriggerParticle(0),
119 fTriggerEta(0.9),
700b9416 120 fTriggerMultiplicity(0),
121 fTriggerMultiplicityEta(0),
38112f3f 122 fTriggerMultiplicityPtMin(0),
e8a8adcd 123 fCountMode(kCountAll),
124 fHeader(0),
125 fRL(0),
777e69b0 126 fkFileName(0),
9fd8e520 127 fFragPhotonInCalo(kFALSE),
ec2c406e 128 fPi0InCalo(kFALSE) ,
90a236ce 129 fPhotonInCalo(kFALSE),
9dfe63b3 130 fEleInEMCAL(kFALSE),
9fd8e520 131 fCheckEMCAL(kFALSE),
132 fCheckPHOS(kFALSE),
90a236ce 133 fCheckPHOSeta(kFALSE),
9fd8e520 134 fFragPhotonOrPi0MinPt(0),
90a236ce 135 fPhotonMinPt(0),
9dfe63b3 136 fElectronMinPt(0),
9fd8e520 137 fPHOSMinPhi(219.),
138 fPHOSMaxPhi(321.),
139 fPHOSEta(0.13),
140 fEMCALMinPhi(79.),
141 fEMCALMaxPhi(191.),
800be8b7 142 fEMCALEta(0.71),
143 fkTuneForDiff(0),
144 fProcDiff(0)
8d2cd130 145{
146// Default Constructor
e7c989e4 147 fEnergyCMS = 5500.;
7cdba479 148 if (!AliPythiaRndm::GetPythiaRandom())
e8a8adcd 149 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 150}
151
152AliGenPythia::AliGenPythia(Int_t npart)
e8a8adcd 153 :AliGenMC(npart),
154 fProcess(kPyCharm),
efe3b1cd 155 fItune(-1),
e8a8adcd 156 fStrucFunc(kCTEQ5L),
e8a8adcd 157 fKineBias(0.),
158 fTrials(0),
159 fTrialsRun(0),
160 fQ(0.),
161 fX1(0.),
162 fX2(0.),
163 fEventTime(0.),
164 fInteractionRate(0.),
165 fTimeWindow(0.),
166 fCurSubEvent(0),
167 fEventsTime(0),
168 fNev(0),
169 fFlavorSelect(0),
170 fXsection(0.),
171 fPythia(0),
172 fPtHardMin(0.),
173 fPtHardMax(1.e4),
174 fYHardMin(-1.e10),
175 fYHardMax(1.e10),
176 fGinit(kTRUE),
177 fGfinal(kTRUE),
178 fHadronisation(kTRUE),
03358a32 179 fPatchOmegaDalitz(0),
e8a8adcd 180 fNpartons(0),
181 fReadFromFile(kFALSE),
182 fQuench(kFALSE),
cd07c39b 183 fQhat(0.),
184 fLength(0.),
e6fe9b82 185 fImpact(0.),
e8a8adcd 186 fPtKick(1.),
187 fFullEvent(kTRUE),
188 fDecayer(new AliDecayerPythia()),
189 fDebugEventFirst(-1),
190 fDebugEventLast(-1),
191 fEtMinJet(0.),
192 fEtMaxJet(1.e4),
193 fEtaMinJet(-20.),
194 fEtaMaxJet(20.),
195 fPhiMinJet(0.),
196 fPhiMaxJet(2.* TMath::Pi()),
197 fJetReconstruction(kCell),
198 fEtaMinGamma(-20.),
199 fEtaMaxGamma(20.),
200 fPhiMinGamma(0.),
201 fPhiMaxGamma(2. * TMath::Pi()),
202 fPycellEtaMax(2.),
203 fPycellNEta(274),
204 fPycellNPhi(432),
205 fPycellThreshold(0.),
206 fPycellEtSeed(4.),
207 fPycellMinEtJet(10.),
208 fPycellMaxRadius(1.),
209 fStackFillOpt(kFlavorSelection),
210 fFeedDownOpt(kTRUE),
211 fFragmentation(kTRUE),
212 fSetNuclei(kFALSE),
213 fNewMIS(kFALSE),
214 fHFoff(kFALSE),
20e47f08 215 fNucPdf(0),
e8a8adcd 216 fTriggerParticle(0),
217 fTriggerEta(0.9),
700b9416 218 fTriggerMultiplicity(0),
219 fTriggerMultiplicityEta(0),
38112f3f 220 fTriggerMultiplicityPtMin(0),
e8a8adcd 221 fCountMode(kCountAll),
222 fHeader(0),
223 fRL(0),
777e69b0 224 fkFileName(0),
9fd8e520 225 fFragPhotonInCalo(kFALSE),
ec2c406e 226 fPi0InCalo(kFALSE) ,
90a236ce 227 fPhotonInCalo(kFALSE),
9dfe63b3 228 fEleInEMCAL(kFALSE),
9fd8e520 229 fCheckEMCAL(kFALSE),
230 fCheckPHOS(kFALSE),
90a236ce 231 fCheckPHOSeta(kFALSE),
9fd8e520 232 fFragPhotonOrPi0MinPt(0),
90a236ce 233 fPhotonMinPt(0),
9dfe63b3 234 fElectronMinPt(0),
9fd8e520 235 fPHOSMinPhi(219.),
236 fPHOSMaxPhi(321.),
237 fPHOSEta(0.13),
238 fEMCALMinPhi(79.),
239 fEMCALMaxPhi(191.),
800be8b7 240 fEMCALEta(0.71),
241 fkTuneForDiff(0),
242 fProcDiff(0)
8d2cd130 243{
244// default charm production at 5. 5 TeV
245// semimuonic decay
246// structure function GRVHO
247//
e7c989e4 248 fEnergyCMS = 5500.;
8d2cd130 249 fName = "Pythia";
250 fTitle= "Particle Generator using PYTHIA";
8d2cd130 251 SetForceDecay();
8d2cd130 252 // Set random number generator
7cdba479 253 if (!AliPythiaRndm::GetPythiaRandom())
254 AliPythiaRndm::SetPythiaRandom(GetRandom());
76d6ba9a 255 }
8d2cd130 256
8d2cd130 257AliGenPythia::~AliGenPythia()
258{
259// Destructor
9d7108a7 260 if(fEventsTime) delete fEventsTime;
261}
262
263void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
264{
265// Generate pileup using user specified rate
266 fInteractionRate = rate;
267 fTimeWindow = timewindow;
268 GeneratePileup();
269}
270
271void AliGenPythia::GeneratePileup()
272{
273// Generate sub events time for pileup
274 fEventsTime = 0;
275 if(fInteractionRate == 0.) {
276 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
277 return;
278 }
279
280 Int_t npart = NumberParticles();
281 if(npart < 0) {
282 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
283 return;
284 }
285
286 if(fEventsTime) delete fEventsTime;
287 fEventsTime = new TArrayF(npart);
288 TArrayF &array = *fEventsTime;
289 for(Int_t ipart = 0; ipart < npart; ipart++)
290 array[ipart] = 0.;
291
292 Float_t eventtime = 0.;
293 while(1)
294 {
295 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
296 if(eventtime > fTimeWindow) break;
297 array.Set(array.GetSize()+1);
298 array[array.GetSize()-1] = eventtime;
299 }
300
301 eventtime = 0.;
302 while(1)
303 {
304 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
305 if(TMath::Abs(eventtime) > fTimeWindow) break;
306 array.Set(array.GetSize()+1);
307 array[array.GetSize()-1] = eventtime;
308 }
309
310 SetNumberParticles(fEventsTime->GetSize());
8d2cd130 311}
312
592f8307 313void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
314 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
315{
316// Set pycell parameters
317 fPycellEtaMax = etamax;
318 fPycellNEta = neta;
319 fPycellNPhi = nphi;
320 fPycellThreshold = thresh;
321 fPycellEtSeed = etseed;
322 fPycellMinEtJet = minet;
323 fPycellMaxRadius = r;
324}
325
326
327
8d2cd130 328void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
329{
330 // Set a range of event numbers, for which a table
331 // of generated particle will be printed
332 fDebugEventFirst = eventFirst;
333 fDebugEventLast = eventLast;
334 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
335}
336
337void AliGenPythia::Init()
338{
339// Initialisation
340
341 SetMC(AliPythia::Instance());
b88f5cea 342 fPythia=(AliPythia*) fMCEvGen;
e2bddf81 343
8d2cd130 344//
345 fParentWeight=1./Float_t(fNpart);
346//
8d2cd130 347
348
349 fPythia->SetCKIN(3,fPtHardMin);
350 fPythia->SetCKIN(4,fPtHardMax);
351 fPythia->SetCKIN(7,fYHardMin);
352 fPythia->SetCKIN(8,fYHardMax);
353
20e47f08 354 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
8d2cd130 355 // Fragmentation?
356 if (fFragmentation) {
357 fPythia->SetMSTP(111,1);
358 } else {
359 fPythia->SetMSTP(111,0);
360 }
361
362
363// initial state radiation
364 fPythia->SetMSTP(61,fGinit);
365// final state radiation
366 fPythia->SetMSTP(71,fGfinal);
367// pt - kick
368 if (fPtKick > 0.) {
369 fPythia->SetMSTP(91,1);
688950ef 370 fPythia->SetPARP(91,fPtKick);
371 fPythia->SetPARP(93, 4. * fPtKick);
8d2cd130 372 } else {
373 fPythia->SetMSTP(91,0);
374 }
375
5fa4b20b 376
377 if (fReadFromFile) {
777e69b0 378 fRL = AliRunLoader::Open(fkFileName, "Partons");
5fa4b20b 379 fRL->LoadKinematics();
380 fRL->LoadHeader();
381 } else {
382 fRL = 0x0;
383 }
fdea4387 384 //
efe3b1cd 385 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
fdea4387 386 // Forward Paramters to the AliPythia object
387 fDecayer->SetForceDecay(fForceDecay);
beac474c 388// Switch off Heavy Flavors on request
389 if (fHFoff) {
51387949 390 // Maximum number of quark flavours used in pdf
beac474c 391 fPythia->SetMSTP(58, 3);
51387949 392 // Maximum number of flavors that can be used in showers
beac474c 393 fPythia->SetMSTJ(45, 3);
51387949 394 // Switch off g->QQbar splitting in decay table
395 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
beac474c 396 }
fdea4387 397
51387949 398 fDecayer->Init();
399
8d2cd130 400
401// Parent and Children Selection
402 switch (fProcess)
403 {
65f2626c 404 case kPyOldUEQ2ordered:
405 case kPyOldUEQ2ordered2:
406 case kPyOldPopcorn:
407 break;
8d2cd130 408 case kPyCharm:
409 case kPyCharmUnforced:
adf4d898 410 case kPyCharmPbPbMNR:
aabc7187 411 case kPyCharmpPbMNR:
e0e89f40 412 case kPyCharmppMNR:
413 case kPyCharmppMNRwmi:
8d2cd130 414 fParentSelect[0] = 411;
415 fParentSelect[1] = 421;
416 fParentSelect[2] = 431;
417 fParentSelect[3] = 4122;
9ccc3d0e 418 fParentSelect[4] = 4232;
419 fParentSelect[5] = 4132;
420 fParentSelect[6] = 4332;
8d2cd130 421 fFlavorSelect = 4;
422 break;
adf4d898 423 case kPyD0PbPbMNR:
424 case kPyD0pPbMNR:
425 case kPyD0ppMNR:
8d2cd130 426 fParentSelect[0] = 421;
427 fFlavorSelect = 4;
428 break;
90d7b703 429 case kPyDPlusPbPbMNR:
430 case kPyDPluspPbMNR:
431 case kPyDPlusppMNR:
432 fParentSelect[0] = 411;
433 fFlavorSelect = 4;
434 break;
e0e89f40 435 case kPyDPlusStrangePbPbMNR:
436 case kPyDPlusStrangepPbMNR:
437 case kPyDPlusStrangeppMNR:
438 fParentSelect[0] = 431;
439 fFlavorSelect = 4;
440 break;
68504d42 441 case kPyLambdacppMNR:
442 fParentSelect[0] = 4122;
443 fFlavorSelect = 4;
444 break;
8d2cd130 445 case kPyBeauty:
9dfe63b3 446 case kPyBeautyJets:
adf4d898 447 case kPyBeautyPbPbMNR:
448 case kPyBeautypPbMNR:
449 case kPyBeautyppMNR:
e0e89f40 450 case kPyBeautyppMNRwmi:
8d2cd130 451 fParentSelect[0]= 511;
452 fParentSelect[1]= 521;
453 fParentSelect[2]= 531;
454 fParentSelect[3]= 5122;
455 fParentSelect[4]= 5132;
456 fParentSelect[5]= 5232;
457 fParentSelect[6]= 5332;
458 fFlavorSelect = 5;
459 break;
460 case kPyBeautyUnforced:
461 fParentSelect[0] = 511;
462 fParentSelect[1] = 521;
463 fParentSelect[2] = 531;
464 fParentSelect[3] = 5122;
465 fParentSelect[4] = 5132;
466 fParentSelect[5] = 5232;
467 fParentSelect[6] = 5332;
468 fFlavorSelect = 5;
469 break;
470 case kPyJpsiChi:
471 case kPyJpsi:
472 fParentSelect[0] = 443;
473 break;
f529e69b 474 case kPyMbDefault:
0bd3d7c5 475 case kPyMbAtlasTuneMC09:
8d2cd130 476 case kPyMb:
04081a91 477 case kPyMbWithDirectPhoton:
8d2cd130 478 case kPyMbNonDiffr:
d7de4a67 479 case kPyMbMSEL1:
8d2cd130 480 case kPyJets:
481 case kPyDirectGamma:
e33acb67 482 case kPyLhwgMb:
8d2cd130 483 break;
589380c6 484 case kPyW:
0f6ee828 485 case kPyZ:
589380c6 486 break;
8d2cd130 487 }
488//
592f8307 489//
490// JetFinder for Trigger
491//
492// Configure detector (EMCAL like)
493//
d7de4a67 494 fPythia->SetPARU(51, fPycellEtaMax);
495 fPythia->SetMSTU(51, fPycellNEta);
496 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 497//
498// Configure Jet Finder
499//
d7de4a67 500 fPythia->SetPARU(58, fPycellThreshold);
501 fPythia->SetPARU(52, fPycellEtSeed);
502 fPythia->SetPARU(53, fPycellMinEtJet);
503 fPythia->SetPARU(54, fPycellMaxRadius);
504 fPythia->SetMSTU(54, 2);
592f8307 505//
8d2cd130 506// This counts the total number of calls to Pyevnt() per run.
507 fTrialsRun = 0;
508 fQ = 0.;
509 fX1 = 0.;
510 fX2 = 0.;
511 fNev = 0 ;
512//
1d568bc2 513//
514//
8d2cd130 515 AliGenMC::Init();
1d568bc2 516//
517//
518//
519 if (fSetNuclei) {
520 fDyBoost = 0;
521 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
522 }
d7de4a67 523
cd07c39b 524 fPythia->SetPARJ(200, 0.0);
525 fPythia->SetPARJ(199, 0.0);
526 fPythia->SetPARJ(198, 0.0);
527 fPythia->SetPARJ(197, 0.0);
528
529 if (fQuench == 1) {
5fa4b20b 530 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 531 }
3a709cfa 532
b976f7d8 533 if (fQuench == 3) {
534 // Nestor's change of the splittings
535 fPythia->SetPARJ(200, 0.8);
536 fPythia->SetMSTJ(41, 1); // QCD radiation only
537 fPythia->SetMSTJ(42, 2); // angular ordering
538 fPythia->SetMSTJ(44, 2); // option to run alpha_s
539 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
540 fPythia->SetMSTJ(50, 0); // No coherence in first branching
541 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
cd07c39b 542 } else if (fQuench == 4) {
543 // Armesto-Cunqueiro-Salgado change of the splittings.
544 AliFastGlauber* glauber = AliFastGlauber::Instance();
545 glauber->Init(2);
546 //read and store transverse almonds corresponding to differnt
547 //impact parameters.
cd07c39b 548 glauber->SetCentralityClass(0.,0.1);
549 fPythia->SetPARJ(200, 1.);
550 fPythia->SetPARJ(198, fQhat);
551 fPythia->SetPARJ(199, fLength);
cd07c39b 552 fPythia->SetMSTJ(42, 2); // angular ordering
553 fPythia->SetMSTJ(44, 2); // option to run alpha_s
cd07c39b 554 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
b976f7d8 555 }
8d2cd130 556}
557
558void AliGenPythia::Generate()
559{
560// Generate one event
19ef8cf0 561 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
8d2cd130 562 fDecayer->ForceDecay();
563
13cca24a 564 Double_t polar[3] = {0,0,0};
565 Double_t origin[3] = {0,0,0};
566 Double_t p[4];
8d2cd130 567// converts from mm/c to s
568 const Float_t kconv=0.001/2.999792458e8;
569//
570 Int_t nt=0;
571 Int_t jev=0;
572 Int_t j, kf;
573 fTrials=0;
f913ec4f 574 fEventTime = 0.;
575
2590ccf9 576
8d2cd130 577
578 // Set collision vertex position
2590ccf9 579 if (fVertexSmear == kPerEvent) Vertex();
580
8d2cd130 581// event loop
582 while(1)
583 {
32d6ef7d 584//
5fa4b20b 585// Produce event
32d6ef7d 586//
32d6ef7d 587//
5fa4b20b 588// Switch hadronisation off
589//
590 fPythia->SetMSTJ(1, 0);
cd07c39b 591
592 if (fQuench ==4){
593 Double_t bimp;
594 // Quenching comes through medium-modified splitting functions.
595 AliFastGlauber::Instance()->GetRandomBHard(bimp);
e6fe9b82 596 fPythia->SetPARJ(197, bimp);
597 fImpact = bimp;
6c43eccb 598 fPythia->Qpygin0();
cd07c39b 599 }
32d6ef7d 600//
5fa4b20b 601// Either produce new event or read partons from file
602//
603 if (!fReadFromFile) {
beac474c 604 if (!fNewMIS) {
605 fPythia->Pyevnt();
606 } else {
607 fPythia->Pyevnw();
608 }
5fa4b20b 609 fNpartons = fPythia->GetN();
610 } else {
33c3c91a 611 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
612 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
5fa4b20b 613 fPythia->SetN(0);
614 LoadEvent(fRL->Stack(), 0 , 1);
615 fPythia->Pyedit(21);
616 }
617
32d6ef7d 618//
619// Run quenching routine
620//
5fa4b20b 621 if (fQuench == 1) {
622 fPythia->Quench();
623 } else if (fQuench == 2){
624 fPythia->Pyquen(208., 0, 0.);
b976f7d8 625 } else if (fQuench == 3) {
626 // Quenching is via multiplicative correction of the splittings
5fa4b20b 627 }
b976f7d8 628
32d6ef7d 629//
5fa4b20b 630// Switch hadronisation on
32d6ef7d 631//
20e47f08 632 if (fHadronisation) {
633 fPythia->SetMSTJ(1, 1);
5fa4b20b 634//
635// .. and perform hadronisation
aea21c57 636// printf("Calling hadronisation %d\n", fPythia->GetN());
03358a32 637
638 if (fPatchOmegaDalitz) {
639 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
640 fPythia->Pyexec();
641 fPythia->DalitzDecays();
642 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
643 }
644 fPythia->Pyexec();
20e47f08 645 }
8d2cd130 646 fTrials++;
8507138f 647 fPythia->ImportParticles(&fParticles,"All");
03358a32 648
df275cfa 649 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
e1cf8975 650 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
f64640b5 651
8d2cd130 652//
653//
654//
655 Int_t i;
656
dbd64db6 657 fNprimaries = 0;
8507138f 658 Int_t np = fParticles.GetEntriesFast();
5fa4b20b 659
7c21f297 660 if (np == 0) continue;
8d2cd130 661//
2590ccf9 662
8d2cd130 663//
664 Int_t* pParent = new Int_t[np];
665 Int_t* pSelected = new Int_t[np];
666 Int_t* trackIt = new Int_t[np];
5fa4b20b 667 for (i = 0; i < np; i++) {
8d2cd130 668 pParent[i] = -1;
669 pSelected[i] = 0;
670 trackIt[i] = 0;
671 }
672
673 Int_t nc = 0; // Total n. of selected particles
674 Int_t nParents = 0; // Selected parents
675 Int_t nTkbles = 0; // Trackable particles
f529e69b 676 if (fProcess != kPyMbDefault &&
677 fProcess != kPyMb &&
6d2ec66d 678 fProcess != kPyMbAtlasTuneMC09 &&
04081a91 679 fProcess != kPyMbWithDirectPhoton &&
f529e69b 680 fProcess != kPyJets &&
8d2cd130 681 fProcess != kPyDirectGamma &&
589380c6 682 fProcess != kPyMbNonDiffr &&
d7de4a67 683 fProcess != kPyMbMSEL1 &&
f529e69b 684 fProcess != kPyW &&
685 fProcess != kPyZ &&
686 fProcess != kPyCharmppMNRwmi &&
9dfe63b3 687 fProcess != kPyBeautyppMNRwmi &&
688 fProcess != kPyBeautyJets) {
8d2cd130 689
5fa4b20b 690 for (i = 0; i < np; i++) {
8507138f 691 TParticle* iparticle = (TParticle *) fParticles.At(i);
8d2cd130 692 Int_t ks = iparticle->GetStatusCode();
693 kf = CheckPDGCode(iparticle->GetPdgCode());
694// No initial state partons
695 if (ks==21) continue;
696//
697// Heavy Flavor Selection
698//
699 // quark ?
700 kf = TMath::Abs(kf);
701 Int_t kfl = kf;
9ff6c04c 702 // Resonance
f913ec4f 703
9ff6c04c 704 if (kfl > 100000) kfl %= 100000;
183a5ca9 705 if (kfl > 10000) kfl %= 10000;
8d2cd130 706 // meson ?
707 if (kfl > 10) kfl/=100;
708 // baryon
709 if (kfl > 10) kfl/=10;
8d2cd130 710 Int_t ipa = iparticle->GetFirstMother()-1;
711 Int_t kfMo = 0;
f913ec4f 712//
713// Establish mother daughter relation between heavy quarks and mesons
714//
715 if (kf >= fFlavorSelect && kf <= 6) {
716 Int_t idau = iparticle->GetFirstDaughter() - 1;
717 if (idau > -1) {
8507138f 718 TParticle* daughter = (TParticle *) fParticles.At(idau);
f913ec4f 719 Int_t pdgD = daughter->GetPdgCode();
720 if (pdgD == 91 || pdgD == 92) {
721 Int_t jmin = daughter->GetFirstDaughter() - 1;
722 Int_t jmax = daughter->GetLastDaughter() - 1;
2ab330c9 723 for (Int_t jp = jmin; jp <= jmax; jp++)
724 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
f913ec4f 725 } // is string or cluster
726 } // has daughter
727 } // heavy quark
8d2cd130 728
f913ec4f 729
8d2cd130 730 if (ipa > -1) {
8507138f 731 TParticle * mother = (TParticle *) fParticles.At(ipa);
8d2cd130 732 kfMo = TMath::Abs(mother->GetPdgCode());
733 }
f913ec4f 734
8d2cd130 735 // What to keep in Stack?
736 Bool_t flavorOK = kFALSE;
737 Bool_t selectOK = kFALSE;
738 if (fFeedDownOpt) {
32d6ef7d 739 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 740 } else {
32d6ef7d 741 if (kfl > fFlavorSelect) {
742 nc = -1;
743 break;
744 }
745 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 746 }
747 switch (fStackFillOpt) {
748 case kFlavorSelection:
32d6ef7d 749 selectOK = kTRUE;
750 break;
8d2cd130 751 case kParentSelection:
32d6ef7d 752 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
753 break;
8d2cd130 754 }
755 if (flavorOK && selectOK) {
756//
757// Heavy flavor hadron or quark
758//
759// Kinematic seletion on final state heavy flavor mesons
760 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
761 {
9ff6c04c 762 continue;
8d2cd130 763 }
764 pSelected[i] = 1;
765 if (ParentSelected(kf)) ++nParents; // Update parent count
766// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
767 } else {
768// Kinematic seletion on decay products
769 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 770 && !KinematicSelection(iparticle, 1))
8d2cd130 771 {
9ff6c04c 772 continue;
8d2cd130 773 }
774//
775// Decay products
776// Select if mother was selected and is not tracked
777
778 if (pSelected[ipa] &&
779 !trackIt[ipa] && // mother will be tracked ?
780 kfMo != 5 && // mother is b-quark, don't store fragments
781 kfMo != 4 && // mother is c-quark, don't store fragments
782 kf != 92) // don't store string
783 {
784//
785// Semi-stable or de-selected: diselect decay products:
786//
787//
788 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
789 {
790 Int_t ipF = iparticle->GetFirstDaughter();
791 Int_t ipL = iparticle->GetLastDaughter();
792 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
793 }
794// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
795 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
796 }
797 }
798 if (pSelected[i] == -1) pSelected[i] = 0;
799 if (!pSelected[i]) continue;
800 // Count quarks only if you did not include fragmentation
801 if (fFragmentation && kf <= 10) continue;
9ff6c04c 802
8d2cd130 803 nc++;
804// Decision on tracking
805 trackIt[i] = 0;
806//
807// Track final state particle
808 if (ks == 1) trackIt[i] = 1;
809// Track semi-stable particles
d25cfd65 810 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 811// Track particles selected by process if undecayed.
812 if (fForceDecay == kNoDecay) {
813 if (ParentSelected(kf)) trackIt[i] = 1;
814 } else {
815 if (ParentSelected(kf)) trackIt[i] = 0;
816 }
817 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
818//
819//
820
821 } // particle selection loop
822 if (nc > 0) {
823 for (i = 0; i<np; i++) {
824 if (!pSelected[i]) continue;
8507138f 825 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 826 kf = CheckPDGCode(iparticle->GetPdgCode());
827 Int_t ks = iparticle->GetStatusCode();
828 p[0] = iparticle->Px();
829 p[1] = iparticle->Py();
830 p[2] = iparticle->Pz();
a920faf9 831 p[3] = iparticle->Energy();
832
2590ccf9 833 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
834 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
835 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
836
8d2cd130 837 Float_t tof = kconv*iparticle->T();
838 Int_t ipa = iparticle->GetFirstMother()-1;
839 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 840
841 PushTrack(fTrackIt*trackIt[i], iparent, kf,
842 p[0], p[1], p[2], p[3],
843 origin[0], origin[1], origin[2], tof,
844 polar[0], polar[1], polar[2],
845 kPPrimary, nt, 1., ks);
8d2cd130 846 pParent[i] = nt;
dbd64db6 847 KeepTrack(nt);
848 fNprimaries++;
642f15cf 849 } // PushTrack loop
8d2cd130 850 }
851 } else {
852 nc = GenerateMB();
853 } // mb ?
f913ec4f 854
855 GetSubEventTime();
8d2cd130 856
234f6d32 857 delete[] pParent;
858 delete[] pSelected;
859 delete[] trackIt;
8d2cd130 860
861 if (nc > 0) {
862 switch (fCountMode) {
863 case kCountAll:
864 // printf(" Count all \n");
865 jev += nc;
866 break;
867 case kCountParents:
868 // printf(" Count parents \n");
869 jev += nParents;
870 break;
871 case kCountTrackables:
872 // printf(" Count trackable \n");
873 jev += nTkbles;
874 break;
875 }
876 if (jev >= fNpart || fNpart == -1) {
877 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 878
8d2cd130 879 fQ += fPythia->GetVINT(51);
880 fX1 += fPythia->GetVINT(41);
881 fX2 += fPythia->GetVINT(42);
882 fTrialsRun += fTrials;
883 fNev++;
884 MakeHeader();
885 break;
886 }
887 }
888 } // event loop
889 SetHighWaterMark(nt);
890// adjust weight due to kinematic selection
b88f5cea 891// AdjustWeights();
8d2cd130 892// get cross-section
893 fXsection=fPythia->GetPARI(1);
894}
895
896Int_t AliGenPythia::GenerateMB()
897{
898//
899// Min Bias selection and other global selections
900//
901 Int_t i, kf, nt, iparent;
902 Int_t nc = 0;
13cca24a 903 Double_t p[4];
904 Double_t polar[3] = {0,0,0};
905 Double_t origin[3] = {0,0,0};
8d2cd130 906// converts from mm/c to s
907 const Float_t kconv=0.001/2.999792458e8;
908
e0e89f40 909
910
8507138f 911 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
aea21c57 912
5fa4b20b 913
e0e89f40 914
8d2cd130 915 Int_t* pParent = new Int_t[np];
916 for (i=0; i< np; i++) pParent[i] = -1;
2f405d65 917 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8507138f 918 TParticle* jet1 = (TParticle *) fParticles.At(6);
919 TParticle* jet2 = (TParticle *) fParticles.At(7);
234f6d32 920 if (!CheckTrigger(jet1, jet2)) {
921 delete [] pParent;
922 return 0;
923 }
8d2cd130 924 }
e0e89f40 925
9fd8e520 926 // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
927 if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
ec2c406e 928
929 Bool_t ok = kFALSE;
930
931 Int_t pdg = 0;
9fd8e520 932 if (fFragPhotonInCalo) pdg = 22 ; // Photon
6d2ec66d 933 else if (fPi0InCalo) pdg = 111 ; // Pi0
ec2c406e 934
935 for (i=0; i< np; i++) {
8507138f 936 TParticle* iparticle = (TParticle *) fParticles.At(i);
ec2c406e 937 if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
9fd8e520 938 iparticle->Pt() > fFragPhotonOrPi0MinPt){
562cbbcf 939 Int_t imother = iparticle->GetFirstMother() - 1;
8507138f 940 TParticle* pmother = (TParticle *) fParticles.At(imother);
9fd8e520 941 if(pdg == 111 ||
82e0bdff 942 (pdg == 22 && pmother->GetStatusCode() != 11)) //No photon from hadron decay
9fd8e520 943 {
944 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
82e0bdff 945 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
9fd8e520 946 if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
947 (fCheckPHOS && IsInPHOS(phi,eta)) )
948 ok =kTRUE;
949 }
ec2c406e 950 }
951 }
1ee02229 952 if(!ok) {
953 delete[] pParent;
954 return 0;
955 }
ec2c406e 956 }
9dfe63b3 957
958 // Select beauty jets with electron in EMCAL
959 if (fProcess == kPyBeautyJets && fEleInEMCAL) {
960
961 Bool_t ok = kFALSE;
962
963 Int_t pdg = 11; //electron
964
70574ff8 965 Float_t pt = 0.;
966 Float_t eta = 0.;
967 Float_t phi = 0.;
9dfe63b3 968 for (i=0; i< np; i++) {
969 TParticle* iparticle = (TParticle *) fParticles.At(i);
970 if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
971 iparticle->Pt() > fElectronMinPt){
972 pt = iparticle->Pt();
973 phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
974 eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
975 if(IsInEMCAL(phi,eta))
976 ok =kTRUE;
977 }
978 }
92847124 979 if(!ok) {
980 delete[] pParent;
9dfe63b3 981 return 0;
92847124 982 }
983
9dfe63b3 984 AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
985 }
800be8b7 986 // Check for diffraction
987 if(fkTuneForDiff) {
988 if(fItune==320) {
989 if(!CheckDiffraction()) {
990 delete [] pParent;
991 return 0;
992 }
993 }
994 }
995
700b9416 996 // Check for minimum multiplicity
997 if (fTriggerMultiplicity > 0) {
998 Int_t multiplicity = 0;
999 for (i = 0; i < np; i++) {
8507138f 1000 TParticle * iparticle = (TParticle *) fParticles.At(i);
700b9416 1001
1002 Int_t statusCode = iparticle->GetStatusCode();
1003
1004 // Initial state particle
e296848e 1005 if (statusCode != 1)
700b9416 1006 continue;
38112f3f 1007 // eta cut
700b9416 1008 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1009 continue;
38112f3f 1010 // pt cut
1011 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1012 continue;
1013
700b9416 1014 TParticlePDG* pdgPart = iparticle->GetPDG();
1015 if (pdgPart && pdgPart->Charge() == 0)
1016 continue;
1017
1018 ++multiplicity;
1019 }
1020
1021 if (multiplicity < fTriggerMultiplicity) {
1022 delete [] pParent;
1023 return 0;
1024 }
38112f3f 1025 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
700b9416 1026 }
90a236ce 1027
1028 // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
1029 if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
1030
1031 Bool_t okd = kFALSE;
1032
1033 Int_t pdg = 22;
1034 Int_t iphcand = -1;
1035 for (i=0; i< np; i++) {
8507138f 1036 TParticle* iparticle = (TParticle *) fParticles.At(i);
90a236ce 1037 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1038 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
1039
1040 if(iparticle->GetStatusCode() == 1
1041 && iparticle->GetPdgCode() == pdg
1042 && iparticle->Pt() > fPhotonMinPt
04a26b0a 1043 && eta < fPHOSEta){
90a236ce 1044
1045 // first check if the photon is in PHOS phi
1046 if(IsInPHOS(phi,eta)){
1047 okd = kTRUE;
1048 break;
1049 }
1050 if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
1051
1052 }
1053 }
1054
1055 if(!okd && iphcand != -1) // execute rotation in phi
1056 RotatePhi(iphcand,okd);
1057
b3305127 1058 if(!okd) {
1059 delete [] pParent;
1060 return 0;
1061 }
90a236ce 1062 }
1063
7ce1321b 1064 if (fTriggerParticle) {
1065 Bool_t triggered = kFALSE;
1066 for (i = 0; i < np; i++) {
8507138f 1067 TParticle * iparticle = (TParticle *) fParticles.At(i);
7ce1321b 1068 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 1069 if (kf != fTriggerParticle) continue;
7ce1321b 1070 if (iparticle->Pt() == 0.) continue;
1071 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1072 triggered = kTRUE;
1073 break;
1074 }
234f6d32 1075 if (!triggered) {
1076 delete [] pParent;
1077 return 0;
1078 }
7ce1321b 1079 }
e0e89f40 1080
1081
1082 // Check if there is a ccbar or bbbar pair with at least one of the two
1083 // in fYMin < y < fYMax
2f405d65 1084
9dfe63b3 1085 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
9ccc3d0e 1086 TParticle *partCheck;
1087 TParticle *mother;
e0e89f40 1088 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 1089 Bool_t theChild=kFALSE;
5d76923e 1090 Bool_t triggered=kFALSE;
9ccc3d0e 1091 Float_t y;
1092 Int_t pdg,mpdg,mpdgUpperFamily;
e0e89f40 1093 for(i=0; i<np; i++) {
9ccc3d0e 1094 partCheck = (TParticle*)fParticles.At(i);
1095 pdg = partCheck->GetPdgCode();
1096 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1097 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1098 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1099 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1100 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1101 }
5d76923e 1102 if(fTriggerParticle) {
1103 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1104 }
9ccc3d0e 1105 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1106 Int_t mi = partCheck->GetFirstMother() - 1;
1107 if(mi<0) continue;
1108 mother = (TParticle*)fParticles.At(mi);
1109 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 1110 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 1111 if ( ParentSelected(mpdg) ||
1112 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1113 if (KinematicSelection(partCheck,1)) {
1114 theChild=kTRUE;
1115 }
1116 }
1117 }
e0e89f40 1118 }
9ccc3d0e 1119 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
234f6d32 1120 delete[] pParent;
e0e89f40 1121 return 0;
1122 }
9ccc3d0e 1123 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1124 delete[] pParent;
1125 return 0;
1126 }
5d76923e 1127 if(fTriggerParticle && !triggered) { // particle requested is not produced
1128 delete[] pParent;
1129 return 0;
1130 }
9ccc3d0e 1131
e0e89f40 1132 }
aea21c57 1133
0f6ee828 1134 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
f529e69b 1135 if ( (fProcess == kPyW ||
1136 fProcess == kPyZ ||
1137 fProcess == kPyMbDefault ||
1138 fProcess == kPyMb ||
6d2ec66d 1139 fProcess == kPyMbAtlasTuneMC09 ||
04081a91 1140 fProcess == kPyMbWithDirectPhoton ||
f529e69b 1141 fProcess == kPyMbNonDiffr)
0f6ee828 1142 && (fCutOnChild == 1) ) {
1143 if ( !CheckKinematicsOnChild() ) {
234f6d32 1144 delete[] pParent;
0f6ee828 1145 return 0;
1146 }
aea21c57 1147 }
1148
1149
f913ec4f 1150 for (i = 0; i < np; i++) {
8d2cd130 1151 Int_t trackIt = 0;
8507138f 1152 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 1153 kf = CheckPDGCode(iparticle->GetPdgCode());
1154 Int_t ks = iparticle->GetStatusCode();
1155 Int_t km = iparticle->GetFirstMother();
1156 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
1157 (ks != 1) ||
9dfe63b3 1158 ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
8d2cd130 1159 nc++;
1160 if (ks == 1) trackIt = 1;
1161 Int_t ipa = iparticle->GetFirstMother()-1;
1162
1163 iparent = (ipa > -1) ? pParent[ipa] : -1;
1164
1165//
1166// store track information
1167 p[0] = iparticle->Px();
1168 p[1] = iparticle->Py();
1169 p[2] = iparticle->Pz();
a920faf9 1170 p[3] = iparticle->Energy();
1406f599 1171
a920faf9 1172
2590ccf9 1173 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1174 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1175 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1176
f913ec4f 1177 Float_t tof = fEventTime + kconv * iparticle->T();
a920faf9 1178
1179 PushTrack(fTrackIt*trackIt, iparent, kf,
1180 p[0], p[1], p[2], p[3],
1181 origin[0], origin[1], origin[2], tof,
1182 polar[0], polar[1], polar[2],
1183 kPPrimary, nt, 1., ks);
dbd64db6 1184 fNprimaries++;
8d2cd130 1185 KeepTrack(nt);
1186 pParent[i] = nt;
f913ec4f 1187 SetHighWaterMark(nt);
1188
8d2cd130 1189 } // select particle
1190 } // particle loop
1191
234f6d32 1192 delete[] pParent;
e0e89f40 1193
f913ec4f 1194 return 1;
8d2cd130 1195}
1196
1197
1198void AliGenPythia::FinishRun()
1199{
1200// Print x-section summary
1201 fPythia->Pystat(1);
2779fc64 1202
1203 if (fNev > 0.) {
1204 fQ /= fNev;
1205 fX1 /= fNev;
1206 fX2 /= fNev;
1207 }
1208
8d2cd130 1209 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1210 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1211}
1212
7184e472 1213void AliGenPythia::AdjustWeights() const
8d2cd130 1214{
1215// Adjust the weights after generation of all events
1216//
e2bddf81 1217 if (gAlice) {
1218 TParticle *part;
1219 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1220 for (Int_t i=0; i<ntrack; i++) {
1221 part= gAlice->GetMCApp()->Particle(i);
1222 part->SetWeight(part->GetWeight()*fKineBias);
1223 }
8d2cd130 1224 }
1225}
1226
20e47f08 1227void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
8d2cd130 1228{
1229// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1230
1a626d4e 1231 fAProjectile = a1;
1232 fATarget = a2;
20e47f08 1233 fNucPdf = pdfset; // 0 EKS98 1 EPS08
1d568bc2 1234 fSetNuclei = kTRUE;
8d2cd130 1235}
1236
1237
1238void AliGenPythia::MakeHeader()
1239{
7184e472 1240//
1241// Make header for the simulated event
1242//
183a5ca9 1243 if (gAlice) {
1244 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1245 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1246 }
1247
8d2cd130 1248// Builds the event header, to be called after each event
e5c87a3d 1249 if (fHeader) delete fHeader;
1250 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1251//
1252// Event type
800be8b7 1253 if(fProcDiff>0){
1254 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1255 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1256 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1257 }
1258 else
e5c87a3d 1259 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1260//
1261// Number of trials
e5c87a3d 1262 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1263//
1264// Event Vertex
d25cfd65 1265 fHeader->SetPrimaryVertex(fVertex);
d07f0af2 1266 fHeader->SetInteractionTime(fEventTime);
dbd64db6 1267//
1268// Number of primaries
1269 fHeader->SetNProduced(fNprimaries);
8d2cd130 1270//
1271// Jets that have triggered
f913ec4f 1272
9dfe63b3 1273 //Need to store jets for b-jet studies too!
2f405d65 1274 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
8d2cd130 1275 {
1276 Int_t ntrig, njet;
1277 Float_t jets[4][10];
1278 GetJets(njet, ntrig, jets);
9ff6c04c 1279
8d2cd130 1280
1281 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1282 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1283 jets[3][i]);
1284 }
1285 }
5fa4b20b 1286//
1287// Copy relevant information from external header, if present.
1288//
1289 Float_t uqJet[4];
1290
1291 if (fRL) {
1292 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1293 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1294 {
1295 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1296
1297
1298 exHeader->TriggerJet(i, uqJet);
1299 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1300 }
1301 }
1302//
1303// Store quenching parameters
1304//
1305 if (fQuench){
b6262a45 1306 Double_t z[4] = {0.};
1307 Double_t xp = 0.;
1308 Double_t yp = 0.;
1309
7c21f297 1310 if (fQuench == 1) {
1311 // Pythia::Quench()
1312 fPythia->GetQuenchingParameters(xp, yp, z);
1bab4b79 1313 } else if (fQuench == 2){
7c21f297 1314 // Pyquen
1315 Double_t r1 = PARIMP.rb1;
1316 Double_t r2 = PARIMP.rb2;
1317 Double_t b = PARIMP.b1;
1318 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1319 Double_t phi = PARIMP.psib1;
1320 xp = r * TMath::Cos(phi);
1321 yp = r * TMath::Sin(phi);
1322
1bab4b79 1323 } else if (fQuench == 4) {
1324 // QPythia
5831e75f 1325 Double_t xy[2];
e9719084 1326 Double_t i0i1[2];
1327 AliFastGlauber::Instance()->GetSavedXY(xy);
1328 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
5831e75f 1329 xp = xy[0];
1330 yp = xy[1];
e6fe9b82 1331 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
7c21f297 1332 }
1bab4b79 1333
7c21f297 1334 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1335 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1bab4b79 1336 }
beac474c 1337//
1338// Store pt^hard
1339 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 1340//
cf57b268 1341// Pass header
5fa4b20b 1342//
cf57b268 1343 AddHeader(fHeader);
4c4eac97 1344 fHeader = 0x0;
8d2cd130 1345}
cf57b268 1346
c2fc9d31 1347Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
8d2cd130 1348{
1349// Check the kinematic trigger condition
1350//
1351 Double_t eta[2];
1352 eta[0] = jet1->Eta();
1353 eta[1] = jet2->Eta();
1354 Double_t phi[2];
1355 phi[0] = jet1->Phi();
1356 phi[1] = jet2->Phi();
1357 Int_t pdg[2];
1358 pdg[0] = jet1->GetPdgCode();
1359 pdg[1] = jet2->GetPdgCode();
1360 Bool_t triggered = kFALSE;
1361
2f405d65 1362 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
8d2cd130 1363 Int_t njets = 0;
1364 Int_t ntrig = 0;
1365 Float_t jets[4][10];
1366//
1367// Use Pythia clustering on parton level to determine jet axis
1368//
1369 GetJets(njets, ntrig, jets);
1370
76d6ba9a 1371 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1372//
1373 } else {
1374 Int_t ij = 0;
1375 Int_t ig = 1;
1376 if (pdg[0] == kGamma) {
1377 ij = 1;
1378 ig = 0;
1379 }
1380 //Check eta range first...
1381 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1382 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1383 {
1384 //Eta is okay, now check phi range
1385 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1386 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1387 {
1388 triggered = kTRUE;
1389 }
1390 }
1391 }
1392 return triggered;
1393}
aea21c57 1394
1395
aea21c57 1396
7184e472 1397Bool_t AliGenPythia::CheckKinematicsOnChild(){
1398//
1399//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1400//
aea21c57 1401 Bool_t checking = kFALSE;
1402 Int_t j, kcode, ks, km;
1403 Int_t nPartAcc = 0; //number of particles in the acceptance range
1404 Int_t numberOfAcceptedParticles = 1;
1405 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1406 Int_t npart = fParticles.GetEntriesFast();
aea21c57 1407
0f6ee828 1408 for (j = 0; j<npart; j++) {
8507138f 1409 TParticle * jparticle = (TParticle *) fParticles.At(j);
aea21c57 1410 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1411 ks = jparticle->GetStatusCode();
1412 km = jparticle->GetFirstMother();
1413
1414 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1415 nPartAcc++;
1416 }
0f6ee828 1417 if( numberOfAcceptedParticles <= nPartAcc){
1418 checking = kTRUE;
1419 break;
1420 }
aea21c57 1421 }
0f6ee828 1422
aea21c57 1423 return checking;
aea21c57 1424}
1425
5fa4b20b 1426void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1427{
1058d9df 1428 //
1429 // Load event into Pythia Common Block
1430 //
1431
1432 Int_t npart = stack -> GetNprimary();
1433 Int_t n0 = 0;
1434
1435 if (!flag) {
1436 (fPythia->GetPyjets())->N = npart;
1437 } else {
1438 n0 = (fPythia->GetPyjets())->N;
1439 (fPythia->GetPyjets())->N = n0 + npart;
1440 }
1441
1442
1443 for (Int_t part = 0; part < npart; part++) {
1444 TParticle *mPart = stack->Particle(part);
32d6ef7d 1445
1058d9df 1446 Int_t kf = mPart->GetPdgCode();
1447 Int_t ks = mPart->GetStatusCode();
1448 Int_t idf = mPart->GetFirstDaughter();
1449 Int_t idl = mPart->GetLastDaughter();
1450
1451 if (reHadr) {
1452 if (ks == 11 || ks == 12) {
1453 ks -= 10;
1454 idf = -1;
1455 idl = -1;
1456 }
32d6ef7d 1457 }
1458
1058d9df 1459 Float_t px = mPart->Px();
1460 Float_t py = mPart->Py();
1461 Float_t pz = mPart->Pz();
1462 Float_t e = mPart->Energy();
1463 Float_t m = mPart->GetCalcMass();
32d6ef7d 1464
1058d9df 1465
1466 (fPythia->GetPyjets())->P[0][part+n0] = px;
1467 (fPythia->GetPyjets())->P[1][part+n0] = py;
1468 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1469 (fPythia->GetPyjets())->P[3][part+n0] = e;
1470 (fPythia->GetPyjets())->P[4][part+n0] = m;
1471
1472 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1473 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1474 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1475 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1476 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1477 }
1478}
1479
c2fc9d31 1480void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1058d9df 1481{
1482 //
1483 // Load event into Pythia Common Block
1484 //
1485
1486 Int_t npart = stack -> GetEntries();
1487 Int_t n0 = 0;
1488
1489 if (!flag) {
1490 (fPythia->GetPyjets())->N = npart;
1491 } else {
1492 n0 = (fPythia->GetPyjets())->N;
1493 (fPythia->GetPyjets())->N = n0 + npart;
1494 }
1495
1496
1497 for (Int_t part = 0; part < npart; part++) {
1498 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
92847124 1499 if (!mPart) continue;
1500
1058d9df 1501 Int_t kf = mPart->GetPdgCode();
1502 Int_t ks = mPart->GetStatusCode();
1503 Int_t idf = mPart->GetFirstDaughter();
1504 Int_t idl = mPart->GetLastDaughter();
1505
1506 if (reHadr) {
92847124 1507 if (ks == 11 || ks == 12) {
1508 ks -= 10;
1509 idf = -1;
1510 idl = -1;
1511 }
8d2cd130 1512 }
1058d9df 1513
1514 Float_t px = mPart->Px();
1515 Float_t py = mPart->Py();
1516 Float_t pz = mPart->Pz();
1517 Float_t e = mPart->Energy();
1518 Float_t m = mPart->GetCalcMass();
1519
1520
1521 (fPythia->GetPyjets())->P[0][part+n0] = px;
1522 (fPythia->GetPyjets())->P[1][part+n0] = py;
1523 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1524 (fPythia->GetPyjets())->P[3][part+n0] = e;
1525 (fPythia->GetPyjets())->P[4][part+n0] = m;
1526
1527 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1528 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1529 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1530 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1531 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1532 }
8d2cd130 1533}
1534
5fa4b20b 1535
014a9521 1536void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1537{
1538//
1539// Calls the Pythia jet finding algorithm to find jets in the current event
1540//
1541//
8d2cd130 1542//
1543// Save jets
1544 Int_t n = fPythia->GetN();
1545
1546//
1547// Run Jet Finder
1548 fPythia->Pycell(njets);
1549 Int_t i;
1550 for (i = 0; i < njets; i++) {
1551 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1552 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1553 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1554 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1555
1556 jets[0][i] = px;
1557 jets[1][i] = py;
1558 jets[2][i] = pz;
1559 jets[3][i] = e;
1560 }
1561}
1562
1563
1564
1565void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1566{
1567//
1568// Calls the Pythia clustering algorithm to find jets in the current event
1569//
1570 Int_t n = fPythia->GetN();
1571 nJets = 0;
1572 nJetsTrig = 0;
1573 if (fJetReconstruction == kCluster) {
1574//
1575// Configure cluster algorithm
1576//
1577 fPythia->SetPARU(43, 2.);
1578 fPythia->SetMSTU(41, 1);
1579//
1580// Call cluster algorithm
1581//
1582 fPythia->Pyclus(nJets);
1583//
1584// Loading jets from common block
1585//
1586 } else {
592f8307 1587
8d2cd130 1588//
1589// Run Jet Finder
1590 fPythia->Pycell(nJets);
1591 }
1592
1593 Int_t i;
1594 for (i = 0; i < nJets; i++) {
1595 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1596 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1597 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1598 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1599 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1600 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1601 Float_t theta = TMath::ATan2(pt,pz);
1602 Float_t et = e * TMath::Sin(theta);
1603 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1604 if (
1605 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1606 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1607 et > fEtMinJet && et < fEtMaxJet
1608 )
1609 {
1610 jets[0][nJetsTrig] = px;
1611 jets[1][nJetsTrig] = py;
1612 jets[2][nJetsTrig] = pz;
1613 jets[3][nJetsTrig] = e;
1614 nJetsTrig++;
5fa4b20b 1615// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1616 } else {
1617// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1618 }
1619 }
1620}
1621
f913ec4f 1622void AliGenPythia::GetSubEventTime()
1623{
1624 // Calculates time of the next subevent
9d7108a7 1625 fEventTime = 0.;
1626 if (fEventsTime) {
1627 TArrayF &array = *fEventsTime;
1628 fEventTime = array[fCurSubEvent++];
1629 }
1630 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1631 return;
f913ec4f 1632}
8d2cd130 1633
c2fc9d31 1634Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
ec2c406e 1635{
1636 // Is particle in EMCAL acceptance?
ec2c406e 1637 // phi in degrees, etamin=-etamax
9fd8e520 1638 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1639 eta < fEMCALEta )
ec2c406e 1640 return kTRUE;
1641 else
1642 return kFALSE;
1643}
1644
c2fc9d31 1645Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta) const
ec2c406e 1646{
1647 // Is particle in PHOS acceptance?
1648 // Acceptance slightly larger considered.
1649 // phi in degrees, etamin=-etamax
9fd8e520 1650 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1651 eta < fPHOSEta )
ec2c406e 1652 return kTRUE;
1653 else
1654 return kFALSE;
1655}
1656
90a236ce 1657void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1658{
1659 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1660 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1661 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
800be8b7 1662 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
90a236ce 1663
1664 //calculate deltaphi
8507138f 1665 TParticle* ph = (TParticle *) fParticles.At(iphcand);
90a236ce 1666 Double_t phphi = ph->Phi();
1667 Double_t deltaphi = phiPHOS - phphi;
1668
1669
1670
1671 //loop for all particles and produce the phi rotation
8507138f 1672 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
90a236ce 1673 Double_t oldphi, newphi;
777e69b0 1674 Double_t newVx, newVy, r, vZ, time;
1675 Double_t newPx, newPy, pt, pz, e;
90a236ce 1676 for(Int_t i=0; i< np; i++) {
8507138f 1677 TParticle* iparticle = (TParticle *) fParticles.At(i);
90a236ce 1678 oldphi = iparticle->Phi();
1679 newphi = oldphi + deltaphi;
1680 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1681 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1682
777e69b0 1683 r = iparticle->R();
1684 newVx = r * TMath::Cos(newphi);
1685 newVy = r * TMath::Sin(newphi);
1686 vZ = iparticle->Vz(); // don't transform
90a236ce 1687 time = iparticle->T(); // don't transform
1688
1689 pt = iparticle->Pt();
777e69b0 1690 newPx = pt * TMath::Cos(newphi);
1691 newPy = pt * TMath::Sin(newphi);
1692 pz = iparticle->Pz(); // don't transform
1693 e = iparticle->Energy(); // don't transform
90a236ce 1694
1695 // apply rotation
777e69b0 1696 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1697 iparticle->SetMomentum(newPx, newPy, pz, e);
90a236ce 1698
1699 } //end particle loop
1700
1701 // now let's check that we put correctly the candidate photon in PHOS
1702 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1703 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1704 if(IsInPHOS(phi,eta))
1705 okdd = kTRUE;
1706}
ec2c406e 1707
1708
589380c6 1709
800be8b7 1710Bool_t AliGenPythia::CheckDiffraction()
1711{
1712 // use this method only with Perugia-0 tune!
1713
1714 // printf("AAA\n");
1715
1716 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1717
1718 Int_t iPart1=-1;
1719 Int_t iPart2=-1;
1720
1721 Double_t y1 = 1e10;
1722 Double_t y2 = -1e10;
1723
1724 const Int_t kNstable=20;
1725 const Int_t pdgStable[20] = {
1726 22, // Photon
1727 11, // Electron
1728 12, // Electron Neutrino
1729 13, // Muon
1730 14, // Muon Neutrino
1731 15, // Tau
1732 16, // Tau Neutrino
1733 211, // Pion
1734 321, // Kaon
1735 311, // K0
1736 130, // K0s
1737 310, // K0l
1738 2212, // Proton
1739 2112, // Neutron
1740 3122, // Lambda_0
1741 3112, // Sigma Minus
1742 3222, // Sigma Plus
1743 3312, // Xsi Minus
1744 3322, // Xsi0
1745 3334 // Omega
1746 };
1747
1748 for (Int_t i = 0; i < np; i++) {
1749 TParticle * part = (TParticle *) fParticles.At(i);
1750
1751 Int_t statusCode = part->GetStatusCode();
1752
1753 // Initial state particle
1754 if (statusCode != 1)
1755 continue;
1756
1757 Int_t pdg = TMath::Abs(part->GetPdgCode());
1758 Bool_t isStable = kFALSE;
1759 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1760 if (pdg == pdgStable[i1]) {
1761 isStable = kTRUE;
1762 break;
1763 }
1764 }
1765 if(!isStable)
1766 continue;
1767
1768 Double_t y = part->Y();
1769
1770 if (y < y1)
1771 {
1772 y1 = y;
1773 iPart1 = i;
1774 }
1775 if (y > y2)
1776 {
1777 y2 = y;
1778 iPart2 = i;
1779 }
1780 }
1781
1782 if(iPart1<0 || iPart2<0) return kFALSE;
1783
1784 y1=TMath::Abs(y1);
1785 y2=TMath::Abs(y2);
1786
1787 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1788 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1789
1790 Int_t pdg1 = part1->GetPdgCode();
1791 Int_t pdg2 = part2->GetPdgCode();
1792
1793
1794 Int_t iPart = -1;
1795 if (pdg1 == 2212 && pdg2 == 2212)
1796 {
1797 if(y1 > y2)
1798 iPart = iPart1;
1799 else if(y1 < y2)
1800 iPart = iPart2;
1801 else {
1802 iPart = iPart1;
1803 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1804 }
1805 }
1806 else if (pdg1 == 2212)
1807 iPart = iPart1;
1808 else if (pdg2 == 2212)
1809 iPart = iPart2;
1810
1811
1812
1813
1814
1815 Double_t M=-1.;
1816 if(iPart>0) {
1817 TParticle * part = (TParticle *) fParticles.At(iPart);
1818 Double_t E= part->Energy();
1819 Double_t P= part->P();
1820 M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
1821 }
1822
1823Int_t nbin=120;
1824Double_t bin[]={
18251.080000, 1.274258, 1.468516, 1.662773, 1.857031, 2.051289,
18262.245547, 2.439805, 2.634062, 2.828320, 3.022578, 3.216836,
18273.411094, 3.605352, 3.799609, 3.993867, 4.188125, 4.382383,
18284.576641, 4.770898, 4.965156, 5.547930, 6.130703, 6.713477,
18297.296250, 7.879023, 8.461797, 9.044570, 9.627344, 10.210117,
183010.792891, 11.375664, 11.958437, 12.541211, 13.123984, 13.706758,
183114.289531, 14.872305, 15.455078, 16.037852, 16.620625, 17.203398,
183217.786172, 18.368945, 18.951719, 19.534492, 20.117266, 20.700039,
183321.282812, 21.865586, 22.448359, 23.031133, 23.613906, 24.196680,
183424.779453, 25.362227, 25.945000, 26.527773, 27.110547, 27.693320,
183528.276094, 28.858867, 29.441641, 30.024414, 30.607187, 31.189961,
183631.772734, 32.355508, 32.938281, 33.521055, 34.103828, 34.686602,
183735.269375, 35.852148, 36.434922, 37.017695, 37.600469, 38.183242,
183838.766016, 39.348789, 39.931562, 40.514336, 41.097109, 41.679883,
183942.262656, 42.845430, 43.428203, 44.010977, 44.593750, 45.176523,
184045.759297, 46.342070, 46.924844, 47.507617, 48.090391, 48.673164,
184149.255937, 49.838711, 50.421484, 57.220508, 64.019531, 70.818555,
184277.617578, 84.416602, 91.215625, 98.014648, 104.813672, 111.612695,
1843118.411719, 125.210742, 132.009766, 138.808789, 145.607812, 152.406836,
1844159.205859, 166.004883, 172.803906, 179.602930, 186.401953, 193.200977,
1845200.000000};
1846Double_t w[]={
18471.000000, 0.275812, 0.267308, 0.263268, 0.263828, 0.263039,
18480.260546, 0.259344, 0.255477, 0.253782, 0.253562, 0.252492,
18490.251076, 0.247862, 0.248933, 0.243599, 0.244255, 0.238506,
18500.239546, 0.237845, 0.237977, 0.229379, 0.224771, 0.217581,
18510.208860, 0.207241, 0.198955, 0.196449, 0.193827, 0.190220,
18520.184715, 0.183067, 0.178325, 0.175887, 0.171299, 0.168718,
18530.167453, 0.165153, 0.163457, 0.159637, 0.156855, 0.155488,
18540.154545, 0.155968, 0.150488, 0.148797, 0.145358, 0.146196,
18550.145891, 0.142752, 0.145072, 0.141265, 0.141857, 0.138462,
18560.142992, 0.141357, 0.139391, 0.139629, 0.135197, 0.135731,
18570.133194, 0.137190, 0.135745, 0.134522, 0.136094, 0.130405,
18580.128371, 0.131680, 0.130591, 0.133539, 0.129370, 0.128254,
18590.128262, 0.131088, 0.128294, 0.130070, 0.124553, 0.131489,
18600.128038, 0.122997, 0.130699, 0.125630, 0.124746, 0.123679,
18610.127864, 0.125776, 0.126272, 0.121492, 0.124929, 0.122221,
18620.127017, 0.123706, 0.122701, 0.123731, 0.122219, 0.120783,
18630.120324, 0.120228, 0.123615, 0.120589, 0.119549, 0.118836,
18640.118146, 0.120175, 0.122031, 0.122076, 0.122849, 0.122627,
18650.126281, 0.127696, 0.129084, 0.130001, 0.134062, 0.134786,
18660.137286, 0.136625, 0.139091, 0.143692, 0.144781, 0.146768};
1867
1868 Double_t wSD=1.;
1869 Double_t wDD=0.178783;
1870 //Double_t wND=0.220200;
1871 Double_t wND=0.220200+0.001;
1872
1873 if(M>-1 && M<bin[0]) return kFALSE;
1874 if(M>bin[nbin]) M=-1;
1875
1876 Int_t procType=fPythia->GetMSTI(1);
1877 Int_t proc0=2;
1878 if(procType== 94) proc0=1;
1879 if(procType== 92 || procType== 93) proc0=0;
1880
1881
1882 // printf("M = %f bin[nbin] = %f\n",M, bin[nbin]);
1883
1884 Int_t proc=2;
1885 if(M>0) proc=0;
1886 else if(proc0==1) proc=1;
1887
1888 if(proc==0 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wSD) return kFALSE;
1889 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1890 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
1891
1892
1893 // if(proc==1 || proc==2) return kFALSE;
1894
1895 if(proc!=0) {
1896 if(proc0!=0) fProcDiff = procType;
1897 else fProcDiff = 95;
1898 return kTRUE;
1899 }
1900
1901 Int_t ibin=-1;
1902 for(Int_t i=0; i<nbin; i++)
1903 if(M>bin[i] && M<=bin[i+1]) {
1904 ibin=i;
1905 // printf("Mi> %f && Mi< %f\n", bin[i], bin[i+1]);
1906 break;
1907 }
1908
1909 // printf("w[ibin] = %f\n", w[ibin]);
1910
1911 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> w[ibin]) return kFALSE;
1912
1913 // printf("iPart = %d\n", iPart);
1914
1915 if(iPart==iPart1) fProcDiff=93;
1916 else if(iPart==iPart2) fProcDiff=92;
1917 else {
1918 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1919
1920 }
1921
1922 return kTRUE;
1923}