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