]> 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
Correction.
[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:
04081a91 447 case kPyMbWithDirectPhoton:
8d2cd130 448 case kPyMbNonDiffr:
d7de4a67 449 case kPyMbMSEL1:
8d2cd130 450 case kPyJets:
451 case kPyDirectGamma:
e33acb67 452 case kPyLhwgMb:
8d2cd130 453 break;
589380c6 454 case kPyW:
0f6ee828 455 case kPyZ:
589380c6 456 break;
8d2cd130 457 }
458//
592f8307 459//
460// JetFinder for Trigger
461//
462// Configure detector (EMCAL like)
463//
d7de4a67 464 fPythia->SetPARU(51, fPycellEtaMax);
465 fPythia->SetMSTU(51, fPycellNEta);
466 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 467//
468// Configure Jet Finder
469//
d7de4a67 470 fPythia->SetPARU(58, fPycellThreshold);
471 fPythia->SetPARU(52, fPycellEtSeed);
472 fPythia->SetPARU(53, fPycellMinEtJet);
473 fPythia->SetPARU(54, fPycellMaxRadius);
474 fPythia->SetMSTU(54, 2);
592f8307 475//
8d2cd130 476// This counts the total number of calls to Pyevnt() per run.
477 fTrialsRun = 0;
478 fQ = 0.;
479 fX1 = 0.;
480 fX2 = 0.;
481 fNev = 0 ;
482//
1d568bc2 483//
484//
8d2cd130 485 AliGenMC::Init();
1d568bc2 486//
487//
488//
489 if (fSetNuclei) {
490 fDyBoost = 0;
491 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
492 }
d7de4a67 493
32d6ef7d 494 if (fQuench) {
5fa4b20b 495 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 496 }
3a709cfa 497 fPythia->SetPARJ(200, 0.0);
498
b976f7d8 499 if (fQuench == 3) {
500 // Nestor's change of the splittings
501 fPythia->SetPARJ(200, 0.8);
502 fPythia->SetMSTJ(41, 1); // QCD radiation only
503 fPythia->SetMSTJ(42, 2); // angular ordering
504 fPythia->SetMSTJ(44, 2); // option to run alpha_s
505 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
506 fPythia->SetMSTJ(50, 0); // No coherence in first branching
507 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
508 }
8d2cd130 509}
510
511void AliGenPythia::Generate()
512{
513// Generate one event
19ef8cf0 514 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
8d2cd130 515 fDecayer->ForceDecay();
516
517 Float_t polar[3] = {0,0,0};
518 Float_t origin[3] = {0,0,0};
a920faf9 519 Float_t p[4];
8d2cd130 520// converts from mm/c to s
521 const Float_t kconv=0.001/2.999792458e8;
522//
523 Int_t nt=0;
524 Int_t jev=0;
525 Int_t j, kf;
526 fTrials=0;
f913ec4f 527 fEventTime = 0.;
528
2590ccf9 529
8d2cd130 530
531 // Set collision vertex position
2590ccf9 532 if (fVertexSmear == kPerEvent) Vertex();
533
8d2cd130 534// event loop
535 while(1)
536 {
32d6ef7d 537//
5fa4b20b 538// Produce event
32d6ef7d 539//
32d6ef7d 540//
5fa4b20b 541// Switch hadronisation off
542//
543 fPythia->SetMSTJ(1, 0);
32d6ef7d 544//
5fa4b20b 545// Either produce new event or read partons from file
546//
547 if (!fReadFromFile) {
beac474c 548 if (!fNewMIS) {
549 fPythia->Pyevnt();
550 } else {
551 fPythia->Pyevnw();
552 }
5fa4b20b 553 fNpartons = fPythia->GetN();
554 } else {
555 printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
556 fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
557 fPythia->SetN(0);
558 LoadEvent(fRL->Stack(), 0 , 1);
559 fPythia->Pyedit(21);
560 }
561
32d6ef7d 562//
563// Run quenching routine
564//
5fa4b20b 565 if (fQuench == 1) {
566 fPythia->Quench();
567 } else if (fQuench == 2){
568 fPythia->Pyquen(208., 0, 0.);
b976f7d8 569 } else if (fQuench == 3) {
570 // Quenching is via multiplicative correction of the splittings
5fa4b20b 571 }
b976f7d8 572
32d6ef7d 573//
5fa4b20b 574// Switch hadronisation on
32d6ef7d 575//
20e47f08 576 if (fHadronisation) {
577 fPythia->SetMSTJ(1, 1);
5fa4b20b 578//
579// .. and perform hadronisation
aea21c57 580// printf("Calling hadronisation %d\n", fPythia->GetN());
20e47f08 581 fPythia->Pyexec();
582 }
8d2cd130 583 fTrials++;
8507138f 584 fPythia->ImportParticles(&fParticles,"All");
1d568bc2 585 Boost();
8d2cd130 586//
587//
588//
589 Int_t i;
590
dbd64db6 591 fNprimaries = 0;
8507138f 592 Int_t np = fParticles.GetEntriesFast();
5fa4b20b 593
7c21f297 594 if (np == 0) continue;
8d2cd130 595//
2590ccf9 596
8d2cd130 597//
598 Int_t* pParent = new Int_t[np];
599 Int_t* pSelected = new Int_t[np];
600 Int_t* trackIt = new Int_t[np];
5fa4b20b 601 for (i = 0; i < np; i++) {
8d2cd130 602 pParent[i] = -1;
603 pSelected[i] = 0;
604 trackIt[i] = 0;
605 }
606
607 Int_t nc = 0; // Total n. of selected particles
608 Int_t nParents = 0; // Selected parents
609 Int_t nTkbles = 0; // Trackable particles
f529e69b 610 if (fProcess != kPyMbDefault &&
611 fProcess != kPyMb &&
04081a91 612 fProcess != kPyMbWithDirectPhoton &&
f529e69b 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
82e0bdff 865 else if (fPi0InCalo) pdg = 111 ; // Pi0
ec2c406e 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 ||
82e0bdff 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
82e0bdff 877 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
9fd8e520 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 ||
04081a91 1002 fProcess == kPyMbWithDirectPhoton ||
f529e69b 1003 fProcess == kPyMbNonDiffr)
0f6ee828 1004 && (fCutOnChild == 1) ) {
1005 if ( !CheckKinematicsOnChild() ) {
234f6d32 1006 delete[] pParent;
0f6ee828 1007 return 0;
1008 }
aea21c57 1009 }
1010
1011
f913ec4f 1012 for (i = 0; i < np; i++) {
8d2cd130 1013 Int_t trackIt = 0;
8507138f 1014 TParticle * iparticle = (TParticle *) fParticles.At(i);
8d2cd130 1015 kf = CheckPDGCode(iparticle->GetPdgCode());
1016 Int_t ks = iparticle->GetStatusCode();
1017 Int_t km = iparticle->GetFirstMother();
1018 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
1019 (ks != 1) ||
1020 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
1021 nc++;
1022 if (ks == 1) trackIt = 1;
1023 Int_t ipa = iparticle->GetFirstMother()-1;
1024
1025 iparent = (ipa > -1) ? pParent[ipa] : -1;
1026
1027//
1028// store track information
1029 p[0] = iparticle->Px();
1030 p[1] = iparticle->Py();
1031 p[2] = iparticle->Pz();
a920faf9 1032 p[3] = iparticle->Energy();
1406f599 1033
a920faf9 1034
2590ccf9 1035 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1036 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1037 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1038
f913ec4f 1039 Float_t tof = fEventTime + kconv * iparticle->T();
a920faf9 1040
1041 PushTrack(fTrackIt*trackIt, iparent, kf,
1042 p[0], p[1], p[2], p[3],
1043 origin[0], origin[1], origin[2], tof,
1044 polar[0], polar[1], polar[2],
1045 kPPrimary, nt, 1., ks);
dbd64db6 1046 fNprimaries++;
8d2cd130 1047 KeepTrack(nt);
1048 pParent[i] = nt;
f913ec4f 1049 SetHighWaterMark(nt);
1050
8d2cd130 1051 } // select particle
1052 } // particle loop
1053
234f6d32 1054 delete[] pParent;
e0e89f40 1055
f913ec4f 1056 return 1;
8d2cd130 1057}
1058
1059
1060void AliGenPythia::FinishRun()
1061{
1062// Print x-section summary
1063 fPythia->Pystat(1);
2779fc64 1064
1065 if (fNev > 0.) {
1066 fQ /= fNev;
1067 fX1 /= fNev;
1068 fX2 /= fNev;
1069 }
1070
8d2cd130 1071 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1072 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 1073}
1074
7184e472 1075void AliGenPythia::AdjustWeights() const
8d2cd130 1076{
1077// Adjust the weights after generation of all events
1078//
e2bddf81 1079 if (gAlice) {
1080 TParticle *part;
1081 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1082 for (Int_t i=0; i<ntrack; i++) {
1083 part= gAlice->GetMCApp()->Particle(i);
1084 part->SetWeight(part->GetWeight()*fKineBias);
1085 }
8d2cd130 1086 }
1087}
1088
20e47f08 1089void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
8d2cd130 1090{
1091// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 1092
1a626d4e 1093 fAProjectile = a1;
1094 fATarget = a2;
20e47f08 1095 fNucPdf = pdfset; // 0 EKS98 1 EPS08
1d568bc2 1096 fSetNuclei = kTRUE;
8d2cd130 1097}
1098
1099
1100void AliGenPythia::MakeHeader()
1101{
7184e472 1102//
1103// Make header for the simulated event
1104//
183a5ca9 1105 if (gAlice) {
1106 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 1107 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 1108 }
1109
8d2cd130 1110// Builds the event header, to be called after each event
e5c87a3d 1111 if (fHeader) delete fHeader;
1112 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 1113//
1114// Event type
e5c87a3d 1115 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 1116//
1117// Number of trials
e5c87a3d 1118 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 1119//
1120// Event Vertex
d25cfd65 1121 fHeader->SetPrimaryVertex(fVertex);
dbd64db6 1122
1123//
1124// Number of primaries
1125 fHeader->SetNProduced(fNprimaries);
8d2cd130 1126//
1127// Jets that have triggered
f913ec4f 1128
19ef8cf0 1129 if (fProcess == kPyJets || fProcess == kPyDirectGamma)
8d2cd130 1130 {
1131 Int_t ntrig, njet;
1132 Float_t jets[4][10];
1133 GetJets(njet, ntrig, jets);
9ff6c04c 1134
8d2cd130 1135
1136 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 1137 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 1138 jets[3][i]);
1139 }
1140 }
5fa4b20b 1141//
1142// Copy relevant information from external header, if present.
1143//
1144 Float_t uqJet[4];
1145
1146 if (fRL) {
1147 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1148 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1149 {
1150 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1151
1152
1153 exHeader->TriggerJet(i, uqJet);
1154 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1155 }
1156 }
1157//
1158// Store quenching parameters
1159//
1160 if (fQuench){
1161 Double_t z[4];
1162 Double_t xp, yp;
7c21f297 1163 if (fQuench == 1) {
1164 // Pythia::Quench()
1165 fPythia->GetQuenchingParameters(xp, yp, z);
1166 } else {
1167 // Pyquen
1168 Double_t r1 = PARIMP.rb1;
1169 Double_t r2 = PARIMP.rb2;
1170 Double_t b = PARIMP.b1;
1171 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1172 Double_t phi = PARIMP.psib1;
1173 xp = r * TMath::Cos(phi);
1174 yp = r * TMath::Sin(phi);
1175
1176 }
1177 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1178 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1179 }
beac474c 1180//
1181// Store pt^hard
1182 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 1183//
cf57b268 1184// Pass header
5fa4b20b 1185//
cf57b268 1186 AddHeader(fHeader);
4c4eac97 1187 fHeader = 0x0;
8d2cd130 1188}
cf57b268 1189
8d2cd130 1190Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1191{
1192// Check the kinematic trigger condition
1193//
1194 Double_t eta[2];
1195 eta[0] = jet1->Eta();
1196 eta[1] = jet2->Eta();
1197 Double_t phi[2];
1198 phi[0] = jet1->Phi();
1199 phi[1] = jet2->Phi();
1200 Int_t pdg[2];
1201 pdg[0] = jet1->GetPdgCode();
1202 pdg[1] = jet2->GetPdgCode();
1203 Bool_t triggered = kFALSE;
1204
1205 if (fProcess == kPyJets) {
1206 Int_t njets = 0;
1207 Int_t ntrig = 0;
1208 Float_t jets[4][10];
1209//
1210// Use Pythia clustering on parton level to determine jet axis
1211//
1212 GetJets(njets, ntrig, jets);
1213
76d6ba9a 1214 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1215//
1216 } else {
1217 Int_t ij = 0;
1218 Int_t ig = 1;
1219 if (pdg[0] == kGamma) {
1220 ij = 1;
1221 ig = 0;
1222 }
1223 //Check eta range first...
1224 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1225 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1226 {
1227 //Eta is okay, now check phi range
1228 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1229 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1230 {
1231 triggered = kTRUE;
1232 }
1233 }
1234 }
1235 return triggered;
1236}
aea21c57 1237
1238
aea21c57 1239
7184e472 1240Bool_t AliGenPythia::CheckKinematicsOnChild(){
1241//
1242//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1243//
aea21c57 1244 Bool_t checking = kFALSE;
1245 Int_t j, kcode, ks, km;
1246 Int_t nPartAcc = 0; //number of particles in the acceptance range
1247 Int_t numberOfAcceptedParticles = 1;
1248 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1249 Int_t npart = fParticles.GetEntriesFast();
aea21c57 1250
0f6ee828 1251 for (j = 0; j<npart; j++) {
8507138f 1252 TParticle * jparticle = (TParticle *) fParticles.At(j);
aea21c57 1253 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1254 ks = jparticle->GetStatusCode();
1255 km = jparticle->GetFirstMother();
1256
1257 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1258 nPartAcc++;
1259 }
0f6ee828 1260 if( numberOfAcceptedParticles <= nPartAcc){
1261 checking = kTRUE;
1262 break;
1263 }
aea21c57 1264 }
0f6ee828 1265
aea21c57 1266 return checking;
aea21c57 1267}
1268
5fa4b20b 1269void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1270{
1271//
1272// Load event into Pythia Common Block
1273//
5fa4b20b 1274
32d6ef7d 1275 Int_t npart = stack -> GetNprimary();
1276 Int_t n0 = 0;
1277
1278 if (!flag) {
1279 (fPythia->GetPyjets())->N = npart;
1280 } else {
1281 n0 = (fPythia->GetPyjets())->N;
1282 (fPythia->GetPyjets())->N = n0 + npart;
1283 }
1284
1285
8d2cd130 1286 for (Int_t part = 0; part < npart; part++) {
7184e472 1287 TParticle *mPart = stack->Particle(part);
32d6ef7d 1288
7184e472 1289 Int_t kf = mPart->GetPdgCode();
1290 Int_t ks = mPart->GetStatusCode();
1291 Int_t idf = mPart->GetFirstDaughter();
1292 Int_t idl = mPart->GetLastDaughter();
5fa4b20b 1293
1294 if (reHadr) {
1295 if (ks == 11 || ks == 12) {
1296 ks -= 10;
1297 idf = -1;
1298 idl = -1;
1299 }
1300 }
32d6ef7d 1301
7184e472 1302 Float_t px = mPart->Px();
1303 Float_t py = mPart->Py();
1304 Float_t pz = mPart->Pz();
1305 Float_t e = mPart->Energy();
1306 Float_t m = mPart->GetCalcMass();
8d2cd130 1307
1308
32d6ef7d 1309 (fPythia->GetPyjets())->P[0][part+n0] = px;
1310 (fPythia->GetPyjets())->P[1][part+n0] = py;
1311 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1312 (fPythia->GetPyjets())->P[3][part+n0] = e;
1313 (fPythia->GetPyjets())->P[4][part+n0] = m;
8d2cd130 1314
32d6ef7d 1315 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1316 (fPythia->GetPyjets())->K[0][part+n0] = ks;
5fa4b20b 1317 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1318 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
7184e472 1319 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
8d2cd130 1320 }
1321}
1322
5fa4b20b 1323
014a9521 1324void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1325{
1326//
1327// Calls the Pythia jet finding algorithm to find jets in the current event
1328//
1329//
8d2cd130 1330//
1331// Save jets
1332 Int_t n = fPythia->GetN();
1333
1334//
1335// Run Jet Finder
1336 fPythia->Pycell(njets);
1337 Int_t i;
1338 for (i = 0; i < njets; i++) {
1339 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1340 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1341 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1342 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1343
1344 jets[0][i] = px;
1345 jets[1][i] = py;
1346 jets[2][i] = pz;
1347 jets[3][i] = e;
1348 }
1349}
1350
1351
1352
1353void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1354{
1355//
1356// Calls the Pythia clustering algorithm to find jets in the current event
1357//
1358 Int_t n = fPythia->GetN();
1359 nJets = 0;
1360 nJetsTrig = 0;
1361 if (fJetReconstruction == kCluster) {
1362//
1363// Configure cluster algorithm
1364//
1365 fPythia->SetPARU(43, 2.);
1366 fPythia->SetMSTU(41, 1);
1367//
1368// Call cluster algorithm
1369//
1370 fPythia->Pyclus(nJets);
1371//
1372// Loading jets from common block
1373//
1374 } else {
592f8307 1375
8d2cd130 1376//
1377// Run Jet Finder
1378 fPythia->Pycell(nJets);
1379 }
1380
1381 Int_t i;
1382 for (i = 0; i < nJets; i++) {
1383 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1384 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1385 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1386 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1387 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1388 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1389 Float_t theta = TMath::ATan2(pt,pz);
1390 Float_t et = e * TMath::Sin(theta);
1391 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1392 if (
1393 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1394 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1395 et > fEtMinJet && et < fEtMaxJet
1396 )
1397 {
1398 jets[0][nJetsTrig] = px;
1399 jets[1][nJetsTrig] = py;
1400 jets[2][nJetsTrig] = pz;
1401 jets[3][nJetsTrig] = e;
1402 nJetsTrig++;
5fa4b20b 1403// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1404 } else {
1405// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1406 }
1407 }
1408}
1409
f913ec4f 1410void AliGenPythia::GetSubEventTime()
1411{
1412 // Calculates time of the next subevent
9d7108a7 1413 fEventTime = 0.;
1414 if (fEventsTime) {
1415 TArrayF &array = *fEventsTime;
1416 fEventTime = array[fCurSubEvent++];
1417 }
1418 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1419 return;
f913ec4f 1420}
8d2cd130 1421
ec2c406e 1422Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1423{
1424 // Is particle in EMCAL acceptance?
ec2c406e 1425 // phi in degrees, etamin=-etamax
9fd8e520 1426 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1427 eta < fEMCALEta )
ec2c406e 1428 return kTRUE;
1429 else
1430 return kFALSE;
1431}
1432
1433Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1434{
1435 // Is particle in PHOS acceptance?
1436 // Acceptance slightly larger considered.
1437 // phi in degrees, etamin=-etamax
9fd8e520 1438 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1439 eta < fPHOSEta )
ec2c406e 1440 return kTRUE;
1441 else
1442 return kFALSE;
1443}
1444
90a236ce 1445void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1446{
1447 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1448 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1449 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1450 Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1451
1452 //calculate deltaphi
8507138f 1453 TParticle* ph = (TParticle *) fParticles.At(iphcand);
90a236ce 1454 Double_t phphi = ph->Phi();
1455 Double_t deltaphi = phiPHOS - phphi;
1456
1457
1458
1459 //loop for all particles and produce the phi rotation
8507138f 1460 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
90a236ce 1461 Double_t oldphi, newphi;
1462 Double_t newVx, newVy, R, Vz, time;
1463 Double_t newPx, newPy, pt, Pz, e;
1464 for(Int_t i=0; i< np; i++) {
8507138f 1465 TParticle* iparticle = (TParticle *) fParticles.At(i);
90a236ce 1466 oldphi = iparticle->Phi();
1467 newphi = oldphi + deltaphi;
1468 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1469 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1470
1471 R = iparticle->R();
1472 newVx = R*TMath::Cos(newphi);
1473 newVy = R*TMath::Sin(newphi);
1474 Vz = iparticle->Vz(); // don't transform
1475 time = iparticle->T(); // don't transform
1476
1477 pt = iparticle->Pt();
1478 newPx = pt*TMath::Cos(newphi);
1479 newPy = pt*TMath::Sin(newphi);
1480 Pz = iparticle->Pz(); // don't transform
1481 e = iparticle->Energy(); // don't transform
1482
1483 // apply rotation
1484 iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1485 iparticle->SetMomentum(newPx, newPy, Pz, e);
1486
1487 } //end particle loop
1488
1489 // now let's check that we put correctly the candidate photon in PHOS
1490 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1491 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1492 if(IsInPHOS(phi,eta))
1493 okdd = kTRUE;
1494}
ec2c406e 1495
1496
8d2cd130 1497#ifdef never
1498void AliGenPythia::Streamer(TBuffer &R__b)
1499{
1500 // Stream an object of class AliGenPythia.
1501
1502 if (R__b.IsReading()) {
1503 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1504 AliGenerator::Streamer(R__b);
1505 R__b >> (Int_t&)fProcess;
1506 R__b >> (Int_t&)fStrucFunc;
1507 R__b >> (Int_t&)fForceDecay;
1508 R__b >> fEnergyCMS;
1509 R__b >> fKineBias;
1510 R__b >> fTrials;
1511 fParentSelect.Streamer(R__b);
1512 fChildSelect.Streamer(R__b);
1513 R__b >> fXsection;
1514// (AliPythia::Instance())->Streamer(R__b);
1515 R__b >> fPtHardMin;
1516 R__b >> fPtHardMax;
1517// if (fDecayer) fDecayer->Streamer(R__b);
1518 } else {
1519 R__b.WriteVersion(AliGenPythia::IsA());
1520 AliGenerator::Streamer(R__b);
1521 R__b << (Int_t)fProcess;
1522 R__b << (Int_t)fStrucFunc;
1523 R__b << (Int_t)fForceDecay;
1524 R__b << fEnergyCMS;
1525 R__b << fKineBias;
1526 R__b << fTrials;
1527 fParentSelect.Streamer(R__b);
1528 fChildSelect.Streamer(R__b);
1529 R__b << fXsection;
1530// R__b << fPythia;
1531 R__b << fPtHardMin;
1532 R__b << fPtHardMax;
1533 // fDecayer->Streamer(R__b);
1534 }
1535}
1536#endif
1537
90d7b703 1538
589380c6 1539