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