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