EffC++ warnings corrected.
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
CommitLineData
8d2cd130 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
7cdba479 16/* $Id$ */
8d2cd130 17
18//
19// Generator using the TPythia interface (via AliPythia)
20// to generate pp collisions.
21// Using SetNuclei() also nuclear modifications to the structure functions
22// can be taken into account. This makes, of course, only sense for the
23// generation of the products of hard processes (heavy flavor, jets ...)
24//
25// andreas.morsch@cern.ch
26//
27
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),
113 fFileName(0)
8d2cd130 114{
115// Default Constructor
2d9f9817 116 SetNuclei(0,0);
7cdba479 117 if (!AliPythiaRndm::GetPythiaRandom())
e8a8adcd 118 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 119}
120
121AliGenPythia::AliGenPythia(Int_t npart)
e8a8adcd 122 :AliGenMC(npart),
123 fProcess(kPyCharm),
124 fStrucFunc(kCTEQ5L),
125 fEnergyCMS(5500.),
126 fKineBias(0.),
127 fTrials(0),
128 fTrialsRun(0),
129 fQ(0.),
130 fX1(0.),
131 fX2(0.),
132 fEventTime(0.),
133 fInteractionRate(0.),
134 fTimeWindow(0.),
135 fCurSubEvent(0),
136 fEventsTime(0),
137 fNev(0),
138 fFlavorSelect(0),
139 fXsection(0.),
140 fPythia(0),
141 fPtHardMin(0.),
142 fPtHardMax(1.e4),
143 fYHardMin(-1.e10),
144 fYHardMax(1.e10),
145 fGinit(kTRUE),
146 fGfinal(kTRUE),
147 fHadronisation(kTRUE),
148 fNpartons(0),
149 fReadFromFile(kFALSE),
150 fQuench(kFALSE),
151 fPtKick(1.),
152 fFullEvent(kTRUE),
153 fDecayer(new AliDecayerPythia()),
154 fDebugEventFirst(-1),
155 fDebugEventLast(-1),
156 fEtMinJet(0.),
157 fEtMaxJet(1.e4),
158 fEtaMinJet(-20.),
159 fEtaMaxJet(20.),
160 fPhiMinJet(0.),
161 fPhiMaxJet(2.* TMath::Pi()),
162 fJetReconstruction(kCell),
163 fEtaMinGamma(-20.),
164 fEtaMaxGamma(20.),
165 fPhiMinGamma(0.),
166 fPhiMaxGamma(2. * TMath::Pi()),
167 fPycellEtaMax(2.),
168 fPycellNEta(274),
169 fPycellNPhi(432),
170 fPycellThreshold(0.),
171 fPycellEtSeed(4.),
172 fPycellMinEtJet(10.),
173 fPycellMaxRadius(1.),
174 fStackFillOpt(kFlavorSelection),
175 fFeedDownOpt(kTRUE),
176 fFragmentation(kTRUE),
177 fSetNuclei(kFALSE),
178 fNewMIS(kFALSE),
179 fHFoff(kFALSE),
180 fTriggerParticle(0),
181 fTriggerEta(0.9),
182 fCountMode(kCountAll),
183 fHeader(0),
184 fRL(0),
185 fFileName(0)
8d2cd130 186{
187// default charm production at 5. 5 TeV
188// semimuonic decay
189// structure function GRVHO
190//
191 fName = "Pythia";
192 fTitle= "Particle Generator using PYTHIA";
8d2cd130 193 SetForceDecay();
8d2cd130 194 // Set random number generator
7cdba479 195 if (!AliPythiaRndm::GetPythiaRandom())
196 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 197 fParticles = new TClonesArray("TParticle",1000);
2d9f9817 198 SetNuclei(0,0);
76d6ba9a 199 }
8d2cd130 200
201AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
e8a8adcd 202 :AliGenMC(Pythia),
203 fProcess(kPyCharm),
204 fStrucFunc(kCTEQ5L),
205 fEnergyCMS(5500.),
206 fKineBias(0.),
207 fTrials(0),
208 fTrialsRun(0),
209 fQ(0.),
210 fX1(0.),
211 fX2(0.),
212 fEventTime(0.),
213 fInteractionRate(0.),
214 fTimeWindow(0.),
215 fCurSubEvent(0),
216 fEventsTime(0),
217 fNev(0),
218 fFlavorSelect(0),
219 fXsection(0.),
220 fPythia(0),
221 fPtHardMin(0.),
222 fPtHardMax(1.e4),
223 fYHardMin(-1.e10),
224 fYHardMax(1.e10),
225 fGinit(kTRUE),
226 fGfinal(kTRUE),
227 fHadronisation(kTRUE),
228 fNpartons(0),
229 fReadFromFile(kFALSE),
230 fQuench(kFALSE),
231 fPtKick(1.),
232 fFullEvent(kTRUE),
233 fDecayer(new AliDecayerPythia()),
234 fDebugEventFirst(-1),
235 fDebugEventLast(-1),
236 fEtMinJet(0.),
237 fEtMaxJet(1.e4),
238 fEtaMinJet(-20.),
239 fEtaMaxJet(20.),
240 fPhiMinJet(0.),
241 fPhiMaxJet(2.* TMath::Pi()),
242 fJetReconstruction(kCell),
243 fEtaMinGamma(-20.),
244 fEtaMaxGamma(20.),
245 fPhiMinGamma(0.),
246 fPhiMaxGamma(2. * TMath::Pi()),
247 fPycellEtaMax(2.),
248 fPycellNEta(274),
249 fPycellNPhi(432),
250 fPycellThreshold(0.),
251 fPycellEtSeed(4.),
252 fPycellMinEtJet(10.),
253 fPycellMaxRadius(1.),
254 fStackFillOpt(kFlavorSelection),
255 fFeedDownOpt(kTRUE),
256 fFragmentation(kTRUE),
257 fSetNuclei(kFALSE),
258 fNewMIS(kFALSE),
259 fHFoff(kFALSE),
260 fTriggerParticle(0),
261 fTriggerEta(0.9),
262 fCountMode(kCountAll),
263 fHeader(0),
264 fRL(0),
265 fFileName(0)
8d2cd130 266{
267// copy constructor
268 Pythia.Copy(*this);
269}
270
271AliGenPythia::~AliGenPythia()
272{
273// Destructor
9d7108a7 274 if(fEventsTime) delete fEventsTime;
275}
276
277void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
278{
279// Generate pileup using user specified rate
280 fInteractionRate = rate;
281 fTimeWindow = timewindow;
282 GeneratePileup();
283}
284
285void AliGenPythia::GeneratePileup()
286{
287// Generate sub events time for pileup
288 fEventsTime = 0;
289 if(fInteractionRate == 0.) {
290 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
291 return;
292 }
293
294 Int_t npart = NumberParticles();
295 if(npart < 0) {
296 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
297 return;
298 }
299
300 if(fEventsTime) delete fEventsTime;
301 fEventsTime = new TArrayF(npart);
302 TArrayF &array = *fEventsTime;
303 for(Int_t ipart = 0; ipart < npart; ipart++)
304 array[ipart] = 0.;
305
306 Float_t eventtime = 0.;
307 while(1)
308 {
309 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
310 if(eventtime > fTimeWindow) break;
311 array.Set(array.GetSize()+1);
312 array[array.GetSize()-1] = eventtime;
313 }
314
315 eventtime = 0.;
316 while(1)
317 {
318 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
319 if(TMath::Abs(eventtime) > fTimeWindow) break;
320 array.Set(array.GetSize()+1);
321 array[array.GetSize()-1] = eventtime;
322 }
323
324 SetNumberParticles(fEventsTime->GetSize());
8d2cd130 325}
326
592f8307 327void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
328 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
329{
330// Set pycell parameters
331 fPycellEtaMax = etamax;
332 fPycellNEta = neta;
333 fPycellNPhi = nphi;
334 fPycellThreshold = thresh;
335 fPycellEtSeed = etseed;
336 fPycellMinEtJet = minet;
337 fPycellMaxRadius = r;
338}
339
340
341
8d2cd130 342void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
343{
344 // Set a range of event numbers, for which a table
345 // of generated particle will be printed
346 fDebugEventFirst = eventFirst;
347 fDebugEventLast = eventLast;
348 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
349}
350
351void AliGenPythia::Init()
352{
353// Initialisation
354
355 SetMC(AliPythia::Instance());
b88f5cea 356 fPythia=(AliPythia*) fMCEvGen;
e2bddf81 357
8d2cd130 358//
359 fParentWeight=1./Float_t(fNpart);
360//
361// Forward Paramters to the AliPythia object
362 fDecayer->SetForceDecay(fForceDecay);
363 fDecayer->Init();
364
365
366 fPythia->SetCKIN(3,fPtHardMin);
367 fPythia->SetCKIN(4,fPtHardMax);
368 fPythia->SetCKIN(7,fYHardMin);
369 fPythia->SetCKIN(8,fYHardMax);
370
1a626d4e 371 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
8d2cd130 372 // Fragmentation?
373 if (fFragmentation) {
374 fPythia->SetMSTP(111,1);
375 } else {
376 fPythia->SetMSTP(111,0);
377 }
378
379
380// initial state radiation
381 fPythia->SetMSTP(61,fGinit);
382// final state radiation
383 fPythia->SetMSTP(71,fGfinal);
384// pt - kick
385 if (fPtKick > 0.) {
386 fPythia->SetMSTP(91,1);
387 fPythia->SetPARP(91,fPtKick);
388 } else {
389 fPythia->SetMSTP(91,0);
390 }
391
5fa4b20b 392
393 if (fReadFromFile) {
394 fRL = AliRunLoader::Open(fFileName, "Partons");
395 fRL->LoadKinematics();
396 fRL->LoadHeader();
397 } else {
398 fRL = 0x0;
399 }
beac474c 400// Switch off Heavy Flavors on request
401 if (fHFoff) {
402 fPythia->SetMSTP(58, 3);
403 fPythia->SetMSTJ(45, 3);
404 for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
405 }
8d2cd130 406 //
407 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
408
409// Parent and Children Selection
410 switch (fProcess)
411 {
65f2626c 412 case kPyOldUEQ2ordered:
413 case kPyOldUEQ2ordered2:
414 case kPyOldPopcorn:
415 break;
8d2cd130 416 case kPyCharm:
417 case kPyCharmUnforced:
adf4d898 418 case kPyCharmPbPbMNR:
aabc7187 419 case kPyCharmpPbMNR:
e0e89f40 420 case kPyCharmppMNR:
421 case kPyCharmppMNRwmi:
8d2cd130 422 fParentSelect[0] = 411;
423 fParentSelect[1] = 421;
424 fParentSelect[2] = 431;
425 fParentSelect[3] = 4122;
426 fFlavorSelect = 4;
427 break;
adf4d898 428 case kPyD0PbPbMNR:
429 case kPyD0pPbMNR:
430 case kPyD0ppMNR:
8d2cd130 431 fParentSelect[0] = 421;
432 fFlavorSelect = 4;
433 break;
90d7b703 434 case kPyDPlusPbPbMNR:
435 case kPyDPluspPbMNR:
436 case kPyDPlusppMNR:
437 fParentSelect[0] = 411;
438 fFlavorSelect = 4;
439 break;
e0e89f40 440 case kPyDPlusStrangePbPbMNR:
441 case kPyDPlusStrangepPbMNR:
442 case kPyDPlusStrangeppMNR:
443 fParentSelect[0] = 431;
444 fFlavorSelect = 4;
445 break;
8d2cd130 446 case kPyBeauty:
adf4d898 447 case kPyBeautyPbPbMNR:
448 case kPyBeautypPbMNR:
449 case kPyBeautyppMNR:
e0e89f40 450 case kPyBeautyppMNRwmi:
8d2cd130 451 fParentSelect[0]= 511;
452 fParentSelect[1]= 521;
453 fParentSelect[2]= 531;
454 fParentSelect[3]= 5122;
455 fParentSelect[4]= 5132;
456 fParentSelect[5]= 5232;
457 fParentSelect[6]= 5332;
458 fFlavorSelect = 5;
459 break;
460 case kPyBeautyUnforced:
461 fParentSelect[0] = 511;
462 fParentSelect[1] = 521;
463 fParentSelect[2] = 531;
464 fParentSelect[3] = 5122;
465 fParentSelect[4] = 5132;
466 fParentSelect[5] = 5232;
467 fParentSelect[6] = 5332;
468 fFlavorSelect = 5;
469 break;
470 case kPyJpsiChi:
471 case kPyJpsi:
472 fParentSelect[0] = 443;
473 break;
474 case kPyMb:
475 case kPyMbNonDiffr:
d7de4a67 476 case kPyMbMSEL1:
8d2cd130 477 case kPyJets:
478 case kPyDirectGamma:
479 break;
589380c6 480 case kPyW:
0f6ee828 481 case kPyZ:
589380c6 482 break;
8d2cd130 483 }
484//
592f8307 485//
486// JetFinder for Trigger
487//
488// Configure detector (EMCAL like)
489//
d7de4a67 490 fPythia->SetPARU(51, fPycellEtaMax);
491 fPythia->SetMSTU(51, fPycellNEta);
492 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 493//
494// Configure Jet Finder
495//
d7de4a67 496 fPythia->SetPARU(58, fPycellThreshold);
497 fPythia->SetPARU(52, fPycellEtSeed);
498 fPythia->SetPARU(53, fPycellMinEtJet);
499 fPythia->SetPARU(54, fPycellMaxRadius);
500 fPythia->SetMSTU(54, 2);
592f8307 501//
8d2cd130 502// This counts the total number of calls to Pyevnt() per run.
503 fTrialsRun = 0;
504 fQ = 0.;
505 fX1 = 0.;
506 fX2 = 0.;
507 fNev = 0 ;
508//
1d568bc2 509//
510//
8d2cd130 511 AliGenMC::Init();
1d568bc2 512//
513//
514//
515 if (fSetNuclei) {
516 fDyBoost = 0;
517 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
518 }
d7de4a67 519
32d6ef7d 520 if (fQuench) {
5fa4b20b 521 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 522 }
8d2cd130 523}
524
525void AliGenPythia::Generate()
526{
527// Generate one event
e5c87a3d 528
8d2cd130 529 fDecayer->ForceDecay();
530
531 Float_t polar[3] = {0,0,0};
532 Float_t origin[3] = {0,0,0};
a920faf9 533 Float_t p[4];
8d2cd130 534// converts from mm/c to s
535 const Float_t kconv=0.001/2.999792458e8;
536//
537 Int_t nt=0;
538 Int_t jev=0;
539 Int_t j, kf;
540 fTrials=0;
f913ec4f 541 fEventTime = 0.;
542
2590ccf9 543
8d2cd130 544
545 // Set collision vertex position
2590ccf9 546 if (fVertexSmear == kPerEvent) Vertex();
547
8d2cd130 548// event loop
549 while(1)
550 {
32d6ef7d 551//
5fa4b20b 552// Produce event
32d6ef7d 553//
32d6ef7d 554//
5fa4b20b 555// Switch hadronisation off
556//
557 fPythia->SetMSTJ(1, 0);
32d6ef7d 558//
5fa4b20b 559// Either produce new event or read partons from file
560//
561 if (!fReadFromFile) {
beac474c 562 if (!fNewMIS) {
563 fPythia->Pyevnt();
564 } else {
565 fPythia->Pyevnw();
566 }
5fa4b20b 567 fNpartons = fPythia->GetN();
568 } else {
569 printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
570 fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
571 fPythia->SetN(0);
572 LoadEvent(fRL->Stack(), 0 , 1);
573 fPythia->Pyedit(21);
574 }
575
32d6ef7d 576//
577// Run quenching routine
578//
5fa4b20b 579 if (fQuench == 1) {
580 fPythia->Quench();
581 } else if (fQuench == 2){
582 fPythia->Pyquen(208., 0, 0.);
583 }
32d6ef7d 584//
5fa4b20b 585// Switch hadronisation on
32d6ef7d 586//
aea21c57 587 fPythia->SetMSTJ(1, 1);
5fa4b20b 588//
589// .. and perform hadronisation
aea21c57 590// printf("Calling hadronisation %d\n", fPythia->GetN());
5fa4b20b 591 fPythia->Pyexec();
8d2cd130 592 fTrials++;
8d2cd130 593 fPythia->ImportParticles(fParticles,"All");
1d568bc2 594 Boost();
8d2cd130 595//
596//
597//
598 Int_t i;
599
5fa4b20b 600
8d2cd130 601 Int_t np = fParticles->GetEntriesFast();
5fa4b20b 602
7c21f297 603 if (np == 0) continue;
8d2cd130 604//
2590ccf9 605
8d2cd130 606//
607 Int_t* pParent = new Int_t[np];
608 Int_t* pSelected = new Int_t[np];
609 Int_t* trackIt = new Int_t[np];
5fa4b20b 610 for (i = 0; i < np; i++) {
8d2cd130 611 pParent[i] = -1;
612 pSelected[i] = 0;
613 trackIt[i] = 0;
614 }
615
616 Int_t nc = 0; // Total n. of selected particles
617 Int_t nParents = 0; // Selected parents
618 Int_t nTkbles = 0; // Trackable particles
619 if (fProcess != kPyMb && fProcess != kPyJets &&
620 fProcess != kPyDirectGamma &&
589380c6 621 fProcess != kPyMbNonDiffr &&
d7de4a67 622 fProcess != kPyMbMSEL1 &&
e0e89f40 623 fProcess != kPyW && fProcess != kPyZ &&
624 fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
8d2cd130 625
5fa4b20b 626 for (i = 0; i < np; i++) {
2590ccf9 627 TParticle* iparticle = (TParticle *) fParticles->At(i);
8d2cd130 628 Int_t ks = iparticle->GetStatusCode();
629 kf = CheckPDGCode(iparticle->GetPdgCode());
630// No initial state partons
631 if (ks==21) continue;
632//
633// Heavy Flavor Selection
634//
635 // quark ?
636 kf = TMath::Abs(kf);
637 Int_t kfl = kf;
9ff6c04c 638 // Resonance
f913ec4f 639
9ff6c04c 640 if (kfl > 100000) kfl %= 100000;
183a5ca9 641 if (kfl > 10000) kfl %= 10000;
8d2cd130 642 // meson ?
643 if (kfl > 10) kfl/=100;
644 // baryon
645 if (kfl > 10) kfl/=10;
8d2cd130 646 Int_t ipa = iparticle->GetFirstMother()-1;
647 Int_t kfMo = 0;
f913ec4f 648//
649// Establish mother daughter relation between heavy quarks and mesons
650//
651 if (kf >= fFlavorSelect && kf <= 6) {
652 Int_t idau = iparticle->GetFirstDaughter() - 1;
653 if (idau > -1) {
654 TParticle* daughter = (TParticle *) fParticles->At(idau);
655 Int_t pdgD = daughter->GetPdgCode();
656 if (pdgD == 91 || pdgD == 92) {
657 Int_t jmin = daughter->GetFirstDaughter() - 1;
658 Int_t jmax = daughter->GetLastDaughter() - 1;
659 for (Int_t j = jmin; j <= jmax; j++)
660 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
661 } // is string or cluster
662 } // has daughter
663 } // heavy quark
8d2cd130 664
f913ec4f 665
8d2cd130 666 if (ipa > -1) {
667 TParticle * mother = (TParticle *) fParticles->At(ipa);
668 kfMo = TMath::Abs(mother->GetPdgCode());
669 }
f913ec4f 670
8d2cd130 671 // What to keep in Stack?
672 Bool_t flavorOK = kFALSE;
673 Bool_t selectOK = kFALSE;
674 if (fFeedDownOpt) {
32d6ef7d 675 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 676 } else {
32d6ef7d 677 if (kfl > fFlavorSelect) {
678 nc = -1;
679 break;
680 }
681 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 682 }
683 switch (fStackFillOpt) {
684 case kFlavorSelection:
32d6ef7d 685 selectOK = kTRUE;
686 break;
8d2cd130 687 case kParentSelection:
32d6ef7d 688 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
689 break;
8d2cd130 690 }
691 if (flavorOK && selectOK) {
692//
693// Heavy flavor hadron or quark
694//
695// Kinematic seletion on final state heavy flavor mesons
696 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
697 {
9ff6c04c 698 continue;
8d2cd130 699 }
700 pSelected[i] = 1;
701 if (ParentSelected(kf)) ++nParents; // Update parent count
702// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
703 } else {
704// Kinematic seletion on decay products
705 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 706 && !KinematicSelection(iparticle, 1))
8d2cd130 707 {
9ff6c04c 708 continue;
8d2cd130 709 }
710//
711// Decay products
712// Select if mother was selected and is not tracked
713
714 if (pSelected[ipa] &&
715 !trackIt[ipa] && // mother will be tracked ?
716 kfMo != 5 && // mother is b-quark, don't store fragments
717 kfMo != 4 && // mother is c-quark, don't store fragments
718 kf != 92) // don't store string
719 {
720//
721// Semi-stable or de-selected: diselect decay products:
722//
723//
724 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
725 {
726 Int_t ipF = iparticle->GetFirstDaughter();
727 Int_t ipL = iparticle->GetLastDaughter();
728 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
729 }
730// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
731 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
732 }
733 }
734 if (pSelected[i] == -1) pSelected[i] = 0;
735 if (!pSelected[i]) continue;
736 // Count quarks only if you did not include fragmentation
737 if (fFragmentation && kf <= 10) continue;
9ff6c04c 738
8d2cd130 739 nc++;
740// Decision on tracking
741 trackIt[i] = 0;
742//
743// Track final state particle
744 if (ks == 1) trackIt[i] = 1;
745// Track semi-stable particles
d25cfd65 746 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 747// Track particles selected by process if undecayed.
748 if (fForceDecay == kNoDecay) {
749 if (ParentSelected(kf)) trackIt[i] = 1;
750 } else {
751 if (ParentSelected(kf)) trackIt[i] = 0;
752 }
753 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
754//
755//
756
757 } // particle selection loop
758 if (nc > 0) {
759 for (i = 0; i<np; i++) {
760 if (!pSelected[i]) continue;
761 TParticle * iparticle = (TParticle *) fParticles->At(i);
762 kf = CheckPDGCode(iparticle->GetPdgCode());
763 Int_t ks = iparticle->GetStatusCode();
764 p[0] = iparticle->Px();
765 p[1] = iparticle->Py();
766 p[2] = iparticle->Pz();
a920faf9 767 p[3] = iparticle->Energy();
768
2590ccf9 769 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
770 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
771 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
772
8d2cd130 773 Float_t tof = kconv*iparticle->T();
774 Int_t ipa = iparticle->GetFirstMother()-1;
775 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 776
777 PushTrack(fTrackIt*trackIt[i], iparent, kf,
778 p[0], p[1], p[2], p[3],
779 origin[0], origin[1], origin[2], tof,
780 polar[0], polar[1], polar[2],
781 kPPrimary, nt, 1., ks);
8d2cd130 782 pParent[i] = nt;
783 KeepTrack(nt);
642f15cf 784 } // PushTrack loop
8d2cd130 785 }
786 } else {
787 nc = GenerateMB();
788 } // mb ?
f913ec4f 789
790 GetSubEventTime();
8d2cd130 791
234f6d32 792 delete[] pParent;
793 delete[] pSelected;
794 delete[] trackIt;
8d2cd130 795
796 if (nc > 0) {
797 switch (fCountMode) {
798 case kCountAll:
799 // printf(" Count all \n");
800 jev += nc;
801 break;
802 case kCountParents:
803 // printf(" Count parents \n");
804 jev += nParents;
805 break;
806 case kCountTrackables:
807 // printf(" Count trackable \n");
808 jev += nTkbles;
809 break;
810 }
811 if (jev >= fNpart || fNpart == -1) {
812 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 813
8d2cd130 814 fQ += fPythia->GetVINT(51);
815 fX1 += fPythia->GetVINT(41);
816 fX2 += fPythia->GetVINT(42);
817 fTrialsRun += fTrials;
818 fNev++;
819 MakeHeader();
820 break;
821 }
822 }
823 } // event loop
824 SetHighWaterMark(nt);
825// adjust weight due to kinematic selection
b88f5cea 826// AdjustWeights();
8d2cd130 827// get cross-section
828 fXsection=fPythia->GetPARI(1);
829}
830
831Int_t AliGenPythia::GenerateMB()
832{
833//
834// Min Bias selection and other global selections
835//
836 Int_t i, kf, nt, iparent;
837 Int_t nc = 0;
bf950da8 838 Float_t p[4];
8d2cd130 839 Float_t polar[3] = {0,0,0};
840 Float_t origin[3] = {0,0,0};
841// converts from mm/c to s
842 const Float_t kconv=0.001/2.999792458e8;
843
e0e89f40 844
845
5fa4b20b 846 Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
aea21c57 847
5fa4b20b 848
e0e89f40 849
8d2cd130 850 Int_t* pParent = new Int_t[np];
851 for (i=0; i< np; i++) pParent[i] = -1;
852 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
853 TParticle* jet1 = (TParticle *) fParticles->At(6);
854 TParticle* jet2 = (TParticle *) fParticles->At(7);
234f6d32 855 if (!CheckTrigger(jet1, jet2)) {
856 delete [] pParent;
857 return 0;
858 }
8d2cd130 859 }
e0e89f40 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
8d2cd130 1186
1187AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
1188{
1189// Assignment operator
014a9521 1190 rhs.Copy(*this);
8d2cd130 1191 return *this;
1192}
1193
5fa4b20b 1194void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1195{
1196//
1197// Load event into Pythia Common Block
1198//
5fa4b20b 1199
32d6ef7d 1200 Int_t npart = stack -> GetNprimary();
1201 Int_t n0 = 0;
1202
1203 if (!flag) {
1204 (fPythia->GetPyjets())->N = npart;
1205 } else {
1206 n0 = (fPythia->GetPyjets())->N;
1207 (fPythia->GetPyjets())->N = n0 + npart;
1208 }
1209
1210
8d2cd130 1211 for (Int_t part = 0; part < npart; part++) {
7184e472 1212 TParticle *mPart = stack->Particle(part);
32d6ef7d 1213
7184e472 1214 Int_t kf = mPart->GetPdgCode();
1215 Int_t ks = mPart->GetStatusCode();
1216 Int_t idf = mPart->GetFirstDaughter();
1217 Int_t idl = mPart->GetLastDaughter();
5fa4b20b 1218
1219 if (reHadr) {
1220 if (ks == 11 || ks == 12) {
1221 ks -= 10;
1222 idf = -1;
1223 idl = -1;
1224 }
1225 }
32d6ef7d 1226
7184e472 1227 Float_t px = mPart->Px();
1228 Float_t py = mPart->Py();
1229 Float_t pz = mPart->Pz();
1230 Float_t e = mPart->Energy();
1231 Float_t m = mPart->GetCalcMass();
8d2cd130 1232
1233
32d6ef7d 1234 (fPythia->GetPyjets())->P[0][part+n0] = px;
1235 (fPythia->GetPyjets())->P[1][part+n0] = py;
1236 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1237 (fPythia->GetPyjets())->P[3][part+n0] = e;
1238 (fPythia->GetPyjets())->P[4][part+n0] = m;
8d2cd130 1239
32d6ef7d 1240 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1241 (fPythia->GetPyjets())->K[0][part+n0] = ks;
5fa4b20b 1242 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1243 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
7184e472 1244 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
8d2cd130 1245 }
1246}
1247
5fa4b20b 1248
014a9521 1249void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1250{
1251//
1252// Calls the Pythia jet finding algorithm to find jets in the current event
1253//
1254//
8d2cd130 1255//
1256// Save jets
1257 Int_t n = fPythia->GetN();
1258
1259//
1260// Run Jet Finder
1261 fPythia->Pycell(njets);
1262 Int_t i;
1263 for (i = 0; i < njets; i++) {
1264 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1265 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1266 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1267 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1268
1269 jets[0][i] = px;
1270 jets[1][i] = py;
1271 jets[2][i] = pz;
1272 jets[3][i] = e;
1273 }
1274}
1275
1276
1277
1278void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1279{
1280//
1281// Calls the Pythia clustering algorithm to find jets in the current event
1282//
1283 Int_t n = fPythia->GetN();
1284 nJets = 0;
1285 nJetsTrig = 0;
1286 if (fJetReconstruction == kCluster) {
1287//
1288// Configure cluster algorithm
1289//
1290 fPythia->SetPARU(43, 2.);
1291 fPythia->SetMSTU(41, 1);
1292//
1293// Call cluster algorithm
1294//
1295 fPythia->Pyclus(nJets);
1296//
1297// Loading jets from common block
1298//
1299 } else {
592f8307 1300
8d2cd130 1301//
1302// Run Jet Finder
1303 fPythia->Pycell(nJets);
1304 }
1305
1306 Int_t i;
1307 for (i = 0; i < nJets; i++) {
1308 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1309 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1310 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1311 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1312 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1313 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1314 Float_t theta = TMath::ATan2(pt,pz);
1315 Float_t et = e * TMath::Sin(theta);
1316 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1317 if (
1318 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1319 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1320 et > fEtMinJet && et < fEtMaxJet
1321 )
1322 {
1323 jets[0][nJetsTrig] = px;
1324 jets[1][nJetsTrig] = py;
1325 jets[2][nJetsTrig] = pz;
1326 jets[3][nJetsTrig] = e;
1327 nJetsTrig++;
5fa4b20b 1328// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1329 } else {
1330// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1331 }
1332 }
1333}
1334
f913ec4f 1335void AliGenPythia::GetSubEventTime()
1336{
1337 // Calculates time of the next subevent
9d7108a7 1338 fEventTime = 0.;
1339 if (fEventsTime) {
1340 TArrayF &array = *fEventsTime;
1341 fEventTime = array[fCurSubEvent++];
1342 }
1343 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1344 return;
f913ec4f 1345}
8d2cd130 1346
1347#ifdef never
1348void AliGenPythia::Streamer(TBuffer &R__b)
1349{
1350 // Stream an object of class AliGenPythia.
1351
1352 if (R__b.IsReading()) {
1353 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1354 AliGenerator::Streamer(R__b);
1355 R__b >> (Int_t&)fProcess;
1356 R__b >> (Int_t&)fStrucFunc;
1357 R__b >> (Int_t&)fForceDecay;
1358 R__b >> fEnergyCMS;
1359 R__b >> fKineBias;
1360 R__b >> fTrials;
1361 fParentSelect.Streamer(R__b);
1362 fChildSelect.Streamer(R__b);
1363 R__b >> fXsection;
1364// (AliPythia::Instance())->Streamer(R__b);
1365 R__b >> fPtHardMin;
1366 R__b >> fPtHardMax;
1367// if (fDecayer) fDecayer->Streamer(R__b);
1368 } else {
1369 R__b.WriteVersion(AliGenPythia::IsA());
1370 AliGenerator::Streamer(R__b);
1371 R__b << (Int_t)fProcess;
1372 R__b << (Int_t)fStrucFunc;
1373 R__b << (Int_t)fForceDecay;
1374 R__b << fEnergyCMS;
1375 R__b << fKineBias;
1376 R__b << fTrials;
1377 fParentSelect.Streamer(R__b);
1378 fChildSelect.Streamer(R__b);
1379 R__b << fXsection;
1380// R__b << fPythia;
1381 R__b << fPtHardMin;
1382 R__b << fPtHardMax;
1383 // fDecayer->Streamer(R__b);
1384 }
1385}
1386#endif
1387
90d7b703 1388
589380c6 1389