]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA6/AliGenPythiaPlus.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
1) decouple the rapidity cut used in AliGenPythia to select the events with c or...
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythiaPlus.cxx
CommitLineData
39d810c8 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
16/* $Id$ */
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
28#include <TClonesArray.h>
29#include <TDatabasePDG.h>
30#include <TParticle.h>
31#include <TPDGCode.h>
32#include <TSystem.h>
33#include <TTree.h>
34#include "AliConst.h"
35#include "AliDecayerPythia.h"
36#include "AliGenPythiaPlus.h"
37#include "AliHeader.h"
38#include "AliGenPythiaEventHeader.h"
39#include "AliPythiaBase.h"
40#include "AliPythiaRndm.h"
41#include "AliRun.h"
42#include "AliStack.h"
43#include "AliRunLoader.h"
44#include "AliMC.h"
45#include "PyquenCommon.h"
46
47ClassImp(AliGenPythiaPlus)
48
49
50AliGenPythiaPlus::AliGenPythiaPlus():
51 AliGenMC(),
52 fPythia(0),
53 fProcess(kPyCharm),
54 fStrucFunc(kCTEQ5L),
39d810c8 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 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()),
d8850d4f 95 fUseYCutHQ(kFALSE),
96 fYMinHQ(-20.),
97 fYMaxHQ(20.),
39d810c8 98 fPycellEtaMax(2.),
99 fPycellNEta(274),
100 fPycellNPhi(432),
101 fPycellThreshold(0.),
102 fPycellEtSeed(4.),
103 fPycellMinEtJet(10.),
104 fPycellMaxRadius(1.),
105 fStackFillOpt(kFlavorSelection),
106 fFeedDownOpt(kTRUE),
107 fFragmentation(kTRUE),
108 fSetNuclei(kFALSE),
109 fNewMIS(kFALSE),
110 fHFoff(kFALSE),
111 fTriggerParticle(0),
112 fTriggerEta(0.9),
113 fCountMode(kCountAll),
114 fHeader(0),
115 fRL(0),
116 fFileName(0),
117 fFragPhotonInCalo(kFALSE),
118 fPi0InCalo(kFALSE) ,
119 fPhotonInCalo(kFALSE),
120 fCheckEMCAL(kFALSE),
121 fCheckPHOS(kFALSE),
122 fCheckPHOSeta(kFALSE),
123 fFragPhotonOrPi0MinPt(0),
124 fPhotonMinPt(0),
125 fPHOSMinPhi(219.),
126 fPHOSMaxPhi(321.),
127 fPHOSEta(0.13),
128 fEMCALMinPhi(79.),
129 fEMCALMaxPhi(191.),
75d4f39e 130 fEMCALEta(0.71),
bffc134e 131 fItune(-1),
132 fInfo(1)
39d810c8 133{
134// Default Constructor
e7c989e4 135 fEnergyCMS = 5500.;
39d810c8 136 SetNuclei(0,0);
137 if (!AliPythiaRndm::GetPythiaRandom())
138 AliPythiaRndm::SetPythiaRandom(GetRandom());
139}
140
141AliGenPythiaPlus::AliGenPythiaPlus(AliPythiaBase* pythia)
142 :AliGenMC(-1),
143 fPythia(pythia),
144 fProcess(kPyCharm),
145 fStrucFunc(kCTEQ5L),
39d810c8 146 fKineBias(0.),
147 fTrials(0),
148 fTrialsRun(0),
149 fQ(0.),
150 fX1(0.),
151 fX2(0.),
152 fEventTime(0.),
153 fInteractionRate(0.),
154 fTimeWindow(0.),
155 fCurSubEvent(0),
156 fEventsTime(0),
157 fNev(0),
158 fFlavorSelect(0),
159 fXsection(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()),
d8850d4f 186 fUseYCutHQ(kFALSE),
187 fYMinHQ(-20.),
188 fYMaxHQ(20.),
39d810c8 189 fPycellEtaMax(2.),
190 fPycellNEta(274),
191 fPycellNPhi(432),
192 fPycellThreshold(0.),
193 fPycellEtSeed(4.),
194 fPycellMinEtJet(10.),
195 fPycellMaxRadius(1.),
196 fStackFillOpt(kFlavorSelection),
197 fFeedDownOpt(kTRUE),
198 fFragmentation(kTRUE),
199 fSetNuclei(kFALSE),
200 fNewMIS(kFALSE),
201 fHFoff(kFALSE),
202 fTriggerParticle(0),
203 fTriggerEta(0.9),
204 fCountMode(kCountAll),
205 fHeader(0),
206 fRL(0),
207 fFileName(0),
208 fFragPhotonInCalo(kFALSE),
209 fPi0InCalo(kFALSE) ,
210 fPhotonInCalo(kFALSE),
211 fCheckEMCAL(kFALSE),
212 fCheckPHOS(kFALSE),
213 fCheckPHOSeta(kFALSE),
214 fFragPhotonOrPi0MinPt(0),
215 fPhotonMinPt(0),
216 fPHOSMinPhi(219.),
217 fPHOSMaxPhi(321.),
218 fPHOSEta(0.13),
219 fEMCALMinPhi(79.),
220 fEMCALMaxPhi(191.),
75d4f39e 221 fEMCALEta(0.71),
bffc134e 222 fItune(-1),
223 fInfo(1)
39d810c8 224{
225// default charm production at 5. 5 TeV
226// semimuonic decay
227// structure function GRVHO
228//
e7c989e4 229 fEnergyCMS = 5500.;
39d810c8 230 fName = "Pythia";
231 fTitle= "Particle Generator using PYTHIA";
232 SetForceDecay();
233 // Set random number generator
234 if (!AliPythiaRndm::GetPythiaRandom())
235 AliPythiaRndm::SetPythiaRandom(GetRandom());
39d810c8 236 SetNuclei(0,0);
237 }
238
239AliGenPythiaPlus::~AliGenPythiaPlus()
240{
241// Destructor
242 if(fEventsTime) delete fEventsTime;
243}
244
245void AliGenPythiaPlus::SetInteractionRate(Float_t rate,Float_t timewindow)
246{
247// Generate pileup using user specified rate
248 fInteractionRate = rate;
249 fTimeWindow = timewindow;
250 GeneratePileup();
251}
252
253void AliGenPythiaPlus::GeneratePileup()
254{
255// Generate sub events time for pileup
256 fEventsTime = 0;
257 if(fInteractionRate == 0.) {
258 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
259 return;
260 }
261
262 Int_t npart = NumberParticles();
263 if(npart < 0) {
264 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
265 return;
266 }
267
268 if(fEventsTime) delete fEventsTime;
269 fEventsTime = new TArrayF(npart);
270 TArrayF &array = *fEventsTime;
271 for(Int_t ipart = 0; ipart < npart; ipart++)
272 array[ipart] = 0.;
273
274 Float_t eventtime = 0.;
275 while(1)
276 {
277 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
278 if(eventtime > fTimeWindow) break;
279 array.Set(array.GetSize()+1);
280 array[array.GetSize()-1] = eventtime;
281 }
282
283 eventtime = 0.;
284 while(1)
285 {
286 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
287 if(TMath::Abs(eventtime) > fTimeWindow) break;
288 array.Set(array.GetSize()+1);
289 array[array.GetSize()-1] = eventtime;
290 }
291
292 SetNumberParticles(fEventsTime->GetSize());
293}
294
295void AliGenPythiaPlus::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
296 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
297{
298// Set pycell parameters
299 fPycellEtaMax = etamax;
300 fPycellNEta = neta;
301 fPycellNPhi = nphi;
302 fPycellThreshold = thresh;
303 fPycellEtSeed = etseed;
304 fPycellMinEtJet = minet;
305 fPycellMaxRadius = r;
306}
307
308
309
310void AliGenPythiaPlus::SetEventListRange(Int_t eventFirst, Int_t eventLast)
311{
312 // Set a range of event numbers, for which a table
313 // of generated particle will be printed
314 fDebugEventFirst = eventFirst;
315 fDebugEventLast = eventLast;
316 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
317}
318
319void AliGenPythiaPlus::Init()
320{
321// Initialisation
322
323// SetMC(AliPythia::Instance());
324// fPythia=(AliPythia*) fMCEvGen;
325
326//
327 fParentWeight=1./Float_t(fNpart);
328//
329
330
331 fPythia->SetPtHardRange(fPtHardMin, fPtHardMax);
332 fPythia->SetYHardRange(fYHardMin, fYHardMax);
333
334 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
335 // Fragmentation?
336 if (fFragmentation) {
337 fPythia->SetFragmentation(1);
338 } else {
339 fPythia->SetFragmentation(0);
340 }
341
342
343// initial state radiation
344 fPythia->SetInitialAndFinalStateRadiation(fGinit, fGfinal);
345
346// pt - kick
347 fPythia->SetIntrinsicKt(fPtKick);
348
349 if (fReadFromFile) {
350 fRL = AliRunLoader::Open(fFileName, "Partons");
351 fRL->LoadKinematics();
352 fRL->LoadHeader();
353 } else {
354 fRL = 0x0;
355 }
356 //
75d4f39e 357 fPythia->ProcInit(fProcess, fEnergyCMS, fStrucFunc, fItune);
39d810c8 358 // Forward Paramters to the AliPythia object
359 fDecayer->SetForceDecay(fForceDecay);
360// Switch off Heavy Flavors on request
361 if (fHFoff) {
362 fPythia->SwitchHFOff();
363 // Switch off g->QQbar splitting in decay table
364 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
365 }
366
367 fDecayer->Init();
368
369
370// Parent and Children Selection
371 switch (fProcess)
372 {
373 case kPyOldUEQ2ordered:
374 case kPyOldUEQ2ordered2:
375 case kPyOldPopcorn:
376 break;
377 case kPyCharm:
378 case kPyCharmUnforced:
379 case kPyCharmPbPbMNR:
380 case kPyCharmpPbMNR:
381 case kPyCharmppMNR:
382 case kPyCharmppMNRwmi:
64da86aa 383 case kPyCharmPWHG:
39d810c8 384 fParentSelect[0] = 411;
385 fParentSelect[1] = 421;
386 fParentSelect[2] = 431;
387 fParentSelect[3] = 4122;
9ccc3d0e 388 fParentSelect[4] = 4232;
389 fParentSelect[5] = 4132;
390 fParentSelect[6] = 4332;
39d810c8 391 fFlavorSelect = 4;
392 break;
393 case kPyD0PbPbMNR:
394 case kPyD0pPbMNR:
395 case kPyD0ppMNR:
396 fParentSelect[0] = 421;
397 fFlavorSelect = 4;
398 break;
399 case kPyDPlusPbPbMNR:
400 case kPyDPluspPbMNR:
401 case kPyDPlusppMNR:
402 fParentSelect[0] = 411;
403 fFlavorSelect = 4;
404 break;
405 case kPyDPlusStrangePbPbMNR:
406 case kPyDPlusStrangepPbMNR:
407 case kPyDPlusStrangeppMNR:
408 fParentSelect[0] = 431;
409 fFlavorSelect = 4;
410 break;
68504d42 411 case kPyLambdacppMNR:
412 fParentSelect[0] = 4122;
413 fFlavorSelect = 4;
414 break;
39d810c8 415 case kPyBeauty:
70574ff8 416 case kPyBeautyJets:
39d810c8 417 case kPyBeautyPbPbMNR:
418 case kPyBeautypPbMNR:
419 case kPyBeautyppMNR:
420 case kPyBeautyppMNRwmi:
64da86aa 421 case kPyBeautyPWHG:
39d810c8 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;
0bd3d7c5 445 case kPyMbAtlasTuneMC09:
39d810c8 446 case kPyMbDefault:
447 case kPyMb:
04081a91 448 case kPyMbWithDirectPhoton:
39d810c8 449 case kPyMbNonDiffr:
450 case kPyMbMSEL1:
451 case kPyJets:
64da86aa 452 case kPyJetsPWHG:
39d810c8 453 case kPyDirectGamma:
454 case kPyLhwgMb:
455 break;
df607629 456 case kPyWPWHG:
39d810c8 457 case kPyW:
458 case kPyZ:
df607629 459 case kPyZgamma:
9a8774a1 460 case kPyMBRSingleDiffraction:
461 case kPyMBRDoubleDiffraction:
462 case kPyMBRCentralDiffraction:
39d810c8 463 break;
464 }
465//
466//
467// JetFinder for Trigger
468//
469// Configure detector (EMCAL like)
470//
471 fPythia->SetPycellParameters(fPycellEtaMax,fPycellNEta, fPycellNPhi,
472 fPycellThreshold, fPycellEtSeed,
473 fPycellMinEtJet, fPycellMaxRadius);
474//
475// This counts the total number of calls to Pyevnt() per run.
476 fTrialsRun = 0;
477 fQ = 0.;
478 fX1 = 0.;
479 fX2 = 0.;
480 fNev = 0 ;
481//
482//
483//
484 AliGenMC::Init();
485//
486//
487//
488 if (fSetNuclei) {
489 fDyBoost = 0;
490 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
491 }
492
493 if (fQuench) {
494 fPythia->InitQuenching(0., 0.1, 0.6e6, 0, 0.97, 30);
495 }
496
497// fPythia->SetPARJ(200, 0.0);
498
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// }
509}
510
511void AliGenPythiaPlus::Generate()
512{
513// Generate one event
514
515 fDecayer->ForceDecay();
516
b2b19810 517 Double_t polar[3] = {0,0,0};
518 Double_t origin[3] = {0,0,0};
519 Double_t p[4];
39d810c8 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;
527 fEventTime = 0.;
528
529
530
531 // Set collision vertex position
532 if (fVertexSmear == kPerEvent) Vertex();
533
534// event loop
535 while(1)
536 {
537//
538// Produce event
539//
540//
541// Switch hadronisation off
542//
543// fPythia->SwitchHadronisationOff();
544//
545// Either produce new event or read partons from file
546//
547 if (!fReadFromFile) {
548 if (!fNewMIS) {
549 fPythia->GenerateEvent();
550 } else {
551 fPythia->GenerateMIEvent();
552 }
553 fNpartons = fPythia->GetNumberOfParticles();
554 } else {
33c3c91a 555 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
556 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
39d810c8 557 fPythia->SetNumberOfParticles(0);
558 fPythia->LoadEvent(fRL->Stack(), 0 , 1);
559 fPythia->EditEventList(21);
560 }
561
562//
563// Run quenching routine
564//
565 if (fQuench == 1) {
566 fPythia->Quench();
567 } else if (fQuench == 2){
568 fPythia->Pyquen(208., 0, 0.);
569 } else if (fQuench == 3) {
570 // Quenching is via multiplicative correction of the splittings
571 }
572
573//
574// Switch hadronisation on
575//
576// fPythia->SwitchHadronisationOn();
577//
578// .. and perform hadronisation
579// printf("Calling hadronisation %d\n", fPythia->GetN());
580// fPythia->HadronizeEvent();
581 fTrials++;
8507138f 582 fPythia->GetParticles(&fParticles);
39d810c8 583 Boost();
584//
585//
586//
587 Int_t i;
588
589 fNprimaries = 0;
8507138f 590 Int_t np = fParticles.GetEntriesFast();
39d810c8 591
592 if (np == 0) continue;
593//
594
595//
596 Int_t* pParent = new Int_t[np];
597 Int_t* pSelected = new Int_t[np];
598 Int_t* trackIt = new Int_t[np];
599 for (i = 0; i < np; i++) {
600 pParent[i] = -1;
601 pSelected[i] = 0;
602 trackIt[i] = 0;
603 }
604
605 Int_t nc = 0; // Total n. of selected particles
606 Int_t nParents = 0; // Selected parents
607 Int_t nTkbles = 0; // Trackable particles
608 if (fProcess != kPyMbDefault &&
609 fProcess != kPyMb &&
04081a91 610 fProcess != kPyMbWithDirectPhoton &&
39d810c8 611 fProcess != kPyJets &&
612 fProcess != kPyDirectGamma &&
613 fProcess != kPyMbNonDiffr &&
614 fProcess != kPyMbMSEL1 &&
615 fProcess != kPyW &&
616 fProcess != kPyZ &&
df607629 617 fProcess != kPyZgamma &&
39d810c8 618 fProcess != kPyCharmppMNRwmi &&
64da86aa 619 fProcess != kPyBeautyppMNRwmi &&
df607629 620 fProcess != kPyWPWHG &&
64da86aa 621 fProcess != kPyJetsPWHG &&
622 fProcess != kPyCharmPWHG &&
623 fProcess != kPyBeautyPWHG) {
39d810c8 624
625 for (i = 0; i < np; i++) {
8507138f 626 TParticle* iparticle = (TParticle *) fParticles.At(i);
39d810c8 627 Int_t ks = iparticle->GetStatusCode();
628 kf = CheckPDGCode(iparticle->GetPdgCode());
629// No initial state partons
630 if (ks==21) continue;
631//
632// Heavy Flavor Selection
633//
634 // quark ?
635 kf = TMath::Abs(kf);
636 Int_t kfl = kf;
637 // Resonance
638
639 if (kfl > 100000) kfl %= 100000;
640 if (kfl > 10000) kfl %= 10000;
641 // meson ?
642 if (kfl > 10) kfl/=100;
643 // baryon
644 if (kfl > 10) kfl/=10;
645 Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
646 Int_t kfMo = 0;
647//
648// Establish mother daughter relation between heavy quarks and mesons
649//
650 if (kf >= fFlavorSelect && kf <= 6) {
651 Int_t idau = (fPythia->Version() == 6) ? (iparticle->GetFirstDaughter() - 1) :(iparticle->GetFirstDaughter());
652 if (idau > -1) {
8507138f 653 TParticle* daughter = (TParticle *) fParticles.At(idau);
39d810c8 654 Int_t pdgD = daughter->GetPdgCode();
655 if (pdgD == 91 || pdgD == 92) {
656 Int_t jmin = (fPythia->Version() == 6) ? (daughter->GetFirstDaughter() - 1) : (daughter->GetFirstDaughter());
657 Int_t jmax = (fPythia->Version() == 6) ? (daughter->GetLastDaughter() - 1) : (daughter->GetLastDaughter());
658
2ab330c9 659 for (Int_t jp = jmin; jp <= jmax; jp++)
660 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
39d810c8 661 } // is string or cluster
662 } // has daughter
663 } // heavy quark
664
665
666 if (ipa > -1) {
8507138f 667 TParticle * mother = (TParticle *) fParticles.At(ipa);
39d810c8 668 kfMo = TMath::Abs(mother->GetPdgCode());
669 }
670
671 // What to keep in Stack?
672 Bool_t flavorOK = kFALSE;
673 Bool_t selectOK = kFALSE;
674 if (fFeedDownOpt) {
675 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
676 } else {
677 if (kfl > fFlavorSelect) {
678 nc = -1;
679 break;
680 }
681 if (kfl == fFlavorSelect) flavorOK = kTRUE;
682 }
683 switch (fStackFillOpt) {
684 case kFlavorSelection:
685 selectOK = kTRUE;
686 break;
687 case kParentSelection:
688 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
689 break;
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 {
698 continue;
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)
706 && !KinematicSelection(iparticle, 1))
707 {
708 continue;
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;
738
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
746 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
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);
39d810c8 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();
767 p[3] = iparticle->Energy();
768
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
21391258 773 Float_t tof = fTime + kconv*iparticle->T();
39d810c8 774 Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
775 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
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);
782 pParent[i] = nt;
783 KeepTrack(nt);
784 fNprimaries++;
785 } // PushTrack loop
786 }
787 } else {
788 nc = GenerateMB();
789 } // mb ?
790
791 GetSubEventTime();
792
793 delete[] pParent;
794 delete[] pSelected;
795 delete[] trackIt;
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);
bffc134e 814 if (fInfo) fPythia->GetXandQ(fX1, fX2, fQ);
39d810c8 815 fTrialsRun += fTrials;
816 fNev++;
817 MakeHeader();
818 break;
819 }
820 }
821 } // event loop
822 SetHighWaterMark(nt);
823// Adjust weight due to kinematic selection
824// AdjustWeights();
825// get cross-section
826 fXsection = fPythia->GetXSection();
827}
828
829Int_t AliGenPythiaPlus::GenerateMB()
830{
831//
832// Min Bias selection and other global selections
833//
834 Int_t i, kf, nt, iparent;
835 Int_t nc = 0;
836 Float_t p[4];
837 Float_t polar[3] = {0,0,0};
838 Float_t origin[3] = {0,0,0};
839// converts from mm/c to s
840 const Float_t kconv = 0.001 / 2.999792458e8;
841
8507138f 842 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
39d810c8 843
844 Int_t* pParent = new Int_t[np];
845 for (i=0; i< np; i++) pParent[i] = -1;
64da86aa 846 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG ) {
8507138f 847 TParticle* jet1 = (TParticle *) fParticles.At(6);
848 TParticle* jet2 = (TParticle *) fParticles.At(7);
39d810c8 849 if (!CheckTrigger(jet1, jet2)) {
850 delete [] pParent;
851 return 0;
852 }
853 }
854
855 // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
64da86aa 856 if ((fProcess == kPyJets || fProcess == kPyJetsPWHG) && (fFragPhotonInCalo || fPi0InCalo) ) {
39d810c8 857
858 Bool_t ok = kFALSE;
859
860 Int_t pdg = 0;
861 if (fFragPhotonInCalo) pdg = 22 ; // Photon
862 else if (fPi0InCalo) pdg = 111 ; // Pi0
863
864 for (i=0; i< np; i++) {
8507138f 865 TParticle* iparticle = (TParticle *) fParticles.At(i);
39d810c8 866 if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
867 iparticle->Pt() > fFragPhotonOrPi0MinPt){
868 Int_t imother = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
8507138f 869 TParticle* pmother = (TParticle *) fParticles.At(imother);
39d810c8 870 if(pdg == 111 ||
871 (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
872 {
873 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
874 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
875 if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
876 (fCheckPHOS && IsInPHOS(phi,eta)) )
877 ok =kTRUE;
878 }
879 }
880 }
1ee02229 881 if(!ok){
882 delete [] pParent;
883 return 0;
884 }
39d810c8 885 }
886
887
888 // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
64da86aa 889 if ((fProcess == kPyJets || fProcess == kPyJetsPWHG || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
39d810c8 890
891 Bool_t okd = kFALSE;
892
893 Int_t pdg = 22;
894 Int_t iphcand = -1;
895 for (i=0; i< np; i++) {
8507138f 896 TParticle* iparticle = (TParticle *) fParticles.At(i);
39d810c8 897 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
898 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
899
900 if(iparticle->GetStatusCode() == 1
901 && iparticle->GetPdgCode() == pdg
902 && iparticle->Pt() > fPhotonMinPt
903 && eta < fPHOSEta){
904
905 // first check if the photon is in PHOS phi
906 if(IsInPHOS(phi,eta)){
907 okd = kTRUE;
908 break;
909 }
910 if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
911
912 }
913 }
914
915 if(!okd && iphcand != -1) // execute rotation in phi
916 RotatePhi(iphcand,okd);
917
92847124 918 if(!okd) {
919 delete[] pParent;
920 return 0;
921 }
39d810c8 922 }
923
924 if (fTriggerParticle) {
925 Bool_t triggered = kFALSE;
926 for (i = 0; i < np; i++) {
8507138f 927 TParticle * iparticle = (TParticle *) fParticles.At(i);
39d810c8 928 kf = CheckPDGCode(iparticle->GetPdgCode());
929 if (kf != fTriggerParticle) continue;
930 if (iparticle->Pt() == 0.) continue;
931 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
932 triggered = kTRUE;
933 break;
934 }
935 if (!triggered) {
936 delete [] pParent;
937 return 0;
938 }
939 }
940
941
942 // Check if there is a ccbar or bbbar pair with at least one of the two
943 // in fYMin < y < fYMax
944 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
9ccc3d0e 945 TParticle *partCheck;
946 TParticle *mother;
39d810c8 947 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
9ccc3d0e 948 Bool_t theChild=kFALSE;
949 Float_t y;
950 Int_t pdg,mpdg,mpdgUpperFamily;
39d810c8 951 for(i=0; i<np; i++) {
9ccc3d0e 952 partCheck = (TParticle*)fParticles.At(i);
953 pdg = partCheck->GetPdgCode();
954 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
955 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
956 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
957 (partCheck->Energy()-partCheck->Pz()+1.e-13));
d8850d4f 958 if(fUseYCutHQ && y>fYMinHQ && y<fYMaxHQ) inYcut=kTRUE;
959 if(!fUseYCutHQ && y>fYMin && y<fYMax) inYcut=kTRUE;
9ccc3d0e 960 }
961
962 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
963 Int_t mi = partCheck->GetFirstMother() - 1;
964 if(mi<0) continue;
965 mother = (TParticle*)fParticles.At(mi);
966 mpdg=TMath::Abs(mother->GetPdgCode());
75fde694 967 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
9ccc3d0e 968 if ( ParentSelected(mpdg) ||
969 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
970 if (KinematicSelection(partCheck,1)) {
971 theChild=kTRUE;
972 }
973 }
974 }
39d810c8 975 }
9ccc3d0e 976 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
39d810c8 977 delete[] pParent;
978 return 0;
979 }
9ccc3d0e 980 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
981 delete[] pParent;
982 return 0;
983 }
984
39d810c8 985 }
986
987 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
e136a92a 988 if ( (
df607629 989 fProcess == kPyWPWHG ||
e136a92a 990 fProcess == kPyW ||
39d810c8 991 fProcess == kPyZ ||
df607629 992 fProcess == kPyZgamma ||
39d810c8 993 fProcess == kPyMbDefault ||
994 fProcess == kPyMb ||
04081a91 995 fProcess == kPyMbWithDirectPhoton ||
39d810c8 996 fProcess == kPyMbNonDiffr)
997 && (fCutOnChild == 1) ) {
998 if ( !CheckKinematicsOnChild() ) {
999 delete[] pParent;
1000 return 0;
1001 }
1002 }
1003
1004
1005 for (i = 0; i < np; i++) {
1006 Int_t trackIt = 0;
8507138f 1007 TParticle * iparticle = (TParticle *) fParticles.At(i);
39d810c8 1008 kf = CheckPDGCode(iparticle->GetPdgCode());
1009 Int_t ks = iparticle->GetStatusCode();
1010 Int_t km = iparticle->GetFirstMother();
1011 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
1012 (ks != 1) ||
64da86aa 1013 ((fProcess == kPyJets || fProcess == kPyJetsPWHG) && ks == 21 && km == 0 && i>1)) {
39d810c8 1014 nc++;
1015 if (ks == 1) trackIt = 1;
1016
1017 Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
1018 iparent = (ipa > -1) ? pParent[ipa] : -1;
1019 if (ipa >= np) fPythia->EventListing();
1020
1021//
1022// store track information
1023 p[0] = iparticle->Px();
1024 p[1] = iparticle->Py();
1025 p[2] = iparticle->Pz();
1026 p[3] = iparticle->Energy();
1027
1028
1029 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1030 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1031 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1032
21391258 1033 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
39d810c8 1034
1035 PushTrack(fTrackIt*trackIt, iparent, kf,
1036 p[0], p[1], p[2], p[3],
1037 origin[0], origin[1], origin[2], tof,
1038 polar[0], polar[1], polar[2],
1039 kPPrimary, nt, 1., ks);
1040 fNprimaries++;
1041
1042
1043 //
1044 // Special Treatment to store color-flow
1045 //
3e54fb8c 1046 /*
39d810c8 1047 if (ks == 3 || ks == 13 || ks == 14) {
1048 TParticle* particle = 0;
1049 if (fStack) {
1050 particle = fStack->Particle(nt);
1051 } else {
33c3c91a 1052 particle = AliRunLoader::Instance()->Stack()->Particle(nt);
39d810c8 1053 }
3e54fb8c 1054 particle->SetFirstDaughter(fPythia->GetK(2, i));
1055 particle->SetLastDaughter(fPythia->GetK(3, i));
39d810c8 1056 }
3e54fb8c 1057 */
39d810c8 1058 KeepTrack(nt);
1059 pParent[i] = nt;
1060 SetHighWaterMark(nt);
1061
1062 } // select particle
1063 } // particle loop
1064
1065 delete[] pParent;
1066
1067 return 1;
1068}
1069
1070
1071void AliGenPythiaPlus::FinishRun()
1072{
1073// Print x-section summary
1074 fPythia->PrintStatistics();
1075
1076 if (fNev > 0.) {
1077 fQ /= fNev;
1078 fX1 /= fNev;
1079 fX2 /= fNev;
1080 }
1081
1082 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1083 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1084}
1085
1086void AliGenPythiaPlus::AdjustWeights() const
1087{
1088// Adjust the weights after generation of all events
1089//
1090 if (gAlice) {
1091 TParticle *part;
1092 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1093 for (Int_t i=0; i<ntrack; i++) {
1094 part= gAlice->GetMCApp()->Particle(i);
1095 part->SetWeight(part->GetWeight()*fKineBias);
1096 }
1097 }
1098}
1099
1100void AliGenPythiaPlus::SetNuclei(Int_t a1, Int_t a2)
1101{
1102// Treat protons as inside nuclei with mass numbers a1 and a2
1103
1104 fAProjectile = a1;
1105 fATarget = a2;
1106 fSetNuclei = kTRUE;
1107}
1108
1109
1110void AliGenPythiaPlus::MakeHeader()
1111{
1112//
1113// Make header for the simulated event
1114//
1115 if (gAlice) {
1116 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1117 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->EventListing();
1118 }
1119
1120// Builds the event header, to be called after each event
1121 if (fHeader) delete fHeader;
1122 fHeader = new AliGenPythiaEventHeader("Pythia");
d8850d4f 1123 fHeader->SetTitle(GetTitle());
1124
39d810c8 1125//
1126// Event type
1127 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->ProcessCode());
1128//
1129// Number of trials
1130 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1131//
1132// Event Vertex
1133 fHeader->SetPrimaryVertex(fVertex);
21391258 1134 fHeader->SetInteractionTime(fTime+fEventTime);
39d810c8 1135//
1136// Number of primaries
1137 fHeader->SetNProduced(fNprimaries);
1138//
1139// Jets that have triggered
1140
64da86aa 1141 if (fProcess == kPyJets || fProcess == kPyJetsPWHG)
39d810c8 1142 {
1143 Int_t ntrig, njet;
1144 Float_t jets[4][10];
1145 GetJets(njet, ntrig, jets);
1146
1147
1148 for (Int_t i = 0; i < ntrig; i++) {
1149 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1150 jets[3][i]);
1151 }
1152 }
1153//
1154// Copy relevant information from external header, if present.
1155//
1156 Float_t uqJet[4];
1157
1158 if (fRL) {
1159 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1160 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1161 {
1162 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1163
1164
1165 exHeader->TriggerJet(i, uqJet);
1166 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1167 }
1168 }
1169//
1170// Store quenching parameters
1171//
1172 if (fQuench){
3e54fb8c 1173 Double_t z[4] = {0.};
1174 Double_t xp = 0.;
1175 Double_t yp = 0.;
39d810c8 1176 if (fQuench == 1) {
1177 // Pythia::Quench()
1178 fPythia->GetQuenchingParameters(xp, yp, z);
1179 } else {
1180 // Pyquen
1181 Double_t r1 = PARIMP.rb1;
1182 Double_t r2 = PARIMP.rb2;
1183 Double_t b = PARIMP.b1;
1184 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1185 Double_t phi = PARIMP.psib1;
1186 xp = r * TMath::Cos(phi);
1187 yp = r * TMath::Sin(phi);
1188
1189 }
1190 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1191 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1192 }
1193//
1194// Store pt^hard
1195 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetPtHard());
1196//
1197// Pass header
1198//
1199 AddHeader(fHeader);
1200 fHeader = 0x0;
1201}
1202
c2fc9d31 1203Bool_t AliGenPythiaPlus::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
39d810c8 1204{
1205// Check the kinematic trigger condition
1206//
1207 Double_t eta[2];
1208 eta[0] = jet1->Eta();
1209 eta[1] = jet2->Eta();
1210 Double_t phi[2];
1211 phi[0] = jet1->Phi();
1212 phi[1] = jet2->Phi();
1213 Int_t pdg[2];
1214 pdg[0] = jet1->GetPdgCode();
1215 pdg[1] = jet2->GetPdgCode();
1216 Bool_t triggered = kFALSE;
1217
64da86aa 1218 if (fProcess == kPyJets || fProcess == kPyJetsPWHG) {
39d810c8 1219 Int_t njets = 0;
1220 Int_t ntrig = 0;
1221 Float_t jets[4][10];
1222//
1223// Use Pythia clustering on parton level to determine jet axis
1224//
1225 GetJets(njets, ntrig, jets);
1226
1227 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1228//
1229 } else {
1230 Int_t ij = 0;
1231 Int_t ig = 1;
1232 if (pdg[0] == kGamma) {
1233 ij = 1;
1234 ig = 0;
1235 }
1236 //Check eta range first...
1237 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1238 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1239 {
1240 //Eta is okay, now check phi range
1241 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1242 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1243 {
1244 triggered = kTRUE;
1245 }
1246 }
1247 }
1248 return triggered;
1249}
1250
1251
1252
1253Bool_t AliGenPythiaPlus::CheckKinematicsOnChild(){
1254//
1255//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1256//
1257 Bool_t checking = kFALSE;
1258 Int_t j, kcode, ks, km;
1259 Int_t nPartAcc = 0; //number of particles in the acceptance range
1260 Int_t numberOfAcceptedParticles = 1;
1261 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
8507138f 1262 Int_t npart = fParticles.GetEntriesFast();
39d810c8 1263
1264 for (j = 0; j<npart; j++) {
8507138f 1265 TParticle * jparticle = (TParticle *) fParticles.At(j);
39d810c8 1266 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1267 ks = jparticle->GetStatusCode();
1268 km = jparticle->GetFirstMother();
1269
1270 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1271 nPartAcc++;
1272 }
1273 if( numberOfAcceptedParticles <= nPartAcc){
1274 checking = kTRUE;
1275 break;
1276 }
1277 }
1278
1279 return checking;
1280}
1281
1282void AliGenPythiaPlus::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1283{
1284//
1285// Calls the Pythia jet finding algorithm to find jets in the current event
1286//
1287//
1288//
1289// Save jets
1290//
1291// Run Jet Finder
1292 fPythia->Pycell(njets);
1293 Int_t i;
1294 for (i = 0; i < njets; i++) {
1295 Float_t px, py, pz, e;
1296 fPythia->GetJet(i, px, py, pz, e);
1297 jets[0][i] = px;
1298 jets[1][i] = py;
1299 jets[2][i] = pz;
1300 jets[3][i] = e;
1301 }
1302}
1303
1304
1305void AliGenPythiaPlus::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1306{
1307//
1308// Calls the Pythia clustering algorithm to find jets in the current event
1309//
1310 nJets = 0;
1311 nJetsTrig = 0;
1312 if (fJetReconstruction == kCluster) {
1313//
1314// Configure cluster algorithm
1315//
1316// fPythia->SetPARU(43, 2.);
1317// fPythia->SetMSTU(41, 1);
1318//
1319// Call cluster algorithm
1320//
1321 fPythia->Pyclus(nJets);
1322//
1323// Loading jets from common block
1324//
1325 } else {
1326
1327//
1328// Run Jet Finder
1329 fPythia->Pycell(nJets);
1330 }
1331
1332 Int_t i;
1333 for (i = 0; i < nJets; i++) {
1334 Float_t px, py, pz, e;
1335 fPythia->GetJet(i, px, py, pz, e);
1336 Float_t pt = TMath::Sqrt(px * px + py * py);
1337 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
1338 Float_t theta = TMath::ATan2(pt,pz);
1339 Float_t et = e * TMath::Sin(theta);
1340 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
1341 if (
1342 eta > fEtaMinJet && eta < fEtaMaxJet &&
1343 phi > fPhiMinJet && phi < fPhiMaxJet &&
1344 et > fEtMinJet && et < fEtMaxJet
1345 )
1346 {
1347 jets[0][nJetsTrig] = px;
1348 jets[1][nJetsTrig] = py;
1349 jets[2][nJetsTrig] = pz;
1350 jets[3][nJetsTrig] = e;
1351 nJetsTrig++;
1352 } else {
1353 }
1354 }
1355}
1356
1357void AliGenPythiaPlus::GetSubEventTime()
1358{
1359 // Calculates time of the next subevent
1360 fEventTime = 0.;
1361 if (fEventsTime) {
1362 TArrayF &array = *fEventsTime;
1363 fEventTime = array[fCurSubEvent++];
1364 }
1365 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1366 return;
1367}
1368
1369
1370
1371
c2fc9d31 1372Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta) const
39d810c8 1373{
1374 // Is particle in EMCAL acceptance?
1375 // phi in degrees, etamin=-etamax
1376 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1377 eta < fEMCALEta )
1378 return kTRUE;
1379 else
1380 return kFALSE;
1381}
1382
c2fc9d31 1383Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta) const
39d810c8 1384{
1385 // Is particle in PHOS acceptance?
1386 // Acceptance slightly larger considered.
1387 // phi in degrees, etamin=-etamax
1388 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1389 eta < fPHOSEta )
1390 return kTRUE;
1391 else
1392 return kFALSE;
1393}
1394
1395void AliGenPythiaPlus::RotatePhi(Int_t iphcand, Bool_t& okdd)
1396{
1397 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1398 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1399 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1400 Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1401
1402 //calculate deltaphi
8507138f 1403 TParticle* ph = (TParticle *) fParticles.At(iphcand);
39d810c8 1404 Double_t phphi = ph->Phi();
1405 Double_t deltaphi = phiPHOS - phphi;
1406
1407
1408
1409 //loop for all particles and produce the phi rotation
8507138f 1410 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
39d810c8 1411 Double_t oldphi, newphi;
1412 Double_t newVx, newVy, R, Vz, time;
1413 Double_t newPx, newPy, pt, Pz, e;
1414 for(Int_t i=0; i< np; i++) {
8507138f 1415 TParticle* iparticle = (TParticle *) fParticles.At(i);
39d810c8 1416 oldphi = iparticle->Phi();
1417 newphi = oldphi + deltaphi;
1418 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1419 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1420
1421 R = iparticle->R();
1422 newVx = R*TMath::Cos(newphi);
1423 newVy = R*TMath::Sin(newphi);
1424 Vz = iparticle->Vz(); // don't transform
1425 time = iparticle->T(); // don't transform
1426
1427 pt = iparticle->Pt();
1428 newPx = pt*TMath::Cos(newphi);
1429 newPy = pt*TMath::Sin(newphi);
1430 Pz = iparticle->Pz(); // don't transform
1431 e = iparticle->Energy(); // don't transform
1432
1433 // apply rotation
1434 iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1435 iparticle->SetMomentum(newPx, newPy, Pz, e);
1436
1437 } //end particle loop
1438
1439 // now let's check that we put correctly the candidate photon in PHOS
1440 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1441 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1442 if(IsInPHOS(phi,eta))
1443 okdd = kTRUE;
1444}
1445
1446
1447#ifdef never
1448void AliGenPythiaPlus::Streamer(TBuffer &R__b)
1449{
1450 // Stream an object of class AliGenPythia.
1451
1452 if (R__b.IsReading()) {
1453 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1454 AliGenerator::Streamer(R__b);
1455 R__b >> (Int_t&)fProcess;
1456 R__b >> (Int_t&)fStrucFunc;
1457 R__b >> (Int_t&)fForceDecay;
1458 R__b >> fEnergyCMS;
1459 R__b >> fKineBias;
1460 R__b >> fTrials;
1461 fParentSelect.Streamer(R__b);
1462 fChildSelect.Streamer(R__b);
1463 R__b >> fXsection;
1464// (AliPythia::Instance())->Streamer(R__b);
1465 R__b >> fPtHardMin;
1466 R__b >> fPtHardMax;
1467// if (fDecayer) fDecayer->Streamer(R__b);
1468 } else {
1469 R__b.WriteVersion(AliGenPythiaPlus::IsA());
1470 AliGenerator::Streamer(R__b);
1471 R__b << (Int_t)fProcess;
1472 R__b << (Int_t)fStrucFunc;
1473 R__b << (Int_t)fForceDecay;
1474 R__b << fEnergyCMS;
1475 R__b << fKineBias;
1476 R__b << fTrials;
1477 fParentSelect.Streamer(R__b);
1478 fChildSelect.Streamer(R__b);
1479 R__b << fXsection;
1480// R__b << fPythia;
1481 R__b << fPtHardMin;
1482 R__b << fPtHardMax;
1483 // fDecayer->Streamer(R__b);
1484 }
1485}
1486#endif
1487
1488
1489