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