]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA6/AliGenPythia.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Macro to calculate the resolution and the efficiency of chamber(s) (Nicolas)
[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
28#include <TDatabasePDG.h>
29#include <TParticle.h>
30#include <TPDGCode.h>
31#include <TSystem.h>
32#include <TTree.h>
8d2cd130 33#include "AliConst.h"
34#include "AliDecayerPythia.h"
35#include "AliGenPythia.h"
5fa4b20b 36#include "AliHeader.h"
8d2cd130 37#include "AliGenPythiaEventHeader.h"
38#include "AliPythia.h"
7cdba479 39#include "AliPythiaRndm.h"
8d2cd130 40#include "AliRun.h"
7ea3ea5b 41#include "AliStack.h"
42#include "AliRunLoader.h"
5d12ce38 43#include "AliMC.h"
7c21f297 44#include "pyquenCommon.h"
8d2cd130 45
014a9521 46ClassImp(AliGenPythia)
8d2cd130 47
e8a8adcd 48
49AliGenPythia::AliGenPythia():
50 AliGenMC(),
51 fProcess(kPyCharm),
52 fStrucFunc(kCTEQ5L),
53 fEnergyCMS(5500.),
54 fKineBias(0.),
55 fTrials(0),
56 fTrialsRun(0),
57 fQ(0.),
58 fX1(0.),
59 fX2(0.),
60 fEventTime(0.),
61 fInteractionRate(0.),
62 fTimeWindow(0.),
63 fCurSubEvent(0),
64 fEventsTime(0),
65 fNev(0),
66 fFlavorSelect(0),
67 fXsection(0.),
68 fPythia(0),
69 fPtHardMin(0.),
70 fPtHardMax(1.e4),
71 fYHardMin(-1.e10),
72 fYHardMax(1.e10),
73 fGinit(1),
74 fGfinal(1),
75 fHadronisation(1),
76 fNpartons(0),
77 fReadFromFile(0),
78 fQuench(0),
79 fPtKick(1.),
80 fFullEvent(kTRUE),
81 fDecayer(new AliDecayerPythia()),
82 fDebugEventFirst(-1),
83 fDebugEventLast(-1),
84 fEtMinJet(0.),
85 fEtMaxJet(1.e4),
86 fEtaMinJet(-20.),
87 fEtaMaxJet(20.),
88 fPhiMinJet(0.),
89 fPhiMaxJet(2.* TMath::Pi()),
90 fJetReconstruction(kCell),
91 fEtaMinGamma(-20.),
92 fEtaMaxGamma(20.),
93 fPhiMinGamma(0.),
94 fPhiMaxGamma(2. * TMath::Pi()),
95 fPycellEtaMax(2.),
96 fPycellNEta(274),
97 fPycellNPhi(432),
98 fPycellThreshold(0.),
99 fPycellEtSeed(4.),
100 fPycellMinEtJet(10.),
101 fPycellMaxRadius(1.),
102 fStackFillOpt(kFlavorSelection),
103 fFeedDownOpt(kTRUE),
104 fFragmentation(kTRUE),
105 fSetNuclei(kFALSE),
106 fNewMIS(kFALSE),
107 fHFoff(kFALSE),
108 fTriggerParticle(0),
109 fTriggerEta(0.9),
110 fCountMode(kCountAll),
111 fHeader(0),
112 fRL(0),
ec2c406e 113 fFileName(0),
9fd8e520 114 fFragPhotonInCalo(kFALSE),
ec2c406e 115 fPi0InCalo(kFALSE) ,
9fd8e520 116 fCheckEMCAL(kFALSE),
117 fCheckPHOS(kFALSE),
118 fFragPhotonOrPi0MinPt(0),
119 fPHOSMinPhi(219.),
120 fPHOSMaxPhi(321.),
121 fPHOSEta(0.13),
122 fEMCALMinPhi(79.),
123 fEMCALMaxPhi(191.),
124 fEMCALEta(0.71)
ec2c406e 125
8d2cd130 126{
127// Default Constructor
2d9f9817 128 SetNuclei(0,0);
7cdba479 129 if (!AliPythiaRndm::GetPythiaRandom())
e8a8adcd 130 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 131}
132
133AliGenPythia::AliGenPythia(Int_t npart)
e8a8adcd 134 :AliGenMC(npart),
135 fProcess(kPyCharm),
136 fStrucFunc(kCTEQ5L),
137 fEnergyCMS(5500.),
138 fKineBias(0.),
139 fTrials(0),
140 fTrialsRun(0),
141 fQ(0.),
142 fX1(0.),
143 fX2(0.),
144 fEventTime(0.),
145 fInteractionRate(0.),
146 fTimeWindow(0.),
147 fCurSubEvent(0),
148 fEventsTime(0),
149 fNev(0),
150 fFlavorSelect(0),
151 fXsection(0.),
152 fPythia(0),
153 fPtHardMin(0.),
154 fPtHardMax(1.e4),
155 fYHardMin(-1.e10),
156 fYHardMax(1.e10),
157 fGinit(kTRUE),
158 fGfinal(kTRUE),
159 fHadronisation(kTRUE),
160 fNpartons(0),
161 fReadFromFile(kFALSE),
162 fQuench(kFALSE),
163 fPtKick(1.),
164 fFullEvent(kTRUE),
165 fDecayer(new AliDecayerPythia()),
166 fDebugEventFirst(-1),
167 fDebugEventLast(-1),
168 fEtMinJet(0.),
169 fEtMaxJet(1.e4),
170 fEtaMinJet(-20.),
171 fEtaMaxJet(20.),
172 fPhiMinJet(0.),
173 fPhiMaxJet(2.* TMath::Pi()),
174 fJetReconstruction(kCell),
175 fEtaMinGamma(-20.),
176 fEtaMaxGamma(20.),
177 fPhiMinGamma(0.),
178 fPhiMaxGamma(2. * TMath::Pi()),
179 fPycellEtaMax(2.),
180 fPycellNEta(274),
181 fPycellNPhi(432),
182 fPycellThreshold(0.),
183 fPycellEtSeed(4.),
184 fPycellMinEtJet(10.),
185 fPycellMaxRadius(1.),
186 fStackFillOpt(kFlavorSelection),
187 fFeedDownOpt(kTRUE),
188 fFragmentation(kTRUE),
189 fSetNuclei(kFALSE),
190 fNewMIS(kFALSE),
191 fHFoff(kFALSE),
192 fTriggerParticle(0),
193 fTriggerEta(0.9),
194 fCountMode(kCountAll),
195 fHeader(0),
196 fRL(0),
ec2c406e 197 fFileName(0),
9fd8e520 198 fFragPhotonInCalo(kFALSE),
ec2c406e 199 fPi0InCalo(kFALSE) ,
9fd8e520 200 fCheckEMCAL(kFALSE),
201 fCheckPHOS(kFALSE),
202 fFragPhotonOrPi0MinPt(0),
203 fPHOSMinPhi(219.),
204 fPHOSMaxPhi(321.),
205 fPHOSEta(0.13),
206 fEMCALMinPhi(79.),
207 fEMCALMaxPhi(191.),
208 fEMCALEta(0.71)
8d2cd130 209{
210// default charm production at 5. 5 TeV
211// semimuonic decay
212// structure function GRVHO
213//
214 fName = "Pythia";
215 fTitle= "Particle Generator using PYTHIA";
8d2cd130 216 SetForceDecay();
8d2cd130 217 // Set random number generator
7cdba479 218 if (!AliPythiaRndm::GetPythiaRandom())
219 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 220 fParticles = new TClonesArray("TParticle",1000);
2d9f9817 221 SetNuclei(0,0);
76d6ba9a 222 }
8d2cd130 223
224AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
e8a8adcd 225 :AliGenMC(Pythia),
226 fProcess(kPyCharm),
227 fStrucFunc(kCTEQ5L),
228 fEnergyCMS(5500.),
229 fKineBias(0.),
230 fTrials(0),
231 fTrialsRun(0),
232 fQ(0.),
233 fX1(0.),
234 fX2(0.),
235 fEventTime(0.),
236 fInteractionRate(0.),
237 fTimeWindow(0.),
238 fCurSubEvent(0),
239 fEventsTime(0),
240 fNev(0),
241 fFlavorSelect(0),
242 fXsection(0.),
243 fPythia(0),
244 fPtHardMin(0.),
245 fPtHardMax(1.e4),
246 fYHardMin(-1.e10),
247 fYHardMax(1.e10),
248 fGinit(kTRUE),
249 fGfinal(kTRUE),
250 fHadronisation(kTRUE),
251 fNpartons(0),
252 fReadFromFile(kFALSE),
253 fQuench(kFALSE),
254 fPtKick(1.),
255 fFullEvent(kTRUE),
256 fDecayer(new AliDecayerPythia()),
257 fDebugEventFirst(-1),
258 fDebugEventLast(-1),
259 fEtMinJet(0.),
260 fEtMaxJet(1.e4),
261 fEtaMinJet(-20.),
262 fEtaMaxJet(20.),
263 fPhiMinJet(0.),
264 fPhiMaxJet(2.* TMath::Pi()),
265 fJetReconstruction(kCell),
266 fEtaMinGamma(-20.),
267 fEtaMaxGamma(20.),
268 fPhiMinGamma(0.),
269 fPhiMaxGamma(2. * TMath::Pi()),
270 fPycellEtaMax(2.),
271 fPycellNEta(274),
272 fPycellNPhi(432),
273 fPycellThreshold(0.),
274 fPycellEtSeed(4.),
275 fPycellMinEtJet(10.),
276 fPycellMaxRadius(1.),
277 fStackFillOpt(kFlavorSelection),
278 fFeedDownOpt(kTRUE),
279 fFragmentation(kTRUE),
280 fSetNuclei(kFALSE),
281 fNewMIS(kFALSE),
282 fHFoff(kFALSE),
283 fTriggerParticle(0),
284 fTriggerEta(0.9),
285 fCountMode(kCountAll),
286 fHeader(0),
287 fRL(0),
ec2c406e 288 fFileName(0),
9fd8e520 289 fFragPhotonInCalo(kFALSE),
ec2c406e 290 fPi0InCalo(kFALSE) ,
9fd8e520 291 fCheckEMCAL(kFALSE),
292 fCheckPHOS(kFALSE),
293 fFragPhotonOrPi0MinPt(0),
294 fPHOSMinPhi(219.),
295 fPHOSMaxPhi(321.),
296 fPHOSEta(0.13),
297 fEMCALMinPhi(79.),
298 fEMCALMaxPhi(191.),
299 fEMCALEta(0.71)
8d2cd130 300{
301// copy constructor
302 Pythia.Copy(*this);
303}
304
305AliGenPythia::~AliGenPythia()
306{
307// Destructor
9d7108a7 308 if(fEventsTime) delete fEventsTime;
309}
310
311void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
312{
313// Generate pileup using user specified rate
314 fInteractionRate = rate;
315 fTimeWindow = timewindow;
316 GeneratePileup();
317}
318
319void AliGenPythia::GeneratePileup()
320{
321// Generate sub events time for pileup
322 fEventsTime = 0;
323 if(fInteractionRate == 0.) {
324 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
325 return;
326 }
327
328 Int_t npart = NumberParticles();
329 if(npart < 0) {
330 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
331 return;
332 }
333
334 if(fEventsTime) delete fEventsTime;
335 fEventsTime = new TArrayF(npart);
336 TArrayF &array = *fEventsTime;
337 for(Int_t ipart = 0; ipart < npart; ipart++)
338 array[ipart] = 0.;
339
340 Float_t eventtime = 0.;
341 while(1)
342 {
343 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
344 if(eventtime > fTimeWindow) break;
345 array.Set(array.GetSize()+1);
346 array[array.GetSize()-1] = eventtime;
347 }
348
349 eventtime = 0.;
350 while(1)
351 {
352 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
353 if(TMath::Abs(eventtime) > fTimeWindow) break;
354 array.Set(array.GetSize()+1);
355 array[array.GetSize()-1] = eventtime;
356 }
357
358 SetNumberParticles(fEventsTime->GetSize());
8d2cd130 359}
360
592f8307 361void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
362 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
363{
364// Set pycell parameters
365 fPycellEtaMax = etamax;
366 fPycellNEta = neta;
367 fPycellNPhi = nphi;
368 fPycellThreshold = thresh;
369 fPycellEtSeed = etseed;
370 fPycellMinEtJet = minet;
371 fPycellMaxRadius = r;
372}
373
374
375
8d2cd130 376void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
377{
378 // Set a range of event numbers, for which a table
379 // of generated particle will be printed
380 fDebugEventFirst = eventFirst;
381 fDebugEventLast = eventLast;
382 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
383}
384
385void AliGenPythia::Init()
386{
387// Initialisation
388
389 SetMC(AliPythia::Instance());
b88f5cea 390 fPythia=(AliPythia*) fMCEvGen;
e2bddf81 391
8d2cd130 392//
393 fParentWeight=1./Float_t(fNpart);
394//
395// Forward Paramters to the AliPythia object
396 fDecayer->SetForceDecay(fForceDecay);
397 fDecayer->Init();
398
399
400 fPythia->SetCKIN(3,fPtHardMin);
401 fPythia->SetCKIN(4,fPtHardMax);
402 fPythia->SetCKIN(7,fYHardMin);
403 fPythia->SetCKIN(8,fYHardMax);
404
1a626d4e 405 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
8d2cd130 406 // Fragmentation?
407 if (fFragmentation) {
408 fPythia->SetMSTP(111,1);
409 } else {
410 fPythia->SetMSTP(111,0);
411 }
412
413
414// initial state radiation
415 fPythia->SetMSTP(61,fGinit);
416// final state radiation
417 fPythia->SetMSTP(71,fGfinal);
418// pt - kick
419 if (fPtKick > 0.) {
420 fPythia->SetMSTP(91,1);
421 fPythia->SetPARP(91,fPtKick);
422 } else {
423 fPythia->SetMSTP(91,0);
424 }
425
5fa4b20b 426
427 if (fReadFromFile) {
428 fRL = AliRunLoader::Open(fFileName, "Partons");
429 fRL->LoadKinematics();
430 fRL->LoadHeader();
431 } else {
432 fRL = 0x0;
433 }
beac474c 434// Switch off Heavy Flavors on request
435 if (fHFoff) {
436 fPythia->SetMSTP(58, 3);
437 fPythia->SetMSTJ(45, 3);
438 for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
439 }
8d2cd130 440 //
441 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
442
443// Parent and Children Selection
444 switch (fProcess)
445 {
65f2626c 446 case kPyOldUEQ2ordered:
447 case kPyOldUEQ2ordered2:
448 case kPyOldPopcorn:
449 break;
8d2cd130 450 case kPyCharm:
451 case kPyCharmUnforced:
adf4d898 452 case kPyCharmPbPbMNR:
aabc7187 453 case kPyCharmpPbMNR:
e0e89f40 454 case kPyCharmppMNR:
455 case kPyCharmppMNRwmi:
8d2cd130 456 fParentSelect[0] = 411;
457 fParentSelect[1] = 421;
458 fParentSelect[2] = 431;
459 fParentSelect[3] = 4122;
460 fFlavorSelect = 4;
461 break;
adf4d898 462 case kPyD0PbPbMNR:
463 case kPyD0pPbMNR:
464 case kPyD0ppMNR:
8d2cd130 465 fParentSelect[0] = 421;
466 fFlavorSelect = 4;
467 break;
90d7b703 468 case kPyDPlusPbPbMNR:
469 case kPyDPluspPbMNR:
470 case kPyDPlusppMNR:
471 fParentSelect[0] = 411;
472 fFlavorSelect = 4;
473 break;
e0e89f40 474 case kPyDPlusStrangePbPbMNR:
475 case kPyDPlusStrangepPbMNR:
476 case kPyDPlusStrangeppMNR:
477 fParentSelect[0] = 431;
478 fFlavorSelect = 4;
479 break;
8d2cd130 480 case kPyBeauty:
adf4d898 481 case kPyBeautyPbPbMNR:
482 case kPyBeautypPbMNR:
483 case kPyBeautyppMNR:
e0e89f40 484 case kPyBeautyppMNRwmi:
8d2cd130 485 fParentSelect[0]= 511;
486 fParentSelect[1]= 521;
487 fParentSelect[2]= 531;
488 fParentSelect[3]= 5122;
489 fParentSelect[4]= 5132;
490 fParentSelect[5]= 5232;
491 fParentSelect[6]= 5332;
492 fFlavorSelect = 5;
493 break;
494 case kPyBeautyUnforced:
495 fParentSelect[0] = 511;
496 fParentSelect[1] = 521;
497 fParentSelect[2] = 531;
498 fParentSelect[3] = 5122;
499 fParentSelect[4] = 5132;
500 fParentSelect[5] = 5232;
501 fParentSelect[6] = 5332;
502 fFlavorSelect = 5;
503 break;
504 case kPyJpsiChi:
505 case kPyJpsi:
506 fParentSelect[0] = 443;
507 break;
508 case kPyMb:
509 case kPyMbNonDiffr:
d7de4a67 510 case kPyMbMSEL1:
8d2cd130 511 case kPyJets:
512 case kPyDirectGamma:
e33acb67 513 case kPyLhwgMb:
8d2cd130 514 break;
589380c6 515 case kPyW:
0f6ee828 516 case kPyZ:
589380c6 517 break;
8d2cd130 518 }
519//
592f8307 520//
521// JetFinder for Trigger
522//
523// Configure detector (EMCAL like)
524//
d7de4a67 525 fPythia->SetPARU(51, fPycellEtaMax);
526 fPythia->SetMSTU(51, fPycellNEta);
527 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 528//
529// Configure Jet Finder
530//
d7de4a67 531 fPythia->SetPARU(58, fPycellThreshold);
532 fPythia->SetPARU(52, fPycellEtSeed);
533 fPythia->SetPARU(53, fPycellMinEtJet);
534 fPythia->SetPARU(54, fPycellMaxRadius);
535 fPythia->SetMSTU(54, 2);
592f8307 536//
8d2cd130 537// This counts the total number of calls to Pyevnt() per run.
538 fTrialsRun = 0;
539 fQ = 0.;
540 fX1 = 0.;
541 fX2 = 0.;
542 fNev = 0 ;
543//
1d568bc2 544//
545//
8d2cd130 546 AliGenMC::Init();
1d568bc2 547//
548//
549//
550 if (fSetNuclei) {
551 fDyBoost = 0;
552 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
553 }
d7de4a67 554
32d6ef7d 555 if (fQuench) {
5fa4b20b 556 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 557 }
3a709cfa 558 fPythia->SetPARJ(200, 0.0);
559
b976f7d8 560 if (fQuench == 3) {
561 // Nestor's change of the splittings
562 fPythia->SetPARJ(200, 0.8);
563 fPythia->SetMSTJ(41, 1); // QCD radiation only
564 fPythia->SetMSTJ(42, 2); // angular ordering
565 fPythia->SetMSTJ(44, 2); // option to run alpha_s
566 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
567 fPythia->SetMSTJ(50, 0); // No coherence in first branching
568 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
569 }
8d2cd130 570}
571
572void AliGenPythia::Generate()
573{
574// Generate one event
e5c87a3d 575
8d2cd130 576 fDecayer->ForceDecay();
577
578 Float_t polar[3] = {0,0,0};
579 Float_t origin[3] = {0,0,0};
a920faf9 580 Float_t p[4];
8d2cd130 581// converts from mm/c to s
582 const Float_t kconv=0.001/2.999792458e8;
583//
584 Int_t nt=0;
585 Int_t jev=0;
586 Int_t j, kf;
587 fTrials=0;
f913ec4f 588 fEventTime = 0.;
589
2590ccf9 590
8d2cd130 591
592 // Set collision vertex position
2590ccf9 593 if (fVertexSmear == kPerEvent) Vertex();
594
8d2cd130 595// event loop
596 while(1)
597 {
32d6ef7d 598//
5fa4b20b 599// Produce event
32d6ef7d 600//
32d6ef7d 601//
5fa4b20b 602// Switch hadronisation off
603//
604 fPythia->SetMSTJ(1, 0);
32d6ef7d 605//
5fa4b20b 606// Either produce new event or read partons from file
607//
608 if (!fReadFromFile) {
beac474c 609 if (!fNewMIS) {
610 fPythia->Pyevnt();
611 } else {
612 fPythia->Pyevnw();
613 }
5fa4b20b 614 fNpartons = fPythia->GetN();
615 } else {
616 printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
617 fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
618 fPythia->SetN(0);
619 LoadEvent(fRL->Stack(), 0 , 1);
620 fPythia->Pyedit(21);
621 }
622
32d6ef7d 623//
624// Run quenching routine
625//
5fa4b20b 626 if (fQuench == 1) {
627 fPythia->Quench();
628 } else if (fQuench == 2){
629 fPythia->Pyquen(208., 0, 0.);
b976f7d8 630 } else if (fQuench == 3) {
631 // Quenching is via multiplicative correction of the splittings
5fa4b20b 632 }
b976f7d8 633
32d6ef7d 634//
5fa4b20b 635// Switch hadronisation on
32d6ef7d 636//
aea21c57 637 fPythia->SetMSTJ(1, 1);
5fa4b20b 638//
639// .. and perform hadronisation
aea21c57 640// printf("Calling hadronisation %d\n", fPythia->GetN());
5fa4b20b 641 fPythia->Pyexec();
8d2cd130 642 fTrials++;
8d2cd130 643 fPythia->ImportParticles(fParticles,"All");
1d568bc2 644 Boost();
8d2cd130 645//
646//
647//
648 Int_t i;
649
5fa4b20b 650
8d2cd130 651 Int_t np = fParticles->GetEntriesFast();
5fa4b20b 652
7c21f297 653 if (np == 0) continue;
8d2cd130 654//
2590ccf9 655
8d2cd130 656//
657 Int_t* pParent = new Int_t[np];
658 Int_t* pSelected = new Int_t[np];
659 Int_t* trackIt = new Int_t[np];
5fa4b20b 660 for (i = 0; i < np; i++) {
8d2cd130 661 pParent[i] = -1;
662 pSelected[i] = 0;
663 trackIt[i] = 0;
664 }
665
666 Int_t nc = 0; // Total n. of selected particles
667 Int_t nParents = 0; // Selected parents
668 Int_t nTkbles = 0; // Trackable particles
669 if (fProcess != kPyMb && fProcess != kPyJets &&
670 fProcess != kPyDirectGamma &&
589380c6 671 fProcess != kPyMbNonDiffr &&
d7de4a67 672 fProcess != kPyMbMSEL1 &&
e0e89f40 673 fProcess != kPyW && fProcess != kPyZ &&
674 fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
8d2cd130 675
5fa4b20b 676 for (i = 0; i < np; i++) {
2590ccf9 677 TParticle* iparticle = (TParticle *) fParticles->At(i);
8d2cd130 678 Int_t ks = iparticle->GetStatusCode();
679 kf = CheckPDGCode(iparticle->GetPdgCode());
680// No initial state partons
681 if (ks==21) continue;
682//
683// Heavy Flavor Selection
684//
685 // quark ?
686 kf = TMath::Abs(kf);
687 Int_t kfl = kf;
9ff6c04c 688 // Resonance
f913ec4f 689
9ff6c04c 690 if (kfl > 100000) kfl %= 100000;
183a5ca9 691 if (kfl > 10000) kfl %= 10000;
8d2cd130 692 // meson ?
693 if (kfl > 10) kfl/=100;
694 // baryon
695 if (kfl > 10) kfl/=10;
8d2cd130 696 Int_t ipa = iparticle->GetFirstMother()-1;
697 Int_t kfMo = 0;
f913ec4f 698//
699// Establish mother daughter relation between heavy quarks and mesons
700//
701 if (kf >= fFlavorSelect && kf <= 6) {
702 Int_t idau = iparticle->GetFirstDaughter() - 1;
703 if (idau > -1) {
704 TParticle* daughter = (TParticle *) fParticles->At(idau);
705 Int_t pdgD = daughter->GetPdgCode();
706 if (pdgD == 91 || pdgD == 92) {
707 Int_t jmin = daughter->GetFirstDaughter() - 1;
708 Int_t jmax = daughter->GetLastDaughter() - 1;
709 for (Int_t j = jmin; j <= jmax; j++)
710 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
711 } // is string or cluster
712 } // has daughter
713 } // heavy quark
8d2cd130 714
f913ec4f 715
8d2cd130 716 if (ipa > -1) {
717 TParticle * mother = (TParticle *) fParticles->At(ipa);
718 kfMo = TMath::Abs(mother->GetPdgCode());
719 }
f913ec4f 720
8d2cd130 721 // What to keep in Stack?
722 Bool_t flavorOK = kFALSE;
723 Bool_t selectOK = kFALSE;
724 if (fFeedDownOpt) {
32d6ef7d 725 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 726 } else {
32d6ef7d 727 if (kfl > fFlavorSelect) {
728 nc = -1;
729 break;
730 }
731 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 732 }
733 switch (fStackFillOpt) {
734 case kFlavorSelection:
32d6ef7d 735 selectOK = kTRUE;
736 break;
8d2cd130 737 case kParentSelection:
32d6ef7d 738 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
739 break;
8d2cd130 740 }
741 if (flavorOK && selectOK) {
742//
743// Heavy flavor hadron or quark
744//
745// Kinematic seletion on final state heavy flavor mesons
746 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
747 {
9ff6c04c 748 continue;
8d2cd130 749 }
750 pSelected[i] = 1;
751 if (ParentSelected(kf)) ++nParents; // Update parent count
752// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
753 } else {
754// Kinematic seletion on decay products
755 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 756 && !KinematicSelection(iparticle, 1))
8d2cd130 757 {
9ff6c04c 758 continue;
8d2cd130 759 }
760//
761// Decay products
762// Select if mother was selected and is not tracked
763
764 if (pSelected[ipa] &&
765 !trackIt[ipa] && // mother will be tracked ?
766 kfMo != 5 && // mother is b-quark, don't store fragments
767 kfMo != 4 && // mother is c-quark, don't store fragments
768 kf != 92) // don't store string
769 {
770//
771// Semi-stable or de-selected: diselect decay products:
772//
773//
774 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
775 {
776 Int_t ipF = iparticle->GetFirstDaughter();
777 Int_t ipL = iparticle->GetLastDaughter();
778 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
779 }
780// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
781 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
782 }
783 }
784 if (pSelected[i] == -1) pSelected[i] = 0;
785 if (!pSelected[i]) continue;
786 // Count quarks only if you did not include fragmentation
787 if (fFragmentation && kf <= 10) continue;
9ff6c04c 788
8d2cd130 789 nc++;
790// Decision on tracking
791 trackIt[i] = 0;
792//
793// Track final state particle
794 if (ks == 1) trackIt[i] = 1;
795// Track semi-stable particles
d25cfd65 796 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 797// Track particles selected by process if undecayed.
798 if (fForceDecay == kNoDecay) {
799 if (ParentSelected(kf)) trackIt[i] = 1;
800 } else {
801 if (ParentSelected(kf)) trackIt[i] = 0;
802 }
803 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
804//
805//
806
807 } // particle selection loop
808 if (nc > 0) {
809 for (i = 0; i<np; i++) {
810 if (!pSelected[i]) continue;
811 TParticle * iparticle = (TParticle *) fParticles->At(i);
812 kf = CheckPDGCode(iparticle->GetPdgCode());
813 Int_t ks = iparticle->GetStatusCode();
814 p[0] = iparticle->Px();
815 p[1] = iparticle->Py();
816 p[2] = iparticle->Pz();
a920faf9 817 p[3] = iparticle->Energy();
818
2590ccf9 819 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
820 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
821 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
822
8d2cd130 823 Float_t tof = kconv*iparticle->T();
824 Int_t ipa = iparticle->GetFirstMother()-1;
825 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 826
827 PushTrack(fTrackIt*trackIt[i], iparent, kf,
828 p[0], p[1], p[2], p[3],
829 origin[0], origin[1], origin[2], tof,
830 polar[0], polar[1], polar[2],
831 kPPrimary, nt, 1., ks);
8d2cd130 832 pParent[i] = nt;
833 KeepTrack(nt);
642f15cf 834 } // PushTrack loop
8d2cd130 835 }
836 } else {
837 nc = GenerateMB();
838 } // mb ?
f913ec4f 839
840 GetSubEventTime();
8d2cd130 841
234f6d32 842 delete[] pParent;
843 delete[] pSelected;
844 delete[] trackIt;
8d2cd130 845
846 if (nc > 0) {
847 switch (fCountMode) {
848 case kCountAll:
849 // printf(" Count all \n");
850 jev += nc;
851 break;
852 case kCountParents:
853 // printf(" Count parents \n");
854 jev += nParents;
855 break;
856 case kCountTrackables:
857 // printf(" Count trackable \n");
858 jev += nTkbles;
859 break;
860 }
861 if (jev >= fNpart || fNpart == -1) {
862 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 863
8d2cd130 864 fQ += fPythia->GetVINT(51);
865 fX1 += fPythia->GetVINT(41);
866 fX2 += fPythia->GetVINT(42);
867 fTrialsRun += fTrials;
868 fNev++;
869 MakeHeader();
870 break;
871 }
872 }
873 } // event loop
874 SetHighWaterMark(nt);
875// adjust weight due to kinematic selection
b88f5cea 876// AdjustWeights();
8d2cd130 877// get cross-section
878 fXsection=fPythia->GetPARI(1);
879}
880
881Int_t AliGenPythia::GenerateMB()
882{
883//
884// Min Bias selection and other global selections
885//
886 Int_t i, kf, nt, iparent;
887 Int_t nc = 0;
bf950da8 888 Float_t p[4];
8d2cd130 889 Float_t polar[3] = {0,0,0};
890 Float_t origin[3] = {0,0,0};
891// converts from mm/c to s
892 const Float_t kconv=0.001/2.999792458e8;
893
e0e89f40 894
895
5fa4b20b 896 Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
aea21c57 897
5fa4b20b 898
e0e89f40 899
8d2cd130 900 Int_t* pParent = new Int_t[np];
901 for (i=0; i< np; i++) pParent[i] = -1;
902 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
903 TParticle* jet1 = (TParticle *) fParticles->At(6);
904 TParticle* jet2 = (TParticle *) fParticles->At(7);
234f6d32 905 if (!CheckTrigger(jet1, jet2)) {
906 delete [] pParent;
907 return 0;
908 }
8d2cd130 909 }
e0e89f40 910
9fd8e520 911 // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
912 if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
ec2c406e 913
914 Bool_t ok = kFALSE;
915
916 Int_t pdg = 0;
9fd8e520 917 if (fFragPhotonInCalo) pdg = 22 ; // Photon
ec2c406e 918 else if (fPi0InCalo) pdg = 111 ; // Pi0
919
920 for (i=0; i< np; i++) {
921 TParticle* iparticle = (TParticle *) fParticles->At(i);
922 if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
9fd8e520 923 iparticle->Pt() > fFragPhotonOrPi0MinPt){
924 if(pdg == 111 ||
925 (pdg == 22 && iparticle->GetFirstMother() != 11))//No photon from hadron decay
926 {
927 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
928 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
929 if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
930 (fCheckPHOS && IsInPHOS(phi,eta)) )
931 ok =kTRUE;
932 }
ec2c406e 933 }
934 }
935 if(!ok)
936 return 0;
937 }
938
7ce1321b 939 if (fTriggerParticle) {
940 Bool_t triggered = kFALSE;
941 for (i = 0; i < np; i++) {
942 TParticle * iparticle = (TParticle *) fParticles->At(i);
943 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 944 if (kf != fTriggerParticle) continue;
7ce1321b 945 if (iparticle->Pt() == 0.) continue;
946 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
947 triggered = kTRUE;
948 break;
949 }
234f6d32 950 if (!triggered) {
951 delete [] pParent;
952 return 0;
953 }
7ce1321b 954 }
e0e89f40 955
956
957 // Check if there is a ccbar or bbbar pair with at least one of the two
958 // in fYMin < y < fYMax
959 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
960 TParticle *hvq;
961 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
962 Float_t yQ;
963 Int_t pdgQ;
964 for(i=0; i<np; i++) {
965 hvq = (TParticle*)fParticles->At(i);
966 pdgQ = hvq->GetPdgCode();
967 if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
968 if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
969 yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
970 (hvq->Energy()-hvq->Pz()+1.e-13));
971 if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
972 }
973 if (!theQ || !theQbar || !inYcut) {
234f6d32 974 delete[] pParent;
e0e89f40 975 return 0;
976 }
977 }
aea21c57 978
0f6ee828 979 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
980 if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
981 && (fCutOnChild == 1) ) {
982 if ( !CheckKinematicsOnChild() ) {
234f6d32 983 delete[] pParent;
0f6ee828 984 return 0;
985 }
aea21c57 986 }
987
988
f913ec4f 989 for (i = 0; i < np; i++) {
8d2cd130 990 Int_t trackIt = 0;
991 TParticle * iparticle = (TParticle *) fParticles->At(i);
992 kf = CheckPDGCode(iparticle->GetPdgCode());
993 Int_t ks = iparticle->GetStatusCode();
994 Int_t km = iparticle->GetFirstMother();
995 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
996 (ks != 1) ||
997 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
998 nc++;
999 if (ks == 1) trackIt = 1;
1000 Int_t ipa = iparticle->GetFirstMother()-1;
1001
1002 iparent = (ipa > -1) ? pParent[ipa] : -1;
1003
1004//
1005// store track information
1006 p[0] = iparticle->Px();
1007 p[1] = iparticle->Py();
1008 p[2] = iparticle->Pz();
a920faf9 1009 p[3] = iparticle->Energy();
1406f599 1010
a920faf9 1011
2590ccf9 1012 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1013 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1014 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1015
f913ec4f 1016 Float_t tof = fEventTime + kconv * iparticle->T();
a920faf9 1017
1018 PushTrack(fTrackIt*trackIt, iparent, kf,
1019 p[0], p[1], p[2], p[3],
1020 origin[0], origin[1], origin[2], tof,
1021 polar[0], polar[1], polar[2],
1022 kPPrimary, nt, 1., ks);
5fa4b20b 1023 //
1024 // Special Treatment to store color-flow
1025 //
1026 if (ks == 3 || ks == 13 || ks == 14) {
1027 TParticle* particle = 0;
1028 if (fStack) {
1029 particle = fStack->Particle(nt);
1030 } else {
1031 particle = gAlice->Stack()->Particle(nt);
1032 }
1033 particle->SetFirstDaughter(fPythia->GetK(2, i));
1034 particle->SetLastDaughter(fPythia->GetK(3, i));
1035 }
1036
8d2cd130 1037 KeepTrack(nt);
1038 pParent[i] = nt;
f913ec4f 1039 SetHighWaterMark(nt);
1040
8d2cd130 1041 } // select particle
1042 } // particle loop
1043
234f6d32 1044 delete[] pParent;
e0e89f40 1045
f913ec4f 1046 return 1;
8d2cd130 1047}
1048
1049
1050void AliGenPythia::FinishRun()
1051{
1052// Print x-section summary
1053 fPythia->Pystat(1);
2779fc64 1054
1055 if (fNev > 0.) {
1056 fQ /= fNev;
1057 fX1 /= fNev;
1058 fX2 /= fNev;
1059 }
1060
8d2cd130 1061 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1062 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1063}
1064
7184e472 1065void AliGenPythia::AdjustWeights() const
8d2cd130 1066{
1067// Adjust the weights after generation of all events
1068//
e2bddf81 1069 if (gAlice) {
1070 TParticle *part;
1071 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1072 for (Int_t i=0; i<ntrack; i++) {
1073 part= gAlice->GetMCApp()->Particle(i);
1074 part->SetWeight(part->GetWeight()*fKineBias);
1075 }
8d2cd130 1076 }
1077}
1078
1079void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
1080{
1081// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1082
1a626d4e 1083 fAProjectile = a1;
1084 fATarget = a2;
1d568bc2 1085 fSetNuclei = kTRUE;
8d2cd130 1086}
1087
1088
1089void AliGenPythia::MakeHeader()
1090{
7184e472 1091//
1092// Make header for the simulated event
1093//
183a5ca9 1094 if (gAlice) {
1095 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1096 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1097 }
1098
8d2cd130 1099// Builds the event header, to be called after each event
e5c87a3d 1100 if (fHeader) delete fHeader;
1101 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1102//
1103// Event type
e5c87a3d 1104 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1105//
1106// Number of trials
e5c87a3d 1107 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1108//
1109// Event Vertex
d25cfd65 1110 fHeader->SetPrimaryVertex(fVertex);
8d2cd130 1111//
1112// Jets that have triggered
f913ec4f 1113
8d2cd130 1114 if (fProcess == kPyJets)
1115 {
1116 Int_t ntrig, njet;
1117 Float_t jets[4][10];
1118 GetJets(njet, ntrig, jets);
9ff6c04c 1119
8d2cd130 1120
1121 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1122 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1123 jets[3][i]);
1124 }
1125 }
5fa4b20b 1126//
1127// Copy relevant information from external header, if present.
1128//
1129 Float_t uqJet[4];
1130
1131 if (fRL) {
1132 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1133 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1134 {
1135 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1136
1137
1138 exHeader->TriggerJet(i, uqJet);
1139 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1140 }
1141 }
1142//
1143// Store quenching parameters
1144//
1145 if (fQuench){
1146 Double_t z[4];
1147 Double_t xp, yp;
7c21f297 1148 if (fQuench == 1) {
1149 // Pythia::Quench()
1150 fPythia->GetQuenchingParameters(xp, yp, z);
1151 } else {
1152 // Pyquen
1153 Double_t r1 = PARIMP.rb1;
1154 Double_t r2 = PARIMP.rb2;
1155 Double_t b = PARIMP.b1;
1156 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1157 Double_t phi = PARIMP.psib1;
1158 xp = r * TMath::Cos(phi);
1159 yp = r * TMath::Sin(phi);
1160
1161 }
1162 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1163 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1164 }
beac474c 1165//
1166// Store pt^hard
1167 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 1168//
cf57b268 1169// Pass header
5fa4b20b 1170//
cf57b268 1171 AddHeader(fHeader);
8d2cd130 1172}
cf57b268 1173
1174void AliGenPythia::AddHeader(AliGenEventHeader* header)
1175{
1176 // Add header to container or runloader
1177 if (fContainer) {
1178 fContainer->AddHeader(header);
1179 } else {
1180 AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
1181 }
1182}
1183
8d2cd130 1184
1185Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1186{
1187// Check the kinematic trigger condition
1188//
1189 Double_t eta[2];
1190 eta[0] = jet1->Eta();
1191 eta[1] = jet2->Eta();
1192 Double_t phi[2];
1193 phi[0] = jet1->Phi();
1194 phi[1] = jet2->Phi();
1195 Int_t pdg[2];
1196 pdg[0] = jet1->GetPdgCode();
1197 pdg[1] = jet2->GetPdgCode();
1198 Bool_t triggered = kFALSE;
1199
1200 if (fProcess == kPyJets) {
1201 Int_t njets = 0;
1202 Int_t ntrig = 0;
1203 Float_t jets[4][10];
1204//
1205// Use Pythia clustering on parton level to determine jet axis
1206//
1207 GetJets(njets, ntrig, jets);
1208
76d6ba9a 1209 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1210//
1211 } else {
1212 Int_t ij = 0;
1213 Int_t ig = 1;
1214 if (pdg[0] == kGamma) {
1215 ij = 1;
1216 ig = 0;
1217 }
1218 //Check eta range first...
1219 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1220 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1221 {
1222 //Eta is okay, now check phi range
1223 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1224 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1225 {
1226 triggered = kTRUE;
1227 }
1228 }
1229 }
1230 return triggered;
1231}
aea21c57 1232
1233
aea21c57 1234
7184e472 1235Bool_t AliGenPythia::CheckKinematicsOnChild(){
1236//
1237//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1238//
aea21c57 1239 Bool_t checking = kFALSE;
1240 Int_t j, kcode, ks, km;
1241 Int_t nPartAcc = 0; //number of particles in the acceptance range
1242 Int_t numberOfAcceptedParticles = 1;
1243 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1244 Int_t npart = fParticles->GetEntriesFast();
1245
0f6ee828 1246 for (j = 0; j<npart; j++) {
aea21c57 1247 TParticle * jparticle = (TParticle *) fParticles->At(j);
1248 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1249 ks = jparticle->GetStatusCode();
1250 km = jparticle->GetFirstMother();
1251
1252 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1253 nPartAcc++;
1254 }
0f6ee828 1255 if( numberOfAcceptedParticles <= nPartAcc){
1256 checking = kTRUE;
1257 break;
1258 }
aea21c57 1259 }
0f6ee828 1260
aea21c57 1261 return checking;
aea21c57 1262}
1263
8d2cd130 1264
1265AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
1266{
1267// Assignment operator
014a9521 1268 rhs.Copy(*this);
8d2cd130 1269 return *this;
1270}
1271
5fa4b20b 1272void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1273{
1274//
1275// Load event into Pythia Common Block
1276//
5fa4b20b 1277
32d6ef7d 1278 Int_t npart = stack -> GetNprimary();
1279 Int_t n0 = 0;
1280
1281 if (!flag) {
1282 (fPythia->GetPyjets())->N = npart;
1283 } else {
1284 n0 = (fPythia->GetPyjets())->N;
1285 (fPythia->GetPyjets())->N = n0 + npart;
1286 }
1287
1288
8d2cd130 1289 for (Int_t part = 0; part < npart; part++) {
7184e472 1290 TParticle *mPart = stack->Particle(part);
32d6ef7d 1291
7184e472 1292 Int_t kf = mPart->GetPdgCode();
1293 Int_t ks = mPart->GetStatusCode();
1294 Int_t idf = mPart->GetFirstDaughter();
1295 Int_t idl = mPart->GetLastDaughter();
5fa4b20b 1296
1297 if (reHadr) {
1298 if (ks == 11 || ks == 12) {
1299 ks -= 10;
1300 idf = -1;
1301 idl = -1;
1302 }
1303 }
32d6ef7d 1304
7184e472 1305 Float_t px = mPart->Px();
1306 Float_t py = mPart->Py();
1307 Float_t pz = mPart->Pz();
1308 Float_t e = mPart->Energy();
1309 Float_t m = mPart->GetCalcMass();
8d2cd130 1310
1311
32d6ef7d 1312 (fPythia->GetPyjets())->P[0][part+n0] = px;
1313 (fPythia->GetPyjets())->P[1][part+n0] = py;
1314 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1315 (fPythia->GetPyjets())->P[3][part+n0] = e;
1316 (fPythia->GetPyjets())->P[4][part+n0] = m;
8d2cd130 1317
32d6ef7d 1318 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1319 (fPythia->GetPyjets())->K[0][part+n0] = ks;
5fa4b20b 1320 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1321 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
7184e472 1322 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
8d2cd130 1323 }
1324}
1325
5fa4b20b 1326
014a9521 1327void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1328{
1329//
1330// Calls the Pythia jet finding algorithm to find jets in the current event
1331//
1332//
8d2cd130 1333//
1334// Save jets
1335 Int_t n = fPythia->GetN();
1336
1337//
1338// Run Jet Finder
1339 fPythia->Pycell(njets);
1340 Int_t i;
1341 for (i = 0; i < njets; i++) {
1342 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1343 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1344 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1345 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1346
1347 jets[0][i] = px;
1348 jets[1][i] = py;
1349 jets[2][i] = pz;
1350 jets[3][i] = e;
1351 }
1352}
1353
1354
1355
1356void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1357{
1358//
1359// Calls the Pythia clustering algorithm to find jets in the current event
1360//
1361 Int_t n = fPythia->GetN();
1362 nJets = 0;
1363 nJetsTrig = 0;
1364 if (fJetReconstruction == kCluster) {
1365//
1366// Configure cluster algorithm
1367//
1368 fPythia->SetPARU(43, 2.);
1369 fPythia->SetMSTU(41, 1);
1370//
1371// Call cluster algorithm
1372//
1373 fPythia->Pyclus(nJets);
1374//
1375// Loading jets from common block
1376//
1377 } else {
592f8307 1378
8d2cd130 1379//
1380// Run Jet Finder
1381 fPythia->Pycell(nJets);
1382 }
1383
1384 Int_t i;
1385 for (i = 0; i < nJets; i++) {
1386 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1387 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1388 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1389 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1390 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1391 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1392 Float_t theta = TMath::ATan2(pt,pz);
1393 Float_t et = e * TMath::Sin(theta);
1394 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1395 if (
1396 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1397 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1398 et > fEtMinJet && et < fEtMaxJet
1399 )
1400 {
1401 jets[0][nJetsTrig] = px;
1402 jets[1][nJetsTrig] = py;
1403 jets[2][nJetsTrig] = pz;
1404 jets[3][nJetsTrig] = e;
1405 nJetsTrig++;
5fa4b20b 1406// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1407 } else {
1408// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1409 }
1410 }
1411}
1412
f913ec4f 1413void AliGenPythia::GetSubEventTime()
1414{
1415 // Calculates time of the next subevent
9d7108a7 1416 fEventTime = 0.;
1417 if (fEventsTime) {
1418 TArrayF &array = *fEventsTime;
1419 fEventTime = array[fCurSubEvent++];
1420 }
1421 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1422 return;
f913ec4f 1423}
8d2cd130 1424
ec2c406e 1425Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1426{
1427 // Is particle in EMCAL acceptance?
ec2c406e 1428 // phi in degrees, etamin=-etamax
9fd8e520 1429 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1430 eta < fEMCALEta )
ec2c406e 1431 return kTRUE;
1432 else
1433 return kFALSE;
1434}
1435
1436Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1437{
1438 // Is particle in PHOS acceptance?
1439 // Acceptance slightly larger considered.
1440 // phi in degrees, etamin=-etamax
9fd8e520 1441 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1442 eta < fPHOSEta )
ec2c406e 1443 return kTRUE;
1444 else
1445 return kFALSE;
1446}
1447
1448
1449
8d2cd130 1450#ifdef never
1451void AliGenPythia::Streamer(TBuffer &R__b)
1452{
1453 // Stream an object of class AliGenPythia.
1454
1455 if (R__b.IsReading()) {
1456 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1457 AliGenerator::Streamer(R__b);
1458 R__b >> (Int_t&)fProcess;
1459 R__b >> (Int_t&)fStrucFunc;
1460 R__b >> (Int_t&)fForceDecay;
1461 R__b >> fEnergyCMS;
1462 R__b >> fKineBias;
1463 R__b >> fTrials;
1464 fParentSelect.Streamer(R__b);
1465 fChildSelect.Streamer(R__b);
1466 R__b >> fXsection;
1467// (AliPythia::Instance())->Streamer(R__b);
1468 R__b >> fPtHardMin;
1469 R__b >> fPtHardMax;
1470// if (fDecayer) fDecayer->Streamer(R__b);
1471 } else {
1472 R__b.WriteVersion(AliGenPythia::IsA());
1473 AliGenerator::Streamer(R__b);
1474 R__b << (Int_t)fProcess;
1475 R__b << (Int_t)fStrucFunc;
1476 R__b << (Int_t)fForceDecay;
1477 R__b << fEnergyCMS;
1478 R__b << fKineBias;
1479 R__b << fTrials;
1480 fParentSelect.Streamer(R__b);
1481 fChildSelect.Streamer(R__b);
1482 R__b << fXsection;
1483// R__b << fPythia;
1484 R__b << fPtHardMin;
1485 R__b << fPtHardMax;
1486 // fDecayer->Streamer(R__b);
1487 }
1488}
1489#endif
1490
90d7b703 1491
589380c6 1492