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