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