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