git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bringing the parametrized TPC tracker up-to-date (A.Dainese)
[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),
114 fBremssInCalo(kFALSE),
115 fPi0InCalo(kFALSE) ,
116 fBremssPi0InEMCAL(kFALSE),
117 fBremssPi0InPHOS(kFALSE),
118 fBremssPi0MinPt(0),
119 fPHOSMinPhi(220.),
120 fPHOSMaxPhi(320.),
121 fPHOSEta(0.12),
122 fEMCALMinPhi(80.),
123 fEMCALMaxPhi(190.),
124 fEMCALEta(0.7)
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),
198 fBremssInCalo(kFALSE),
199 fPi0InCalo(kFALSE) ,
200 fBremssPi0InEMCAL(kFALSE),
201 fBremssPi0InPHOS(kFALSE),
202 fBremssPi0MinPt(0),
203 fPHOSMinPhi(220.),
204 fPHOSMaxPhi(320.),
205 fPHOSEta(0.12),
206 fEMCALMinPhi(80.),
207 fEMCALMaxPhi(190.),
208 fEMCALEta(0.7)
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),
289 fBremssInCalo(kFALSE),
290 fPi0InCalo(kFALSE) ,
291 fBremssPi0InEMCAL(kFALSE),
292 fBremssPi0InPHOS(kFALSE),
293 fBremssPi0MinPt(0),
294 fPHOSMinPhi(220.),
295 fPHOSMaxPhi(320.),
296 fPHOSEta(0.12),
297 fEMCALMinPhi(80.),
298 fEMCALMaxPhi(190.),
299 fEMCALEta(0.7)
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
ec2c406e 911 // Select jets with bremsstrahlung or pi0 going to PHOS or EMCAL
912 if (fProcess == kPyJets && (fBremssInCalo || fPi0InCalo) ) {
913
914 Bool_t ok = kFALSE;
915
916 Int_t pdg = 0;
917 if (fBremssInCalo) pdg = 22 ; // Photon
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 &&
923 iparticle->Pt() > fBremssPi0MinPt){
924 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
925 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
926 if((fBremssPi0InEMCAL && IsInEMCAL(phi,eta)) ||
927 fBremssPi0InPHOS && IsInPHOS(phi,eta))
928 ok =kTRUE;
929 }
930 }
931 if(!ok)
932 return 0;
933 }
934
7ce1321b 935 if (fTriggerParticle) {
936 Bool_t triggered = kFALSE;
937 for (i = 0; i < np; i++) {
938 TParticle * iparticle = (TParticle *) fParticles->At(i);
939 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 940 if (kf != fTriggerParticle) continue;
7ce1321b 941 if (iparticle->Pt() == 0.) continue;
942 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
943 triggered = kTRUE;
944 break;
945 }
234f6d32 946 if (!triggered) {
947 delete [] pParent;
948 return 0;
949 }
7ce1321b 950 }
e0e89f40 951
952
953 // Check if there is a ccbar or bbbar pair with at least one of the two
954 // in fYMin < y < fYMax
955 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
956 TParticle *hvq;
957 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
958 Float_t yQ;
959 Int_t pdgQ;
960 for(i=0; i<np; i++) {
961 hvq = (TParticle*)fParticles->At(i);
962 pdgQ = hvq->GetPdgCode();
963 if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
964 if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
965 yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
966 (hvq->Energy()-hvq->Pz()+1.e-13));
967 if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
968 }
969 if (!theQ || !theQbar || !inYcut) {
234f6d32 970 delete[] pParent;
e0e89f40 971 return 0;
972 }
973 }
aea21c57 974
0f6ee828 975 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
976 if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
977 && (fCutOnChild == 1) ) {
978 if ( !CheckKinematicsOnChild() ) {
234f6d32 979 delete[] pParent;
0f6ee828 980 return 0;
981 }
aea21c57 982 }
983
984
f913ec4f 985 for (i = 0; i < np; i++) {
8d2cd130 986 Int_t trackIt = 0;
987 TParticle * iparticle = (TParticle *) fParticles->At(i);
988 kf = CheckPDGCode(iparticle->GetPdgCode());
989 Int_t ks = iparticle->GetStatusCode();
990 Int_t km = iparticle->GetFirstMother();
991 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
992 (ks != 1) ||
993 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
994 nc++;
995 if (ks == 1) trackIt = 1;
996 Int_t ipa = iparticle->GetFirstMother()-1;
997
998 iparent = (ipa > -1) ? pParent[ipa] : -1;
999
1000//
1001// store track information
1002 p[0] = iparticle->Px();
1003 p[1] = iparticle->Py();
1004 p[2] = iparticle->Pz();
a920faf9 1005 p[3] = iparticle->Energy();
1406f599 1006
a920faf9 1007
2590ccf9 1008 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1009 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1010 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1011
f913ec4f 1012 Float_t tof = fEventTime + kconv * iparticle->T();
a920faf9 1013
1014 PushTrack(fTrackIt*trackIt, iparent, kf,
1015 p[0], p[1], p[2], p[3],
1016 origin[0], origin[1], origin[2], tof,
1017 polar[0], polar[1], polar[2],
1018 kPPrimary, nt, 1., ks);
5fa4b20b 1019 //
1020 // Special Treatment to store color-flow
1021 //
1022 if (ks == 3 || ks == 13 || ks == 14) {
1023 TParticle* particle = 0;
1024 if (fStack) {
1025 particle = fStack->Particle(nt);
1026 } else {
1027 particle = gAlice->Stack()->Particle(nt);
1028 }
1029 particle->SetFirstDaughter(fPythia->GetK(2, i));
1030 particle->SetLastDaughter(fPythia->GetK(3, i));
1031 }
1032
8d2cd130 1033 KeepTrack(nt);
1034 pParent[i] = nt;
f913ec4f 1035 SetHighWaterMark(nt);
1036
8d2cd130 1037 } // select particle
1038 } // particle loop
1039
234f6d32 1040 delete[] pParent;
e0e89f40 1041
f913ec4f 1042 return 1;
8d2cd130 1043}
1044
1045
1046void AliGenPythia::FinishRun()
1047{
1048// Print x-section summary
1049 fPythia->Pystat(1);
2779fc64 1050
1051 if (fNev > 0.) {
1052 fQ /= fNev;
1053 fX1 /= fNev;
1054 fX2 /= fNev;
1055 }
1056
8d2cd130 1057 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1058 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1059}
1060
7184e472 1061void AliGenPythia::AdjustWeights() const
8d2cd130 1062{
1063// Adjust the weights after generation of all events
1064//
e2bddf81 1065 if (gAlice) {
1066 TParticle *part;
1067 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1068 for (Int_t i=0; i<ntrack; i++) {
1069 part= gAlice->GetMCApp()->Particle(i);
1070 part->SetWeight(part->GetWeight()*fKineBias);
1071 }
8d2cd130 1072 }
1073}
1074
1075void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
1076{
1077// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1078
1a626d4e 1079 fAProjectile = a1;
1080 fATarget = a2;
1d568bc2 1081 fSetNuclei = kTRUE;
8d2cd130 1082}
1083
1084
1085void AliGenPythia::MakeHeader()
1086{
7184e472 1087//
1088// Make header for the simulated event
1089//
183a5ca9 1090 if (gAlice) {
1091 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1092 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1093 }
1094
8d2cd130 1095// Builds the event header, to be called after each event
e5c87a3d 1096 if (fHeader) delete fHeader;
1097 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1098//
1099// Event type
e5c87a3d 1100 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1101//
1102// Number of trials
e5c87a3d 1103 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1104//
1105// Event Vertex
d25cfd65 1106 fHeader->SetPrimaryVertex(fVertex);
8d2cd130 1107//
1108// Jets that have triggered
f913ec4f 1109
8d2cd130 1110 if (fProcess == kPyJets)
1111 {
1112 Int_t ntrig, njet;
1113 Float_t jets[4][10];
1114 GetJets(njet, ntrig, jets);
9ff6c04c 1115
8d2cd130 1116
1117 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1118 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1119 jets[3][i]);
1120 }
1121 }
5fa4b20b 1122//
1123// Copy relevant information from external header, if present.
1124//
1125 Float_t uqJet[4];
1126
1127 if (fRL) {
1128 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1129 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1130 {
1131 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1132
1133
1134 exHeader->TriggerJet(i, uqJet);
1135 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1136 }
1137 }
1138//
1139// Store quenching parameters
1140//
1141 if (fQuench){
1142 Double_t z[4];
1143 Double_t xp, yp;
7c21f297 1144 if (fQuench == 1) {
1145 // Pythia::Quench()
1146 fPythia->GetQuenchingParameters(xp, yp, z);
1147 } else {
1148 // Pyquen
1149 Double_t r1 = PARIMP.rb1;
1150 Double_t r2 = PARIMP.rb2;
1151 Double_t b = PARIMP.b1;
1152 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1153 Double_t phi = PARIMP.psib1;
1154 xp = r * TMath::Cos(phi);
1155 yp = r * TMath::Sin(phi);
1156
1157 }
1158 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1159 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1160 }
beac474c 1161//
1162// Store pt^hard
1163 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 1164//
cf57b268 1165// Pass header
5fa4b20b 1166//
cf57b268 1167 AddHeader(fHeader);
8d2cd130 1168}
cf57b268 1169
1170void AliGenPythia::AddHeader(AliGenEventHeader* header)
1171{
1172 // Add header to container or runloader
1173 if (fContainer) {
1174 fContainer->AddHeader(header);
1175 } else {
1176 AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
1177 }
1178}
1179
8d2cd130 1180
1181Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1182{
1183// Check the kinematic trigger condition
1184//
1185 Double_t eta[2];
1186 eta[0] = jet1->Eta();
1187 eta[1] = jet2->Eta();
1188 Double_t phi[2];
1189 phi[0] = jet1->Phi();
1190 phi[1] = jet2->Phi();
1191 Int_t pdg[2];
1192 pdg[0] = jet1->GetPdgCode();
1193 pdg[1] = jet2->GetPdgCode();
1194 Bool_t triggered = kFALSE;
1195
1196 if (fProcess == kPyJets) {
1197 Int_t njets = 0;
1198 Int_t ntrig = 0;
1199 Float_t jets[4][10];
1200//
1201// Use Pythia clustering on parton level to determine jet axis
1202//
1203 GetJets(njets, ntrig, jets);
1204
76d6ba9a 1205 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1206//
1207 } else {
1208 Int_t ij = 0;
1209 Int_t ig = 1;
1210 if (pdg[0] == kGamma) {
1211 ij = 1;
1212 ig = 0;
1213 }
1214 //Check eta range first...
1215 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1216 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1217 {
1218 //Eta is okay, now check phi range
1219 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1220 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1221 {
1222 triggered = kTRUE;
1223 }
1224 }
1225 }
1226 return triggered;
1227}
aea21c57 1228
1229
aea21c57 1230
7184e472 1231Bool_t AliGenPythia::CheckKinematicsOnChild(){
1232//
1233//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1234//
aea21c57 1235 Bool_t checking = kFALSE;
1236 Int_t j, kcode, ks, km;
1237 Int_t nPartAcc = 0; //number of particles in the acceptance range
1238 Int_t numberOfAcceptedParticles = 1;
1239 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1240 Int_t npart = fParticles->GetEntriesFast();
1241
0f6ee828 1242 for (j = 0; j<npart; j++) {
aea21c57 1243 TParticle * jparticle = (TParticle *) fParticles->At(j);
1244 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1245 ks = jparticle->GetStatusCode();
1246 km = jparticle->GetFirstMother();
1247
1248 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1249 nPartAcc++;
1250 }
0f6ee828 1251 if( numberOfAcceptedParticles <= nPartAcc){
1252 checking = kTRUE;
1253 break;
1254 }
aea21c57 1255 }
0f6ee828 1256
aea21c57 1257 return checking;
aea21c57 1258}
1259
8d2cd130 1260
1261AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
1262{
1263// Assignment operator
014a9521 1264 rhs.Copy(*this);
8d2cd130 1265 return *this;
1266}
1267
5fa4b20b 1268void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1269{
1270//
1271// Load event into Pythia Common Block
1272//
5fa4b20b 1273
32d6ef7d 1274 Int_t npart = stack -> GetNprimary();
1275 Int_t n0 = 0;
1276
1277 if (!flag) {
1278 (fPythia->GetPyjets())->N = npart;
1279 } else {
1280 n0 = (fPythia->GetPyjets())->N;
1281 (fPythia->GetPyjets())->N = n0 + npart;
1282 }
1283
1284
8d2cd130 1285 for (Int_t part = 0; part < npart; part++) {
7184e472 1286 TParticle *mPart = stack->Particle(part);
32d6ef7d 1287
7184e472 1288 Int_t kf = mPart->GetPdgCode();
1289 Int_t ks = mPart->GetStatusCode();
1290 Int_t idf = mPart->GetFirstDaughter();
1291 Int_t idl = mPart->GetLastDaughter();
5fa4b20b 1292
1293 if (reHadr) {
1294 if (ks == 11 || ks == 12) {
1295 ks -= 10;
1296 idf = -1;
1297 idl = -1;
1298 }
1299 }
32d6ef7d 1300
7184e472 1301 Float_t px = mPart->Px();
1302 Float_t py = mPart->Py();
1303 Float_t pz = mPart->Pz();
1304 Float_t e = mPart->Energy();
1305 Float_t m = mPart->GetCalcMass();
8d2cd130 1306
1307
32d6ef7d 1308 (fPythia->GetPyjets())->P[0][part+n0] = px;
1309 (fPythia->GetPyjets())->P[1][part+n0] = py;
1310 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1311 (fPythia->GetPyjets())->P[3][part+n0] = e;
1312 (fPythia->GetPyjets())->P[4][part+n0] = m;
8d2cd130 1313
32d6ef7d 1314 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1315 (fPythia->GetPyjets())->K[0][part+n0] = ks;
5fa4b20b 1316 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1317 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
7184e472 1318 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
8d2cd130 1319 }
1320}
1321
5fa4b20b 1322
014a9521 1323void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1324{
1325//
1326// Calls the Pythia jet finding algorithm to find jets in the current event
1327//
1328//
8d2cd130 1329//
1330// Save jets
1331 Int_t n = fPythia->GetN();
1332
1333//
1334// Run Jet Finder
1335 fPythia->Pycell(njets);
1336 Int_t i;
1337 for (i = 0; i < njets; i++) {
1338 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1339 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1340 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1341 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1342
1343 jets[0][i] = px;
1344 jets[1][i] = py;
1345 jets[2][i] = pz;
1346 jets[3][i] = e;
1347 }
1348}
1349
1350
1351
1352void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1353{
1354//
1355// Calls the Pythia clustering algorithm to find jets in the current event
1356//
1357 Int_t n = fPythia->GetN();
1358 nJets = 0;
1359 nJetsTrig = 0;
1360 if (fJetReconstruction == kCluster) {
1361//
1362// Configure cluster algorithm
1363//
1364 fPythia->SetPARU(43, 2.);
1365 fPythia->SetMSTU(41, 1);
1366//
1367// Call cluster algorithm
1368//
1369 fPythia->Pyclus(nJets);
1370//
1371// Loading jets from common block
1372//
1373 } else {
592f8307 1374
8d2cd130 1375//
1376// Run Jet Finder
1377 fPythia->Pycell(nJets);
1378 }
1379
1380 Int_t i;
1381 for (i = 0; i < nJets; i++) {
1382 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1383 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1384 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1385 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1386 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1387 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1388 Float_t theta = TMath::ATan2(pt,pz);
1389 Float_t et = e * TMath::Sin(theta);
1390 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1391 if (
1392 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1393 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1394 et > fEtMinJet && et < fEtMaxJet
1395 )
1396 {
1397 jets[0][nJetsTrig] = px;
1398 jets[1][nJetsTrig] = py;
1399 jets[2][nJetsTrig] = pz;
1400 jets[3][nJetsTrig] = e;
1401 nJetsTrig++;
5fa4b20b 1402// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1403 } else {
1404// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1405 }
1406 }
1407}
1408
f913ec4f 1409void AliGenPythia::GetSubEventTime()
1410{
1411 // Calculates time of the next subevent
9d7108a7 1412 fEventTime = 0.;
1413 if (fEventsTime) {
1414 TArrayF &array = *fEventsTime;
1415 fEventTime = array[fCurSubEvent++];
1416 }
1417 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1418 return;
f913ec4f 1419}
8d2cd130 1420
ec2c406e 1421Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1422{
1423 // Is particle in EMCAL acceptance?
1424 // Acceptance slightly larger considered.
1425 // phi in degrees, etamin=-etamax
1426 if(phi > fEMCALMinPhi-0.1 && phi < fEMCALMaxPhi+0.1 &&
1427 eta < fEMCALEta+0.01 )
1428 return kTRUE;
1429 else
1430 return kFALSE;
1431}
1432
1433Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1434{
1435 // Is particle in PHOS acceptance?
1436 // Acceptance slightly larger considered.
1437 // phi in degrees, etamin=-etamax
1438 if(phi > fPHOSMinPhi-0.1 && phi < fPHOSMaxPhi+0.1 &&
1439 eta < fPHOSEta+0.01 )
1440 return kTRUE;
1441 else
1442 return kFALSE;
1443}
1444
1445
1446
8d2cd130 1447#ifdef never
1448void AliGenPythia::Streamer(TBuffer &R__b)
1449{
1450 // Stream an object of class AliGenPythia.
1451
1452 if (R__b.IsReading()) {
1453 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1454 AliGenerator::Streamer(R__b);
1455 R__b >> (Int_t&)fProcess;
1456 R__b >> (Int_t&)fStrucFunc;
1457 R__b >> (Int_t&)fForceDecay;
1458 R__b >> fEnergyCMS;
1459 R__b >> fKineBias;
1460 R__b >> fTrials;
1461 fParentSelect.Streamer(R__b);
1462 fChildSelect.Streamer(R__b);
1463 R__b >> fXsection;
1464// (AliPythia::Instance())->Streamer(R__b);
1465 R__b >> fPtHardMin;
1466 R__b >> fPtHardMax;
1467// if (fDecayer) fDecayer->Streamer(R__b);
1468 } else {
1469 R__b.WriteVersion(AliGenPythia::IsA());
1470 AliGenerator::Streamer(R__b);
1471 R__b << (Int_t)fProcess;
1472 R__b << (Int_t)fStrucFunc;
1473 R__b << (Int_t)fForceDecay;
1474 R__b << fEnergyCMS;
1475 R__b << fKineBias;
1476 R__b << fTrials;
1477 fParentSelect.Streamer(R__b);
1478 fChildSelect.Streamer(R__b);
1479 R__b << fXsection;
1480// R__b << fPythia;
1481 R__b << fPtHardMin;
1482 R__b << fPtHardMax;
1483 // fDecayer->Streamer(R__b);
1484 }
1485}
1486#endif
1487
90d7b703 1488
589380c6 1489