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