]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PYTHIA6/AliGenPythia.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Test script updated
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
... / ...
CommitLineData
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 <TObjArray.h>
33#include <TSystem.h>
34#include <TTree.h>
35#include "AliConst.h"
36#include "AliDecayerPythia.h"
37#include "AliGenPythia.h"
38#include "AliFastGlauber.h"
39#include "AliHeader.h"
40#include "AliGenPythiaEventHeader.h"
41#include "AliPythia.h"
42#include "AliPythiaRndm.h"
43#include "AliRun.h"
44#include "AliStack.h"
45#include "AliRunLoader.h"
46#include "AliMC.h"
47#include "AliLog.h"
48#include "PyquenCommon.h"
49#include "AliLog.h"
50
51ClassImp(AliGenPythia)
52
53
54AliGenPythia::AliGenPythia():
55 AliGenMC(),
56 fProcess(kPyCharm),
57 fItune(-1),
58 fStrucFunc(kCTEQ5L),
59 fKineBias(0.),
60 fTrials(0),
61 fTrialsRun(0),
62 fQ(0.),
63 fX1(0.),
64 fX2(0.),
65 fEventTime(0.),
66 fInteractionRate(0.),
67 fTimeWindow(0.),
68 fCurSubEvent(0),
69 fEventsTime(0),
70 fNev(0),
71 fFlavorSelect(0),
72 fXsection(0.),
73 fPythia(0),
74 fPtHardMin(0.),
75 fPtHardMax(1.e4),
76 fYHardMin(-1.e10),
77 fYHardMax(1.e10),
78 fGinit(1),
79 fGfinal(1),
80 fHadronisation(1),
81 fPatchOmegaDalitz(0),
82 fNpartons(0),
83 fReadFromFile(0),
84 fQuench(0),
85 fQhat(0.),
86 fLength(0.),
87 fpyquenT(1.),
88 fpyquenTau(0.1),
89 fpyquenNf(0),
90 fpyquenEloss(0),
91 fpyquenAngle(0),
92 fImpact(0.),
93 fPtKick(1.),
94 fFullEvent(kTRUE),
95 fDecayer(new AliDecayerPythia()),
96 fDebugEventFirst(-1),
97 fDebugEventLast(-1),
98 fEtMinJet(0.),
99 fEtMaxJet(1.e4),
100 fEtaMinJet(-20.),
101 fEtaMaxJet(20.),
102 fPhiMinJet(0.),
103 fPhiMaxJet(2.* TMath::Pi()),
104 fJetReconstruction(kCell),
105 fEtaMinGamma(-20.),
106 fEtaMaxGamma(20.),
107 fPhiMinGamma(0.),
108 fPhiMaxGamma(2. * TMath::Pi()),
109 fPycellEtaMax(2.),
110 fPycellNEta(274),
111 fPycellNPhi(432),
112 fPycellThreshold(0.),
113 fPycellEtSeed(4.),
114 fPycellMinEtJet(10.),
115 fPycellMaxRadius(1.),
116 fStackFillOpt(kFlavorSelection),
117 fFeedDownOpt(kTRUE),
118 fFragmentation(kTRUE),
119 fSetNuclei(kFALSE),
120 fNewMIS(kFALSE),
121 fHFoff(kFALSE),
122 fNucPdf(0),
123 fTriggerParticle(0),
124 fTriggerEta(0.9),
125 fTriggerMinPt(-1),
126 fTriggerMaxPt(1000),
127 fTriggerMultiplicity(0),
128 fTriggerMultiplicityEta(0),
129 fTriggerMultiplicityPtMin(0),
130 fCountMode(kCountAll),
131 fHeader(0),
132 fRL(0),
133 fkFileName(0),
134 fFragPhotonInCalo(kFALSE),
135 fHadronInCalo(kFALSE) ,
136 fPi0InCalo(kFALSE) ,
137 fEtaInCalo(kFALSE) ,
138 fPhotonInCalo(kFALSE), // not in use
139 fDecayPhotonInCalo(kFALSE),
140 fForceNeutralMeson2PhotonDecay(kFALSE),
141 fEleInCalo(kFALSE),
142 fEleInEMCAL(kFALSE), // not in use
143 fCheckBarrel(kFALSE),
144 fCheckEMCAL(kFALSE),
145 fCheckPHOS(kFALSE),
146 fCheckPHOSeta(kFALSE),
147 fPHOSRotateCandidate(-1),
148 fTriggerParticleMinPt(0),
149 fPhotonMinPt(0), // not in use
150 fElectronMinPt(0), // not in use
151 fPHOSMinPhi(219.),
152 fPHOSMaxPhi(321.),
153 fPHOSEta(0.13),
154 fEMCALMinPhi(79.),
155 fEMCALMaxPhi(191.),
156 fEMCALEta(0.71),
157 fkTuneForDiff(0),
158 fProcDiff(0)
159{
160// Default Constructor
161 fEnergyCMS = 5500.;
162 if (!AliPythiaRndm::GetPythiaRandom())
163 AliPythiaRndm::SetPythiaRandom(GetRandom());
164}
165
166AliGenPythia::AliGenPythia(Int_t npart)
167 :AliGenMC(npart),
168 fProcess(kPyCharm),
169 fItune(-1),
170 fStrucFunc(kCTEQ5L),
171 fKineBias(0.),
172 fTrials(0),
173 fTrialsRun(0),
174 fQ(0.),
175 fX1(0.),
176 fX2(0.),
177 fEventTime(0.),
178 fInteractionRate(0.),
179 fTimeWindow(0.),
180 fCurSubEvent(0),
181 fEventsTime(0),
182 fNev(0),
183 fFlavorSelect(0),
184 fXsection(0.),
185 fPythia(0),
186 fPtHardMin(0.),
187 fPtHardMax(1.e4),
188 fYHardMin(-1.e10),
189 fYHardMax(1.e10),
190 fGinit(kTRUE),
191 fGfinal(kTRUE),
192 fHadronisation(kTRUE),
193 fPatchOmegaDalitz(0),
194 fNpartons(0),
195 fReadFromFile(kFALSE),
196 fQuench(kFALSE),
197 fQhat(0.),
198 fLength(0.),
199 fpyquenT(1.),
200 fpyquenTau(0.1),
201 fpyquenNf(0),
202 fpyquenEloss(0),
203 fpyquenAngle(0),
204 fImpact(0.),
205 fPtKick(1.),
206 fFullEvent(kTRUE),
207 fDecayer(new AliDecayerPythia()),
208 fDebugEventFirst(-1),
209 fDebugEventLast(-1),
210 fEtMinJet(0.),
211 fEtMaxJet(1.e4),
212 fEtaMinJet(-20.),
213 fEtaMaxJet(20.),
214 fPhiMinJet(0.),
215 fPhiMaxJet(2.* TMath::Pi()),
216 fJetReconstruction(kCell),
217 fEtaMinGamma(-20.),
218 fEtaMaxGamma(20.),
219 fPhiMinGamma(0.),
220 fPhiMaxGamma(2. * TMath::Pi()),
221 fPycellEtaMax(2.),
222 fPycellNEta(274),
223 fPycellNPhi(432),
224 fPycellThreshold(0.),
225 fPycellEtSeed(4.),
226 fPycellMinEtJet(10.),
227 fPycellMaxRadius(1.),
228 fStackFillOpt(kFlavorSelection),
229 fFeedDownOpt(kTRUE),
230 fFragmentation(kTRUE),
231 fSetNuclei(kFALSE),
232 fNewMIS(kFALSE),
233 fHFoff(kFALSE),
234 fNucPdf(0),
235 fTriggerParticle(0),
236 fTriggerEta(0.9),
237 fTriggerMinPt(-1),
238 fTriggerMaxPt(1000),
239 fTriggerMultiplicity(0),
240 fTriggerMultiplicityEta(0),
241 fTriggerMultiplicityPtMin(0),
242 fCountMode(kCountAll),
243 fHeader(0),
244 fRL(0),
245 fkFileName(0),
246 fFragPhotonInCalo(kFALSE),
247 fHadronInCalo(kFALSE) ,
248 fPi0InCalo(kFALSE) ,
249 fEtaInCalo(kFALSE) ,
250 fPhotonInCalo(kFALSE), // not in use
251 fDecayPhotonInCalo(kFALSE),
252 fForceNeutralMeson2PhotonDecay(kFALSE),
253 fEleInCalo(kFALSE),
254 fEleInEMCAL(kFALSE), // not in use
255 fCheckBarrel(kFALSE),
256 fCheckEMCAL(kFALSE),
257 fCheckPHOS(kFALSE),
258 fCheckPHOSeta(kFALSE),
259 fPHOSRotateCandidate(-1),
260 fTriggerParticleMinPt(0),
261 fPhotonMinPt(0), // not in use
262 fElectronMinPt(0), // not in use
263 fPHOSMinPhi(219.),
264 fPHOSMaxPhi(321.),
265 fPHOSEta(0.13),
266 fEMCALMinPhi(79.),
267 fEMCALMaxPhi(191.),
268 fEMCALEta(0.71),
269 fkTuneForDiff(0),
270 fProcDiff(0)
271{
272// default charm production at 5. 5 TeV
273// semimuonic decay
274// structure function GRVHO
275//
276 fEnergyCMS = 5500.;
277 fName = "Pythia";
278 fTitle= "Particle Generator using PYTHIA";
279 SetForceDecay();
280 // Set random number generator
281 if (!AliPythiaRndm::GetPythiaRandom())
282 AliPythiaRndm::SetPythiaRandom(GetRandom());
283 }
284
285AliGenPythia::~AliGenPythia()
286{
287// Destructor
288 if(fEventsTime) delete fEventsTime;
289}
290
291void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
292{
293// Generate pileup using user specified rate
294 fInteractionRate = rate;
295 fTimeWindow = timewindow;
296 GeneratePileup();
297}
298
299void AliGenPythia::GeneratePileup()
300{
301// Generate sub events time for pileup
302 fEventsTime = 0;
303 if(fInteractionRate == 0.) {
304 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
305 return;
306 }
307
308 Int_t npart = NumberParticles();
309 if(npart < 0) {
310 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
311 return;
312 }
313
314 if(fEventsTime) delete fEventsTime;
315 fEventsTime = new TArrayF(npart);
316 TArrayF &array = *fEventsTime;
317 for(Int_t ipart = 0; ipart < npart; ipart++)
318 array[ipart] = 0.;
319
320 Float_t eventtime = 0.;
321 while(1)
322 {
323 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
324 if(eventtime > fTimeWindow) break;
325 array.Set(array.GetSize()+1);
326 array[array.GetSize()-1] = eventtime;
327 }
328
329 eventtime = 0.;
330 while(1)
331 {
332 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
333 if(TMath::Abs(eventtime) > fTimeWindow) break;
334 array.Set(array.GetSize()+1);
335 array[array.GetSize()-1] = eventtime;
336 }
337
338 SetNumberParticles(fEventsTime->GetSize());
339}
340
341void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
342 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
343{
344// Set pycell parameters
345 fPycellEtaMax = etamax;
346 fPycellNEta = neta;
347 fPycellNPhi = nphi;
348 fPycellThreshold = thresh;
349 fPycellEtSeed = etseed;
350 fPycellMinEtJet = minet;
351 fPycellMaxRadius = r;
352}
353
354
355
356void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
357{
358 // Set a range of event numbers, for which a table
359 // of generated particle will be printed
360 fDebugEventFirst = eventFirst;
361 fDebugEventLast = eventLast;
362 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
363}
364
365void AliGenPythia::Init()
366{
367// Initialisation
368
369 SetMC(AliPythia::Instance());
370 fPythia=(AliPythia*) fMCEvGen;
371
372//
373 fParentWeight=1./Float_t(fNpart);
374//
375
376
377 fPythia->SetCKIN(3,fPtHardMin);
378 fPythia->SetCKIN(4,fPtHardMax);
379 fPythia->SetCKIN(7,fYHardMin);
380 fPythia->SetCKIN(8,fYHardMax);
381
382 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
383 // Fragmentation?
384 if (fFragmentation) {
385 fPythia->SetMSTP(111,1);
386 } else {
387 fPythia->SetMSTP(111,0);
388 }
389
390
391// initial state radiation
392 fPythia->SetMSTP(61,fGinit);
393// final state radiation
394 fPythia->SetMSTP(71,fGfinal);
395// pt - kick
396 if (fPtKick > 0.) {
397 fPythia->SetMSTP(91,1);
398 fPythia->SetPARP(91,fPtKick);
399 fPythia->SetPARP(93, 4. * fPtKick);
400 } else {
401 fPythia->SetMSTP(91,0);
402 }
403
404
405 if (fReadFromFile) {
406 fRL = AliRunLoader::Open(fkFileName, "Partons");
407 fRL->LoadKinematics();
408 fRL->LoadHeader();
409 } else {
410 fRL = 0x0;
411 }
412 //
413 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
414 // Forward Paramters to the AliPythia object
415 fDecayer->SetForceDecay(fForceDecay);
416// Switch off Heavy Flavors on request
417 if (fHFoff) {
418 // Maximum number of quark flavours used in pdf
419 fPythia->SetMSTP(58, 3);
420 // Maximum number of flavors that can be used in showers
421 fPythia->SetMSTJ(45, 3);
422 // Switch off g->QQbar splitting in decay table
423 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
424 }
425
426 fDecayer->Init();
427
428
429// Parent and Children Selection
430 switch (fProcess)
431 {
432 case kPyOldUEQ2ordered:
433 case kPyOldUEQ2ordered2:
434 case kPyOldPopcorn:
435 break;
436 case kPyCharm:
437 case kPyCharmUnforced:
438 case kPyCharmPbPbMNR:
439 case kPyCharmpPbMNR:
440 case kPyCharmppMNR:
441 case kPyCharmppMNRwmi:
442 fParentSelect[0] = 411;
443 fParentSelect[1] = 421;
444 fParentSelect[2] = 431;
445 fParentSelect[3] = 4122;
446 fParentSelect[4] = 4232;
447 fParentSelect[5] = 4132;
448 fParentSelect[6] = 4332;
449 fFlavorSelect = 4;
450 break;
451 case kPyD0PbPbMNR:
452 case kPyD0pPbMNR:
453 case kPyD0ppMNR:
454 fParentSelect[0] = 421;
455 fFlavorSelect = 4;
456 break;
457 case kPyDPlusPbPbMNR:
458 case kPyDPluspPbMNR:
459 case kPyDPlusppMNR:
460 fParentSelect[0] = 411;
461 fFlavorSelect = 4;
462 break;
463 case kPyDPlusStrangePbPbMNR:
464 case kPyDPlusStrangepPbMNR:
465 case kPyDPlusStrangeppMNR:
466 fParentSelect[0] = 431;
467 fFlavorSelect = 4;
468 break;
469 case kPyLambdacppMNR:
470 fParentSelect[0] = 4122;
471 fFlavorSelect = 4;
472 break;
473 case kPyBeauty:
474 case kPyBeautyJets:
475 case kPyBeautyPbPbMNR:
476 case kPyBeautypPbMNR:
477 case kPyBeautyppMNR:
478 case kPyBeautyppMNRwmi:
479 fParentSelect[0]= 511;
480 fParentSelect[1]= 521;
481 fParentSelect[2]= 531;
482 fParentSelect[3]= 5122;
483 fParentSelect[4]= 5132;
484 fParentSelect[5]= 5232;
485 fParentSelect[6]= 5332;
486 fFlavorSelect = 5;
487 break;
488 case kPyBeautyUnforced:
489 fParentSelect[0] = 511;
490 fParentSelect[1] = 521;
491 fParentSelect[2] = 531;
492 fParentSelect[3] = 5122;
493 fParentSelect[4] = 5132;
494 fParentSelect[5] = 5232;
495 fParentSelect[6] = 5332;
496 fFlavorSelect = 5;
497 break;
498 case kPyJpsiChi:
499 case kPyJpsi:
500 fParentSelect[0] = 443;
501 break;
502 case kPyMbDefault:
503 case kPyMbAtlasTuneMC09:
504 case kPyMb:
505 case kPyMbWithDirectPhoton:
506 case kPyMbNonDiffr:
507 case kPyMbMSEL1:
508 case kPyJets:
509 case kPyDirectGamma:
510 case kPyLhwgMb:
511 break;
512 case kPyW:
513 case kPyZ:
514 break;
515 }
516//
517//
518// JetFinder for Trigger
519//
520// Configure detector (EMCAL like)
521//
522 fPythia->SetPARU(51, fPycellEtaMax);
523 fPythia->SetMSTU(51, fPycellNEta);
524 fPythia->SetMSTU(52, fPycellNPhi);
525//
526// Configure Jet Finder
527//
528 fPythia->SetPARU(58, fPycellThreshold);
529 fPythia->SetPARU(52, fPycellEtSeed);
530 fPythia->SetPARU(53, fPycellMinEtJet);
531 fPythia->SetPARU(54, fPycellMaxRadius);
532 fPythia->SetMSTU(54, 2);
533//
534// This counts the total number of calls to Pyevnt() per run.
535 fTrialsRun = 0;
536 fQ = 0.;
537 fX1 = 0.;
538 fX2 = 0.;
539 fNev = 0 ;
540//
541//
542//
543 AliGenMC::Init();
544//
545//
546//
547 if (fSetNuclei) {
548 fDyBoost = 0;
549 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
550 }
551
552 fPythia->SetPARJ(200, 0.0);
553 fPythia->SetPARJ(199, 0.0);
554 fPythia->SetPARJ(198, 0.0);
555 fPythia->SetPARJ(197, 0.0);
556
557 if (fQuench == 1) {
558 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
559 }
560
561 if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
562
563 if (fQuench == 3) {
564 // Nestor's change of the splittings
565 fPythia->SetPARJ(200, 0.8);
566 fPythia->SetMSTJ(41, 1); // QCD radiation only
567 fPythia->SetMSTJ(42, 2); // angular ordering
568 fPythia->SetMSTJ(44, 2); // option to run alpha_s
569 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
570 fPythia->SetMSTJ(50, 0); // No coherence in first branching
571 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
572 } else if (fQuench == 4) {
573 // Armesto-Cunqueiro-Salgado change of the splittings.
574 AliFastGlauber* glauber = AliFastGlauber::Instance();
575 glauber->Init(2);
576 //read and store transverse almonds corresponding to differnt
577 //impact parameters.
578 glauber->SetCentralityClass(0.,0.1);
579 fPythia->SetPARJ(200, 1.);
580 fPythia->SetPARJ(198, fQhat);
581 fPythia->SetPARJ(199, fLength);
582 fPythia->SetMSTJ(42, 2); // angular ordering
583 fPythia->SetMSTJ(44, 2); // option to run alpha_s
584 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
585 }
586}
587
588void AliGenPythia::Generate()
589{
590// Generate one event
591 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
592 fDecayer->ForceDecay();
593
594 Double_t polar[3] = {0,0,0};
595 Double_t origin[3] = {0,0,0};
596 Double_t p[4];
597// converts from mm/c to s
598 const Float_t kconv=0.001/2.999792458e8;
599//
600 Int_t nt=0;
601 Int_t jev=0;
602 Int_t j, kf;
603 fTrials=0;
604 fEventTime = 0.;
605
606
607
608 // Set collision vertex position
609 if (fVertexSmear == kPerEvent) Vertex();
610
611// event loop
612 while(1)
613 {
614//
615// Produce event
616//
617//
618// Switch hadronisation off
619//
620 fPythia->SetMSTJ(1, 0);
621
622 if (fQuench ==4){
623 Double_t bimp;
624 // Quenching comes through medium-modified splitting functions.
625 AliFastGlauber::Instance()->GetRandomBHard(bimp);
626 fPythia->SetPARJ(197, bimp);
627 fImpact = bimp;
628 fPythia->Qpygin0();
629 }
630//
631// Either produce new event or read partons from file
632//
633 if (!fReadFromFile) {
634 if (!fNewMIS) {
635 fPythia->Pyevnt();
636 } else {
637 fPythia->Pyevnw();
638 }
639 fNpartons = fPythia->GetN();
640 } else {
641 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
642 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
643 fPythia->SetN(0);
644 LoadEvent(fRL->Stack(), 0 , 1);
645 fPythia->Pyedit(21);
646 }
647
648//
649// Run quenching routine
650//
651 if (fQuench == 1) {
652 fPythia->Quench();
653 } else if (fQuench == 2){
654 fPythia->Pyquen(208., 0, 0.);
655 } else if (fQuench == 3) {
656 // Quenching is via multiplicative correction of the splittings
657 }
658
659//
660// Switch hadronisation on
661//
662 if (fHadronisation) {
663 fPythia->SetMSTJ(1, 1);
664//
665// .. and perform hadronisation
666// printf("Calling hadronisation %d\n", fPythia->GetN());
667
668 if (fPatchOmegaDalitz) {
669 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
670 fPythia->Pyexec();
671 fPythia->DalitzDecays();
672 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
673 }
674 fPythia->Pyexec();
675 }
676 fTrials++;
677 fPythia->ImportParticles(&fParticles,"All");
678
679 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
680 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
681
682//
683//
684//
685 Int_t i;
686
687 fNprimaries = 0;
688 Int_t np = fParticles.GetEntriesFast();
689
690 if (np == 0) continue;
691//
692
693//
694 Int_t* pParent = new Int_t[np];
695 Int_t* pSelected = new Int_t[np];
696 Int_t* trackIt = new Int_t[np];
697 for (i = 0; i < np; i++) {
698 pParent[i] = -1;
699 pSelected[i] = 0;
700 trackIt[i] = 0;
701 }
702
703 Int_t nc = 0; // Total n. of selected particles
704 Int_t nParents = 0; // Selected parents
705 Int_t nTkbles = 0; // Trackable particles
706 if (fProcess != kPyMbDefault &&
707 fProcess != kPyMb &&
708 fProcess != kPyMbAtlasTuneMC09 &&
709 fProcess != kPyMbWithDirectPhoton &&
710 fProcess != kPyJets &&
711 fProcess != kPyDirectGamma &&
712 fProcess != kPyMbNonDiffr &&
713 fProcess != kPyMbMSEL1 &&
714 fProcess != kPyW &&
715 fProcess != kPyZ &&
716 fProcess != kPyCharmppMNRwmi &&
717 fProcess != kPyBeautyppMNRwmi &&
718 fProcess != kPyBeautyJets) {
719
720 for (i = 0; i < np; i++) {
721 TParticle* iparticle = (TParticle *) fParticles.At(i);
722 Int_t ks = iparticle->GetStatusCode();
723 kf = CheckPDGCode(iparticle->GetPdgCode());
724// No initial state partons
725 if (ks==21) continue;
726//
727// Heavy Flavor Selection
728//
729 // quark ?
730 kf = TMath::Abs(kf);
731 Int_t kfl = kf;
732 // Resonance
733
734 if (kfl > 100000) kfl %= 100000;
735 if (kfl > 10000) kfl %= 10000;
736 // meson ?
737 if (kfl > 10) kfl/=100;
738 // baryon
739 if (kfl > 10) kfl/=10;
740 Int_t ipa = iparticle->GetFirstMother()-1;
741 Int_t kfMo = 0;
742//
743// Establish mother daughter relation between heavy quarks and mesons
744//
745 if (kf >= fFlavorSelect && kf <= 6) {
746 Int_t idau = iparticle->GetFirstDaughter() - 1;
747 if (idau > -1) {
748 TParticle* daughter = (TParticle *) fParticles.At(idau);
749 Int_t pdgD = daughter->GetPdgCode();
750 if (pdgD == 91 || pdgD == 92) {
751 Int_t jmin = daughter->GetFirstDaughter() - 1;
752 Int_t jmax = daughter->GetLastDaughter() - 1;
753 for (Int_t jp = jmin; jp <= jmax; jp++)
754 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
755 } // is string or cluster
756 } // has daughter
757 } // heavy quark
758
759
760 if (ipa > -1) {
761 TParticle * mother = (TParticle *) fParticles.At(ipa);
762 kfMo = TMath::Abs(mother->GetPdgCode());
763 }
764
765 // What to keep in Stack?
766 Bool_t flavorOK = kFALSE;
767 Bool_t selectOK = kFALSE;
768 if (fFeedDownOpt) {
769 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
770 } else {
771 if (kfl > fFlavorSelect) {
772 nc = -1;
773 break;
774 }
775 if (kfl == fFlavorSelect) flavorOK = kTRUE;
776 }
777 switch (fStackFillOpt) {
778 case kHeavyFlavor:
779 case kFlavorSelection:
780 selectOK = kTRUE;
781 break;
782 case kParentSelection:
783 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
784 break;
785 }
786 if (flavorOK && selectOK) {
787//
788// Heavy flavor hadron or quark
789//
790// Kinematic seletion on final state heavy flavor mesons
791 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
792 {
793 continue;
794 }
795 pSelected[i] = 1;
796 if (ParentSelected(kf)) ++nParents; // Update parent count
797// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
798 } else {
799// Kinematic seletion on decay products
800 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
801 && !KinematicSelection(iparticle, 1))
802 {
803 continue;
804 }
805//
806// Decay products
807// Select if mother was selected and is not tracked
808
809 if (pSelected[ipa] &&
810 !trackIt[ipa] && // mother will be tracked ?
811 kfMo != 5 && // mother is b-quark, don't store fragments
812 kfMo != 4 && // mother is c-quark, don't store fragments
813 kf != 92) // don't store string
814 {
815//
816// Semi-stable or de-selected: diselect decay products:
817//
818//
819 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
820 {
821 Int_t ipF = iparticle->GetFirstDaughter();
822 Int_t ipL = iparticle->GetLastDaughter();
823 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
824 }
825// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
826 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
827 }
828 }
829 if (pSelected[i] == -1) pSelected[i] = 0;
830 if (!pSelected[i]) continue;
831 // Count quarks only if you did not include fragmentation
832 if (fFragmentation && kf <= 10) continue;
833
834 nc++;
835// Decision on tracking
836 trackIt[i] = 0;
837//
838// Track final state particle
839 if (ks == 1) trackIt[i] = 1;
840// Track semi-stable particles
841 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
842// Track particles selected by process if undecayed.
843 if (fForceDecay == kNoDecay) {
844 if (ParentSelected(kf)) trackIt[i] = 1;
845 } else {
846 if (ParentSelected(kf)) trackIt[i] = 0;
847 }
848 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
849//
850//
851
852 } // particle selection loop
853 if (nc > 0) {
854 for (i = 0; i<np; i++) {
855 if (!pSelected[i]) continue;
856 TParticle * iparticle = (TParticle *) fParticles.At(i);
857 kf = CheckPDGCode(iparticle->GetPdgCode());
858 Int_t ks = iparticle->GetStatusCode();
859 p[0] = iparticle->Px();
860 p[1] = iparticle->Py();
861 p[2] = iparticle->Pz();
862 p[3] = iparticle->Energy();
863
864 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
865 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
866 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
867
868 Float_t tof = fTime + kconv*iparticle->T();
869 Int_t ipa = iparticle->GetFirstMother()-1;
870 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
871
872 PushTrack(fTrackIt*trackIt[i], iparent, kf,
873 p[0], p[1], p[2], p[3],
874 origin[0], origin[1], origin[2], tof,
875 polar[0], polar[1], polar[2],
876 kPPrimary, nt, 1., ks);
877 pParent[i] = nt;
878 KeepTrack(nt);
879 fNprimaries++;
880 } // PushTrack loop
881 }
882 } else {
883 nc = GenerateMB();
884 } // mb ?
885
886 GetSubEventTime();
887
888 delete[] pParent;
889 delete[] pSelected;
890 delete[] trackIt;
891
892 if (nc > 0) {
893 switch (fCountMode) {
894 case kCountAll:
895 // printf(" Count all \n");
896 jev += nc;
897 break;
898 case kCountParents:
899 // printf(" Count parents \n");
900 jev += nParents;
901 break;
902 case kCountTrackables:
903 // printf(" Count trackable \n");
904 jev += nTkbles;
905 break;
906 }
907 if (jev >= fNpart || fNpart == -1) {
908 fKineBias=Float_t(fNpart)/Float_t(fTrials);
909
910 fQ += fPythia->GetVINT(51);
911 fX1 += fPythia->GetVINT(41);
912 fX2 += fPythia->GetVINT(42);
913 fTrialsRun += fTrials;
914 fNev++;
915 MakeHeader();
916 break;
917 }
918 }
919 } // event loop
920 SetHighWaterMark(nt);
921// adjust weight due to kinematic selection
922// AdjustWeights();
923// get cross-section
924 fXsection=fPythia->GetPARI(1);
925}
926
927Int_t AliGenPythia::GenerateMB()
928{
929//
930// Min Bias selection and other global selections
931//
932 Int_t i, kf, nt, iparent;
933 Int_t nc = 0;
934 Double_t p[4];
935 Double_t polar[3] = {0,0,0};
936 Double_t origin[3] = {0,0,0};
937// converts from mm/c to s
938 const Float_t kconv=0.001/2.999792458e8;
939
940
941 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
942
943
944
945 Int_t* pParent = new Int_t[np];
946 for (i=0; i< np; i++) pParent[i] = -1;
947 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
948 TParticle* jet1 = (TParticle *) fParticles.At(6);
949 TParticle* jet2 = (TParticle *) fParticles.At(7);
950 if (!CheckTrigger(jet1, jet2)) {
951 delete [] pParent;
952 return 0;
953 }
954 }
955
956
957 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
958 // implemented primaryly for kPyJets, but extended to any kind of process.
959 if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
960 (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
961 Bool_t ok = TriggerOnSelectedParticles(np);
962
963 if(!ok) {
964 delete[] pParent;
965 return 0;
966 }
967 }
968
969 // Check for diffraction
970 if(fkTuneForDiff) {
971 if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
972 if(!CheckDiffraction()) {
973 delete [] pParent;
974 return 0;
975 }
976 }
977 }
978
979 // Check for minimum multiplicity
980 if (fTriggerMultiplicity > 0) {
981 Int_t multiplicity = 0;
982 for (i = 0; i < np; i++) {
983 TParticle * iparticle = (TParticle *) fParticles.At(i);
984
985 Int_t statusCode = iparticle->GetStatusCode();
986
987 // Initial state particle
988 if (statusCode != 1)
989 continue;
990 // eta cut
991 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
992 continue;
993 // pt cut
994 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
995 continue;
996
997 TParticlePDG* pdgPart = iparticle->GetPDG();
998 if (pdgPart && pdgPart->Charge() == 0)
999 continue;
1000
1001 ++multiplicity;
1002 }
1003
1004 if (multiplicity < fTriggerMultiplicity) {
1005 delete [] pParent;
1006 return 0;
1007 }
1008 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
1009 }
1010
1011
1012 if (fTriggerParticle) {
1013 Bool_t triggered = kFALSE;
1014 for (i = 0; i < np; i++) {
1015 TParticle * iparticle = (TParticle *) fParticles.At(i);
1016 kf = CheckPDGCode(iparticle->GetPdgCode());
1017 if (kf != fTriggerParticle) continue;
1018 if (iparticle->Pt() == 0.) continue;
1019 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1020 if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
1021 triggered = kTRUE;
1022 break;
1023 }
1024 if (!triggered) {
1025 delete [] pParent;
1026 return 0;
1027 }
1028 }
1029
1030
1031 // Check if there is a ccbar or bbbar pair with at least one of the two
1032 // in fYMin < y < fYMax
1033
1034 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1035 TParticle *partCheck;
1036 TParticle *mother;
1037 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1038 Bool_t theChild=kFALSE;
1039 Bool_t triggered=kFALSE;
1040 Float_t y;
1041 Int_t pdg,mpdg,mpdgUpperFamily;
1042 for(i=0; i<np; i++) {
1043 partCheck = (TParticle*)fParticles.At(i);
1044 pdg = partCheck->GetPdgCode();
1045 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1046 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1047 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1048 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1049 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1050 }
1051 if(fTriggerParticle) {
1052 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1053 }
1054 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1055 Int_t mi = partCheck->GetFirstMother() - 1;
1056 if(mi<0) continue;
1057 mother = (TParticle*)fParticles.At(mi);
1058 mpdg=TMath::Abs(mother->GetPdgCode());
1059 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1060 if ( ParentSelected(mpdg) ||
1061 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1062 if (KinematicSelection(partCheck,1)) {
1063 theChild=kTRUE;
1064 }
1065 }
1066 }
1067 }
1068 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1069 delete[] pParent;
1070 return 0;
1071 }
1072 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1073 delete[] pParent;
1074 return 0;
1075 }
1076 if(fTriggerParticle && !triggered) { // particle requested is not produced
1077 delete[] pParent;
1078 return 0;
1079 }
1080
1081 }
1082
1083 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1084 if ( (fProcess == kPyW ||
1085 fProcess == kPyZ ||
1086 fProcess == kPyMbDefault ||
1087 fProcess == kPyMb ||
1088 fProcess == kPyMbAtlasTuneMC09 ||
1089 fProcess == kPyMbWithDirectPhoton ||
1090 fProcess == kPyMbNonDiffr)
1091 && (fCutOnChild == 1) ) {
1092 if ( !CheckKinematicsOnChild() ) {
1093 delete[] pParent;
1094 return 0;
1095 }
1096 }
1097
1098
1099 for (i = 0; i < np; i++) {
1100 Int_t trackIt = 0;
1101 TParticle * iparticle = (TParticle *) fParticles.At(i);
1102 kf = CheckPDGCode(iparticle->GetPdgCode());
1103 Int_t ks = iparticle->GetStatusCode();
1104 Int_t km = iparticle->GetFirstMother();
1105 if (
1106 (((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
1107 ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
1108 )
1109 {
1110 nc++;
1111 if (ks == 1) trackIt = 1;
1112 Int_t ipa = iparticle->GetFirstMother()-1;
1113
1114 iparent = (ipa > -1) ? pParent[ipa] : -1;
1115
1116//
1117// store track information
1118 p[0] = iparticle->Px();
1119 p[1] = iparticle->Py();
1120 p[2] = iparticle->Pz();
1121 p[3] = iparticle->Energy();
1122
1123
1124 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1125 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1126 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1127
1128 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
1129
1130 PushTrack(fTrackIt*trackIt, iparent, kf,
1131 p[0], p[1], p[2], p[3],
1132 origin[0], origin[1], origin[2], tof,
1133 polar[0], polar[1], polar[2],
1134 kPPrimary, nt, 1., ks);
1135 fNprimaries++;
1136 KeepTrack(nt);
1137 pParent[i] = nt;
1138 SetHighWaterMark(nt);
1139
1140 } // select particle
1141 } // particle loop
1142
1143 delete[] pParent;
1144
1145 return 1;
1146}
1147
1148
1149void AliGenPythia::FinishRun()
1150{
1151// Print x-section summary
1152 fPythia->Pystat(1);
1153
1154 if (fNev > 0.) {
1155 fQ /= fNev;
1156 fX1 /= fNev;
1157 fX2 /= fNev;
1158 }
1159
1160 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1161 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1162}
1163
1164void AliGenPythia::AdjustWeights() const
1165{
1166// Adjust the weights after generation of all events
1167//
1168 if (gAlice) {
1169 TParticle *part;
1170 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1171 for (Int_t i=0; i<ntrack; i++) {
1172 part= gAlice->GetMCApp()->Particle(i);
1173 part->SetWeight(part->GetWeight()*fKineBias);
1174 }
1175 }
1176}
1177
1178void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1179{
1180// Treat protons as inside nuclei with mass numbers a1 and a2
1181
1182 fAProjectile = a1;
1183 fATarget = a2;
1184 fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
1185 fSetNuclei = kTRUE;
1186}
1187
1188
1189void AliGenPythia::MakeHeader()
1190{
1191//
1192// Make header for the simulated event
1193//
1194 if (gAlice) {
1195 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1196 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1197 }
1198
1199// Builds the event header, to be called after each event
1200 if (fHeader) delete fHeader;
1201 fHeader = new AliGenPythiaEventHeader("Pythia");
1202//
1203// Event type
1204 if(fProcDiff>0){
1205 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1206 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1207 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1208 }
1209 else
1210 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1211//
1212// Number of trials
1213 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1214//
1215// Event Vertex
1216 fHeader->SetPrimaryVertex(fVertex);
1217 fHeader->SetInteractionTime(fTime+fEventTime);
1218//
1219// Number of primaries
1220 fHeader->SetNProduced(fNprimaries);
1221//
1222// Jets that have triggered
1223
1224 //Need to store jets for b-jet studies too!
1225 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
1226 {
1227 Int_t ntrig, njet;
1228 Float_t jets[4][10];
1229 GetJets(njet, ntrig, jets);
1230
1231
1232 for (Int_t i = 0; i < ntrig; i++) {
1233 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1234 jets[3][i]);
1235 }
1236 }
1237//
1238// Copy relevant information from external header, if present.
1239//
1240 Float_t uqJet[4];
1241
1242 if (fRL) {
1243 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1244 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1245 {
1246 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1247
1248
1249 exHeader->TriggerJet(i, uqJet);
1250 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1251 }
1252 }
1253//
1254// Store quenching parameters
1255//
1256 if (fQuench){
1257 Double_t z[4] = {0.};
1258 Double_t xp = 0.;
1259 Double_t yp = 0.;
1260
1261 if (fQuench == 1) {
1262 // Pythia::Quench()
1263 fPythia->GetQuenchingParameters(xp, yp, z);
1264 } else if (fQuench == 2){
1265 // Pyquen
1266 Double_t r1 = PARIMP.rb1;
1267 Double_t r2 = PARIMP.rb2;
1268 Double_t b = PARIMP.b1;
1269 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1270 Double_t phi = PARIMP.psib1;
1271 xp = r * TMath::Cos(phi);
1272 yp = r * TMath::Sin(phi);
1273
1274 } else if (fQuench == 4) {
1275 // QPythia
1276 Double_t xy[2];
1277 Double_t i0i1[2];
1278 AliFastGlauber::Instance()->GetSavedXY(xy);
1279 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1280 xp = xy[0];
1281 yp = xy[1];
1282 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1283 }
1284
1285 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1286 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1287 }
1288//
1289// Store pt^hard
1290 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1291//
1292// Pass header
1293//
1294 AddHeader(fHeader);
1295 fHeader = 0x0;
1296}
1297
1298Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1299{
1300// Check the kinematic trigger condition
1301//
1302 Double_t eta[2];
1303 eta[0] = jet1->Eta();
1304 eta[1] = jet2->Eta();
1305 Double_t phi[2];
1306 phi[0] = jet1->Phi();
1307 phi[1] = jet2->Phi();
1308 Int_t pdg[2];
1309 pdg[0] = jet1->GetPdgCode();
1310 pdg[1] = jet2->GetPdgCode();
1311 Bool_t triggered = kFALSE;
1312
1313 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
1314 Int_t njets = 0;
1315 Int_t ntrig = 0;
1316 Float_t jets[4][10];
1317//
1318// Use Pythia clustering on parton level to determine jet axis
1319//
1320 GetJets(njets, ntrig, jets);
1321
1322 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1323//
1324 } else {
1325 Int_t ij = 0;
1326 Int_t ig = 1;
1327 if (pdg[0] == kGamma) {
1328 ij = 1;
1329 ig = 0;
1330 }
1331 //Check eta range first...
1332 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1333 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1334 {
1335 //Eta is okay, now check phi range
1336 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1337 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1338 {
1339 triggered = kTRUE;
1340 }
1341 }
1342 }
1343 return triggered;
1344}
1345
1346
1347
1348Bool_t AliGenPythia::CheckKinematicsOnChild(){
1349//
1350//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1351//
1352 Bool_t checking = kFALSE;
1353 Int_t j, kcode, ks, km;
1354 Int_t nPartAcc = 0; //number of particles in the acceptance range
1355 Int_t numberOfAcceptedParticles = 1;
1356 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1357 Int_t npart = fParticles.GetEntriesFast();
1358
1359 for (j = 0; j<npart; j++) {
1360 TParticle * jparticle = (TParticle *) fParticles.At(j);
1361 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1362 ks = jparticle->GetStatusCode();
1363 km = jparticle->GetFirstMother();
1364
1365 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1366 nPartAcc++;
1367 }
1368 if( numberOfAcceptedParticles <= nPartAcc){
1369 checking = kTRUE;
1370 break;
1371 }
1372 }
1373
1374 return checking;
1375}
1376
1377void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1378{
1379 //
1380 // Load event into Pythia Common Block
1381 //
1382
1383 Int_t npart = stack -> GetNprimary();
1384 Int_t n0 = 0;
1385
1386 if (!flag) {
1387 (fPythia->GetPyjets())->N = npart;
1388 } else {
1389 n0 = (fPythia->GetPyjets())->N;
1390 (fPythia->GetPyjets())->N = n0 + npart;
1391 }
1392
1393
1394 for (Int_t part = 0; part < npart; part++) {
1395 TParticle *mPart = stack->Particle(part);
1396
1397 Int_t kf = mPart->GetPdgCode();
1398 Int_t ks = mPart->GetStatusCode();
1399 Int_t idf = mPart->GetFirstDaughter();
1400 Int_t idl = mPart->GetLastDaughter();
1401
1402 if (reHadr) {
1403 if (ks == 11 || ks == 12) {
1404 ks -= 10;
1405 idf = -1;
1406 idl = -1;
1407 }
1408 }
1409
1410 Float_t px = mPart->Px();
1411 Float_t py = mPart->Py();
1412 Float_t pz = mPart->Pz();
1413 Float_t e = mPart->Energy();
1414 Float_t m = mPart->GetCalcMass();
1415
1416
1417 (fPythia->GetPyjets())->P[0][part+n0] = px;
1418 (fPythia->GetPyjets())->P[1][part+n0] = py;
1419 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1420 (fPythia->GetPyjets())->P[3][part+n0] = e;
1421 (fPythia->GetPyjets())->P[4][part+n0] = m;
1422
1423 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1424 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1425 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1426 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1427 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1428 }
1429}
1430
1431void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1432{
1433 //
1434 // Load event into Pythia Common Block
1435 //
1436
1437 Int_t npart = stack -> GetEntries();
1438 Int_t n0 = 0;
1439
1440 if (!flag) {
1441 (fPythia->GetPyjets())->N = npart;
1442 } else {
1443 n0 = (fPythia->GetPyjets())->N;
1444 (fPythia->GetPyjets())->N = n0 + npart;
1445 }
1446
1447
1448 for (Int_t part = 0; part < npart; part++) {
1449 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1450 if (!mPart) continue;
1451
1452 Int_t kf = mPart->GetPdgCode();
1453 Int_t ks = mPart->GetStatusCode();
1454 Int_t idf = mPart->GetFirstDaughter();
1455 Int_t idl = mPart->GetLastDaughter();
1456
1457 if (reHadr) {
1458 if (ks == 11 || ks == 12) {
1459 ks -= 10;
1460 idf = -1;
1461 idl = -1;
1462 }
1463 }
1464
1465 Float_t px = mPart->Px();
1466 Float_t py = mPart->Py();
1467 Float_t pz = mPart->Pz();
1468 Float_t e = mPart->Energy();
1469 Float_t m = mPart->GetCalcMass();
1470
1471
1472 (fPythia->GetPyjets())->P[0][part+n0] = px;
1473 (fPythia->GetPyjets())->P[1][part+n0] = py;
1474 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1475 (fPythia->GetPyjets())->P[3][part+n0] = e;
1476 (fPythia->GetPyjets())->P[4][part+n0] = m;
1477
1478 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1479 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1480 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1481 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1482 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1483 }
1484}
1485
1486
1487void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1488{
1489//
1490// Calls the Pythia jet finding algorithm to find jets in the current event
1491//
1492//
1493//
1494// Save jets
1495 Int_t n = fPythia->GetN();
1496
1497//
1498// Run Jet Finder
1499 fPythia->Pycell(njets);
1500 Int_t i;
1501 for (i = 0; i < njets; i++) {
1502 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1503 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1504 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1505 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1506
1507 jets[0][i] = px;
1508 jets[1][i] = py;
1509 jets[2][i] = pz;
1510 jets[3][i] = e;
1511 }
1512}
1513
1514
1515
1516void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1517{
1518//
1519// Calls the Pythia clustering algorithm to find jets in the current event
1520//
1521 Int_t n = fPythia->GetN();
1522 nJets = 0;
1523 nJetsTrig = 0;
1524 if (fJetReconstruction == kCluster) {
1525//
1526// Configure cluster algorithm
1527//
1528 fPythia->SetPARU(43, 2.);
1529 fPythia->SetMSTU(41, 1);
1530//
1531// Call cluster algorithm
1532//
1533 fPythia->Pyclus(nJets);
1534//
1535// Loading jets from common block
1536//
1537 } else {
1538
1539//
1540// Run Jet Finder
1541 fPythia->Pycell(nJets);
1542 }
1543
1544 Int_t i;
1545 for (i = 0; i < nJets; i++) {
1546 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1547 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1548 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1549 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1550 Float_t pt = TMath::Sqrt(px * px + py * py);
1551 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
1552 Float_t theta = TMath::ATan2(pt,pz);
1553 Float_t et = e * TMath::Sin(theta);
1554 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
1555 if (
1556 eta > fEtaMinJet && eta < fEtaMaxJet &&
1557 phi > fPhiMinJet && phi < fPhiMaxJet &&
1558 et > fEtMinJet && et < fEtMaxJet
1559 )
1560 {
1561 jets[0][nJetsTrig] = px;
1562 jets[1][nJetsTrig] = py;
1563 jets[2][nJetsTrig] = pz;
1564 jets[3][nJetsTrig] = e;
1565 nJetsTrig++;
1566// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1567 } else {
1568// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1569 }
1570 }
1571}
1572
1573void AliGenPythia::GetSubEventTime()
1574{
1575 // Calculates time of the next subevent
1576 fEventTime = 0.;
1577 if (fEventsTime) {
1578 TArrayF &array = *fEventsTime;
1579 fEventTime = array[fCurSubEvent++];
1580 }
1581 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1582 return;
1583}
1584
1585Bool_t AliGenPythia::IsInBarrel(const Float_t eta) const
1586{
1587 // Is particle in Central Barrel acceptance?
1588 // etamin=-etamax
1589 if( eta < fTriggerEta )
1590 return kTRUE;
1591 else
1592 return kFALSE;
1593}
1594
1595Bool_t AliGenPythia::IsInEMCAL(const Float_t phi, const Float_t eta) const
1596{
1597 // Is particle in EMCAL acceptance?
1598 // phi in degrees, etamin=-etamax
1599 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1600 eta < fEMCALEta )
1601 return kTRUE;
1602 else
1603 return kFALSE;
1604}
1605
1606Bool_t AliGenPythia::IsInPHOS(const Float_t phi, const Float_t eta, const Int_t iparticle)
1607{
1608 // Is particle in PHOS acceptance?
1609 // Acceptance slightly larger considered.
1610 // phi in degrees, etamin=-etamax
1611 // iparticle is the index of the particle to be checked, for PHOS rotation case
1612
1613 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1614 eta < fPHOSEta )
1615 return kTRUE;
1616 else
1617 {
1618 if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1619
1620 return kFALSE;
1621 }
1622}
1623
1624void AliGenPythia::RotatePhi(Bool_t& okdd)
1625{
1626 //Rotate event in phi to enhance events in PHOS acceptance
1627
1628 if(fPHOSRotateCandidate < 0) return ;
1629
1630 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1631 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1632 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1633 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
1634
1635 //calculate deltaphi
1636 TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
1637 Double_t phphi = ph->Phi();
1638 Double_t deltaphi = phiPHOS - phphi;
1639
1640
1641
1642 //loop for all particles and produce the phi rotation
1643 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1644 Double_t oldphi, newphi;
1645 Double_t newVx, newVy, r, vZ, time;
1646 Double_t newPx, newPy, pt, pz, e;
1647 for(Int_t i=0; i< np; i++) {
1648 TParticle* iparticle = (TParticle *) fParticles.At(i);
1649 oldphi = iparticle->Phi();
1650 newphi = oldphi + deltaphi;
1651 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1652 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1653
1654 r = iparticle->R();
1655 newVx = r * TMath::Cos(newphi);
1656 newVy = r * TMath::Sin(newphi);
1657 vZ = iparticle->Vz(); // don't transform
1658 time = iparticle->T(); // don't transform
1659
1660 pt = iparticle->Pt();
1661 newPx = pt * TMath::Cos(newphi);
1662 newPy = pt * TMath::Sin(newphi);
1663 pz = iparticle->Pz(); // don't transform
1664 e = iparticle->Energy(); // don't transform
1665
1666 // apply rotation
1667 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1668 iparticle->SetMomentum(newPx, newPy, pz, e);
1669
1670 } //end particle loop
1671
1672 // now let's check that we put correctly the candidate photon in PHOS
1673 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1674 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1675 if(IsInPHOS(phi,eta,-1))
1676 okdd = kTRUE;
1677
1678 // reset the value for next event
1679 fPHOSRotateCandidate = -1;
1680
1681}
1682
1683
1684Bool_t AliGenPythia::CheckDiffraction()
1685{
1686 // use this method only with Perugia-0 tune!
1687
1688 // printf("AAA\n");
1689
1690 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1691
1692 Int_t iPart1=-1;
1693 Int_t iPart2=-1;
1694
1695 Double_t y1 = 1e10;
1696 Double_t y2 = -1e10;
1697
1698 const Int_t kNstable=20;
1699 const Int_t pdgStable[20] = {
1700 22, // Photon
1701 11, // Electron
1702 12, // Electron Neutrino
1703 13, // Muon
1704 14, // Muon Neutrino
1705 15, // Tau
1706 16, // Tau Neutrino
1707 211, // Pion
1708 321, // Kaon
1709 311, // K0
1710 130, // K0s
1711 310, // K0l
1712 2212, // Proton
1713 2112, // Neutron
1714 3122, // Lambda_0
1715 3112, // Sigma Minus
1716 3222, // Sigma Plus
1717 3312, // Xsi Minus
1718 3322, // Xsi0
1719 3334 // Omega
1720 };
1721
1722 for (Int_t i = 0; i < np; i++) {
1723 TParticle * part = (TParticle *) fParticles.At(i);
1724
1725 Int_t statusCode = part->GetStatusCode();
1726
1727 // Initial state particle
1728 if (statusCode != 1)
1729 continue;
1730
1731 Int_t pdg = TMath::Abs(part->GetPdgCode());
1732 Bool_t isStable = kFALSE;
1733 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1734 if (pdg == pdgStable[i1]) {
1735 isStable = kTRUE;
1736 break;
1737 }
1738 }
1739 if(!isStable)
1740 continue;
1741
1742 Double_t y = part->Y();
1743
1744 if (y < y1)
1745 {
1746 y1 = y;
1747 iPart1 = i;
1748 }
1749 if (y > y2)
1750 {
1751 y2 = y;
1752 iPart2 = i;
1753 }
1754 }
1755
1756 if(iPart1<0 || iPart2<0) return kFALSE;
1757
1758 y1=TMath::Abs(y1);
1759 y2=TMath::Abs(y2);
1760
1761 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1762 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1763
1764 Int_t pdg1 = part1->GetPdgCode();
1765 Int_t pdg2 = part2->GetPdgCode();
1766
1767
1768 Int_t iPart = -1;
1769 if (pdg1 == 2212 && pdg2 == 2212)
1770 {
1771 if(y1 > y2)
1772 iPart = iPart1;
1773 else if(y1 < y2)
1774 iPart = iPart2;
1775 else {
1776 iPart = iPart1;
1777 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1778 }
1779 }
1780 else if (pdg1 == 2212)
1781 iPart = iPart1;
1782 else if (pdg2 == 2212)
1783 iPart = iPart2;
1784
1785
1786
1787
1788
1789 Double_t M=-1.;
1790 if(iPart>0) {
1791 TParticle * part = (TParticle *) fParticles.At(iPart);
1792 Double_t E= part->Energy();
1793 Double_t P= part->P();
1794 M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
1795 }
1796
1797 Double_t Mmin, Mmax, wSD, wDD, wND;
1798 if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1799
1800 if(M>-1 && M<Mmin) return kFALSE;
1801 if(M>Mmax) M=-1;
1802
1803 Int_t procType=fPythia->GetMSTI(1);
1804 Int_t proc0=2;
1805 if(procType== 94) proc0=1;
1806 if(procType== 92 || procType== 93) proc0=0;
1807
1808 Int_t proc=2;
1809 if(M>0) proc=0;
1810 else if(proc0==1) proc=1;
1811
1812 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1813 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
1814
1815
1816 // if(proc==1 || proc==2) return kFALSE;
1817
1818 if(proc!=0) {
1819 if(proc0!=0) fProcDiff = procType;
1820 else fProcDiff = 95;
1821 return kTRUE;
1822 }
1823
1824 if(wSD<0) AliError("wSD<0 ! \n");
1825
1826 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1827
1828 // printf("iPart = %d\n", iPart);
1829
1830 if(iPart==iPart1) fProcDiff=93;
1831 else if(iPart==iPart2) fProcDiff=92;
1832 else {
1833 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1834
1835 }
1836
1837 return kTRUE;
1838}
1839
1840
1841
1842Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1843 Double_t &wSD, Double_t &wDD, Double_t &wND)
1844{
1845
1846 // 900 GeV
1847 if(TMath::Abs(fEnergyCMS-900)<1 ){
1848
1849const Int_t nbin=400;
1850Double_t bin[]={
18511.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
18524.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
18537.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
185410.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
185513.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
185615.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
185718.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
185821.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
185924.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
186027.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
186130.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
186233.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
186336.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
186439.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
186542.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
186645.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
186748.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
186851.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
186954.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
187057.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
187160.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
187263.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
187366.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
187469.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
187572.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
187675.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
187778.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
187881.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
187984.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
188087.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
188190.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
188293.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
188396.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
188499.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1885102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1886105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1887108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1888111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1889114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1890117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1891120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1892123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1893126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
1894129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
1895132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
1896135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
1897138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
1898141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
1899144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
1900147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
1901150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
1902153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
1903156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
1904159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
1905162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
1906165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
1907168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
1908171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
1909174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
1910177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
1911180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
1912183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
1913186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
1914189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
1915192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
1916195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
1917198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1918Double_t w[]={
19191.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
19200.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
19210.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
19220.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
19230.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
19240.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
19250.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
19260.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
19270.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
19280.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
19290.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
19300.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
19310.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
19320.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
19330.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
19340.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
19350.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
19360.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
19370.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
19380.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
19390.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
19400.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
19410.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
19420.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
19430.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
19440.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
19450.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
19460.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
19470.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
19480.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
19490.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
19500.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
19510.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
19520.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
19530.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
19540.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
19550.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
19560.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
19570.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
19580.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
19590.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
19600.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
19610.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
19620.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
19630.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
19640.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
19650.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
19660.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
19670.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
19680.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
19690.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
19700.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
19710.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
19720.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
19730.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
19740.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
19750.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
19760.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
19770.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
19780.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
19790.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
19800.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
19810.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
19820.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
19830.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
19840.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
19850.386112, 0.364314, 0.434375, 0.334629};
1986wDD = 0.379611;
1987wND = 0.496961;
1988wSD = -1;
1989
1990 Mmin = bin[0];
1991 Mmax = bin[nbin];
1992 if(M<Mmin || M>Mmax) return kTRUE;
1993
1994 Int_t ibin=nbin-1;
1995 for(Int_t i=1; i<=nbin; i++)
1996 if(M<=bin[i]) {
1997 ibin=i-1;
1998 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
1999 break;
2000 }
2001 wSD=w[ibin];
2002 return kTRUE;
2003 }
2004 else if(TMath::Abs(fEnergyCMS-2760)<1 ){
2005
2006const Int_t nbin=400;
2007Double_t bin[]={
20081.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
20094.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
20107.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
201110.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
201213.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
201315.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
201418.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
201521.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
201624.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
201727.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
201830.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
201933.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
202036.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
202139.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
202242.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
202345.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
202448.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
202551.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
202654.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
202757.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
202860.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
202963.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
203066.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
203169.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
203272.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
203375.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
203478.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
203581.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
203684.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
203787.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
203890.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
203993.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
204096.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
204199.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2042102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2043105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2044108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2045111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2046114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2047117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2048120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2049123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2050126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2051129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2052132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2053135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2054138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2055141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2056144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2057147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2058150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2059153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2060156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2061159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2062162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2063165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2064168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2065171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2066174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2067177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2068180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2069183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2070186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2071189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2072192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2073195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2074198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2075Double_t w[]={
20761.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
20770.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
20780.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
20790.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
20800.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
20810.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
20820.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
20830.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
20840.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
20850.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
20860.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
20870.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
20880.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
20890.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
20900.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
20910.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
20920.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
20930.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
20940.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
20950.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
20960.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
20970.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
20980.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
20990.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
21000.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
21010.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
21020.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
21030.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
21040.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
21050.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
21060.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
21070.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
21080.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
21090.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
21100.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
21110.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
21120.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
21130.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
21140.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
21150.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
21160.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
21170.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
21180.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
21190.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
21200.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
21210.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
21220.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
21230.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
21240.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
21250.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
21260.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
21270.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
21280.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
21290.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
21300.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
21310.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
21320.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
21330.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
21340.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
21350.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
21360.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
21370.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
21380.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
21390.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
21400.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
21410.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
21420.201712, 0.242225, 0.255565, 0.258738};
2143wDD = 0.512813;
2144wND = 0.518820;
2145wSD = -1;
2146
2147 Mmin = bin[0];
2148 Mmax = bin[nbin];
2149 if(M<Mmin || M>Mmax) return kTRUE;
2150
2151 Int_t ibin=nbin-1;
2152 for(Int_t i=1; i<=nbin; i++)
2153 if(M<=bin[i]) {
2154 ibin=i-1;
2155 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2156 break;
2157 }
2158 wSD=w[ibin];
2159 return kTRUE;
2160 }
2161
2162
2163 else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2164const Int_t nbin=400;
2165Double_t bin[]={
21661.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
21674.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
21687.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
216910.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
217013.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
217115.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
217218.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
217321.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
217424.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
217527.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
217630.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
217733.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
217836.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
217939.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
218042.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
218145.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
218248.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
218351.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
218454.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
218557.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
218660.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
218763.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
218866.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
218969.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
219072.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
219175.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
219278.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
219381.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
219484.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
219587.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
219690.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
219793.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
219896.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
219999.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2200102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2201105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2202108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2203111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2204114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2205117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2206120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2207123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2208126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2209129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2210132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2211135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2212138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2213141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2214144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2215147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2216150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2217153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2218156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2219159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2220162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2221165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2222168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2223171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2224174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2225177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2226180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2227183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2228186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2229189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2230192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2231195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2232198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2233Double_t w[]={
22341.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
22350.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
22360.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
22370.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
22380.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
22390.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
22400.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
22410.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
22420.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
22430.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
22440.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
22450.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
22460.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
22470.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
22480.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
22490.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
22500.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
22510.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
22520.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
22530.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
22540.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
22550.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
22560.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
22570.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
22580.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
22590.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
22600.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
22610.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
22620.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
22630.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
22640.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
22650.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
22660.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
22670.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
22680.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
22690.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
22700.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
22710.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
22720.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
22730.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
22740.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
22750.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
22760.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
22770.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
22780.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
22790.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
22800.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
22810.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
22820.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
22830.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
22840.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
22850.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
22860.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
22870.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
22880.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
22890.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
22900.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
22910.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
22920.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
22930.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
22940.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
22950.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
22960.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
22970.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
22980.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
22990.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
23000.175006, 0.223482, 0.202706, 0.218108};
2301wDD = 0.207705;
2302wND = 0.289628;
2303wSD = -1;
2304
2305 Mmin = bin[0];
2306 Mmax = bin[nbin];
2307
2308 if(M<Mmin || M>Mmax) return kTRUE;
2309
2310 Int_t ibin=nbin-1;
2311 for(Int_t i=1; i<=nbin; i++)
2312 if(M<=bin[i]) {
2313 ibin=i-1;
2314 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2315 break;
2316 }
2317 wSD=w[ibin];
2318 return kTRUE;
2319 }
2320
2321 return kFALSE;
2322}
2323
2324Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2325{
2326// Check if this is a heavy flavor decay product
2327 TParticle * part = (TParticle *) fParticles.At(ipart);
2328 Int_t mpdg = TMath::Abs(part->GetPdgCode());
2329 Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2330 //
2331 // Light hadron
2332 if (mfl >= 4 && mfl < 6) return kTRUE;
2333
2334 Int_t imo = part->GetFirstMother()-1;
2335 TParticle* pm = part;
2336 //
2337 // Heavy flavor hadron produced by generator
2338 while (imo > -1) {
2339 pm = (TParticle*)fParticles.At(imo);
2340 mpdg = TMath::Abs(pm->GetPdgCode());
2341 mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2342 if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2343 imo = pm->GetFirstMother()-1;
2344 }
2345 return kFALSE;
2346}
2347
2348Bool_t AliGenPythia::CheckDetectorAcceptance(const Float_t phi, const Float_t eta, const Int_t iparticle)
2349{
2350 // check the eta/phi correspond to the detectors acceptance
2351 // iparticle is the index of the particle to be checked, for PHOS rotation case
2352 if (fCheckPHOS && IsInPHOS (phi,eta,iparticle)) return kTRUE;
2353 else if(fCheckEMCAL && IsInEMCAL (phi,eta)) return kTRUE;
2354 else if(fCheckBarrel && IsInBarrel( eta)) return kTRUE;
2355 else return kFALSE;
2356}
2357
2358Bool_t AliGenPythia::TriggerOnSelectedParticles(const Int_t np)
2359{
2360 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2361 // implemented primaryly for kPyJets, but extended to any kind of process.
2362
2363 //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2364 // fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
2365
2366 Bool_t ok = kFALSE;
2367 for (Int_t i=0; i< np; i++) {
2368
2369 TParticle* iparticle = (TParticle *) fParticles.At(i);
2370
2371 Int_t pdg = iparticle->GetPdgCode();
2372 Int_t status = iparticle->GetStatusCode();
2373 Int_t imother = iparticle->GetFirstMother() - 1;
2374
2375 TParticle* pmother = 0x0;
2376 Int_t momStatus = -1;
2377 Int_t momPdg = -1;
2378 if(imother > 0 ){
2379 pmother = (TParticle *) fParticles.At(imother);
2380 momStatus = pmother->GetStatusCode();
2381 momPdg = pmother->GetPdgCode();
2382 }
2383
2384 ok = kFALSE;
2385
2386 //
2387 // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2388 //
2389 // Hadron
2390 if (fHadronInCalo && status == 1)
2391 {
2392 if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
2393 // (in case neutral mesons were declared stable)
2394 ok = kTRUE;
2395 }
2396 //Electron
2397 else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2398 {
2399 ok = kTRUE;
2400 }
2401 //Fragmentation photon
2402 else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2403 {
2404 if(momStatus != 11) ok = kTRUE ; // No photon from hadron decay
2405 }
2406 // Decay photon
2407 else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2408 {
2409 if( momStatus == 11)
2410 {
2411 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2412 // pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2413 ok = kTRUE ; // photon from hadron decay
2414
2415 //In case only decays from pi0 or eta requested
2416 if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2417 if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2418 }
2419
2420 }
2421 // Pi0 or Eta particle
2422 else if ((fPi0InCalo || fEtaInCalo))
2423 {
2424 if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2425
2426 if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2427 {
2428 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2429 ok = kTRUE;
2430 }
2431 else if (fEtaInCalo && pdg == 221)
2432 {
2433 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2434 ok = kTRUE;
2435 }
2436
2437 }// pi0 or eta
2438
2439 //
2440 // Check that the selected particle is in the calorimeter acceptance
2441 //
2442 if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2443 {
2444 //Just check if the selected particle falls in the acceptance
2445 if(!fForceNeutralMeson2PhotonDecay )
2446 {
2447 //printf("\t Check acceptance! \n");
2448 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2449 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2450
2451 if(CheckDetectorAcceptance(phi,eta,i))
2452 {
2453 ok =kTRUE;
2454 AliDebug(1,Form("Selected trigger pdg %d, status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2455 //printf("\t Accept \n");
2456 break;
2457 }
2458 else ok = kFALSE;
2459 }
2460 //Mesons have several decay modes, select only those decaying into 2 photons
2461 else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2462 {
2463 // In case we want the pi0/eta trigger,
2464 // check the decay mode (2 photons)
2465
2466 //printf("\t Force decay 2 gamma\n");
2467
2468 Int_t ndaughters = iparticle->GetNDaughters();
2469 if(ndaughters != 2){
2470 ok=kFALSE;
2471 continue;
2472 }
2473
2474 TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2475 TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2476 if(!d1 || !d2) {
2477 ok=kFALSE;
2478 continue;
2479 }
2480
2481 //iparticle->Print();
2482 //d1->Print();
2483 //d2->Print();
2484
2485 Int_t pdgD1 = d1->GetPdgCode();
2486 Int_t pdgD2 = d2->GetPdgCode();
2487 //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2488 //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2489
2490 if(pdgD1 != 22 || pdgD2 != 22){
2491 ok = kFALSE;
2492 continue;
2493 }
2494
2495 //printf("\t accept decay\n");
2496
2497 //Trigger on the meson, not on the daughter
2498 if(!fDecayPhotonInCalo){
2499
2500 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2501 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2502
2503 if(CheckDetectorAcceptance(phi,eta,i))
2504 {
2505 //printf("\t Accept meson pdg %d\n",pdg);
2506 ok =kTRUE;
2507 AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2508 break;
2509 } else {
2510 ok=kFALSE;
2511 continue;
2512 }
2513 }
2514
2515 //printf("Check daughters acceptance\n");
2516
2517 //Trigger on the meson daughters
2518 //Photon 1
2519 Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2520 Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
2521 if(d1->Pt() > fTriggerParticleMinPt)
2522 {
2523 //printf("\t Check acceptance photon 1! \n");
2524 if(CheckDetectorAcceptance(phi,eta,i))
2525 {
2526 //printf("\t Accept Photon 1\n");
2527 ok =kTRUE;
2528 AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2529 break;
2530 }
2531 else ok = kFALSE;
2532 } // pt cut
2533 else ok = kFALSE;
2534
2535 //Photon 2
2536 phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2537 eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2538
2539 if(d2->Pt() > fTriggerParticleMinPt)
2540 {
2541 //printf("\t Check acceptance photon 2! \n");
2542 if(CheckDetectorAcceptance(phi,eta,i))
2543 {
2544 //printf("\t Accept Photon 2\n");
2545 ok =kTRUE;
2546 AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2547 break;
2548 }
2549 else ok = kFALSE;
2550 } // pt cut
2551 else ok = kFALSE;
2552 } // force 2 photon daughters in pi0/eta decays
2553 else ok = kFALSE;
2554 } else ok = kFALSE; // check acceptance
2555 } // primary loop
2556
2557 //
2558 // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2559 // A particle passing all trigger conditions except phi position in PHOS, is used as reference
2560 //
2561 if(fCheckPHOSeta)
2562 {
2563 RotatePhi(ok);
2564 }
2565
2566 return ok;
2567}
2568