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