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