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