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