1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 ...)
25 // andreas.morsch@cern.ch
28 #include <TClonesArray.h>
29 #include <TDatabasePDG.h>
30 #include <TParticle.h>
32 #include <TObjArray.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"
45 #include "AliRunLoader.h"
48 #include "PyquenCommon.h"
51 ClassImp(AliGenPythia)
54 AliGenPythia::AliGenPythia():
96 fDecayer(new AliDecayerPythia()),
104 fPhiMaxJet(2.* TMath::Pi()),
105 fJetReconstruction(kCell),
109 fPhiMaxGamma(2. * TMath::Pi()),
113 fPycellThreshold(0.),
115 fPycellMinEtJet(10.),
116 fPycellMaxRadius(1.),
117 fStackFillOpt(kFlavorSelection),
119 fFragmentation(kTRUE),
121 fUseNuclearPDF(kFALSE),
122 fUseLorentzBoost(kTRUE),
130 fTriggerMultiplicity(0),
131 fTriggerMultiplicityEta(0),
132 fTriggerMultiplicityPtMin(0),
133 fCountMode(kCountAll),
138 fFragPhotonInCalo(kFALSE),
139 fHadronInCalo(kFALSE) ,
142 fPhotonInCalo(kFALSE), // not in use
143 fDecayPhotonInCalo(kFALSE),
144 fForceNeutralMeson2PhotonDecay(kFALSE),
146 fEleInEMCAL(kFALSE), // not in use
147 fCheckBarrel(kFALSE),
150 fCheckPHOSeta(kFALSE),
151 fPHOSRotateCandidate(-1),
152 fTriggerParticleMinPt(0),
153 fPhotonMinPt(0), // not in use
154 fElectronMinPt(0), // not in use
164 // Default Constructor
166 if (!AliPythiaRndm::GetPythiaRandom())
167 AliPythiaRndm::SetPythiaRandom(GetRandom());
170 AliGenPythia::AliGenPythia(Int_t npart)
182 fInteractionRate(0.),
196 fHadronisation(kTRUE),
197 fPatchOmegaDalitz(0),
199 fReadFromFile(kFALSE),
212 fDecayer(new AliDecayerPythia()),
213 fDebugEventFirst(-1),
220 fPhiMaxJet(2.* TMath::Pi()),
221 fJetReconstruction(kCell),
225 fPhiMaxGamma(2. * TMath::Pi()),
229 fPycellThreshold(0.),
231 fPycellMinEtJet(10.),
232 fPycellMaxRadius(1.),
233 fStackFillOpt(kFlavorSelection),
235 fFragmentation(kTRUE),
237 fUseNuclearPDF(kFALSE),
238 fUseLorentzBoost(kTRUE),
246 fTriggerMultiplicity(0),
247 fTriggerMultiplicityEta(0),
248 fTriggerMultiplicityPtMin(0),
249 fCountMode(kCountAll),
254 fFragPhotonInCalo(kFALSE),
255 fHadronInCalo(kFALSE) ,
258 fPhotonInCalo(kFALSE), // not in use
259 fDecayPhotonInCalo(kFALSE),
260 fForceNeutralMeson2PhotonDecay(kFALSE),
262 fEleInEMCAL(kFALSE), // not in use
263 fCheckBarrel(kFALSE),
266 fCheckPHOSeta(kFALSE),
267 fPHOSRotateCandidate(-1),
268 fTriggerParticleMinPt(0),
269 fPhotonMinPt(0), // not in use
270 fElectronMinPt(0), // not in use
280 // default charm production at 5. 5 TeV
282 // structure function GRVHO
286 fTitle= "Particle Generator using PYTHIA";
288 // Set random number generator
289 if (!AliPythiaRndm::GetPythiaRandom())
290 AliPythiaRndm::SetPythiaRandom(GetRandom());
293 AliGenPythia::~AliGenPythia()
296 if(fEventsTime) delete fEventsTime;
299 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
301 // Generate pileup using user specified rate
302 fInteractionRate = rate;
303 fTimeWindow = timewindow;
307 void AliGenPythia::GeneratePileup()
309 // Generate sub events time for pileup
311 if(fInteractionRate == 0.) {
312 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
316 Int_t npart = NumberParticles();
318 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
322 if(fEventsTime) delete fEventsTime;
323 fEventsTime = new TArrayF(npart);
324 TArrayF &array = *fEventsTime;
325 for(Int_t ipart = 0; ipart < npart; ipart++)
328 Float_t eventtime = 0.;
331 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
332 if(eventtime > fTimeWindow) break;
333 array.Set(array.GetSize()+1);
334 array[array.GetSize()-1] = eventtime;
340 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
341 if(TMath::Abs(eventtime) > fTimeWindow) break;
342 array.Set(array.GetSize()+1);
343 array[array.GetSize()-1] = eventtime;
346 SetNumberParticles(fEventsTime->GetSize());
349 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
350 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
352 // Set pycell parameters
353 fPycellEtaMax = etamax;
356 fPycellThreshold = thresh;
357 fPycellEtSeed = etseed;
358 fPycellMinEtJet = minet;
359 fPycellMaxRadius = r;
364 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
366 // Set a range of event numbers, for which a table
367 // of generated particle will be printed
368 fDebugEventFirst = eventFirst;
369 fDebugEventLast = eventLast;
370 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
373 void AliGenPythia::Init()
377 SetMC(AliPythia::Instance());
378 fPythia=(AliPythia*) fMCEvGen;
381 fParentWeight=1./Float_t(fNpart);
385 fPythia->SetCKIN(3,fPtHardMin);
386 fPythia->SetCKIN(4,fPtHardMax);
387 fPythia->SetCKIN(7,fYHardMin);
388 fPythia->SetCKIN(8,fYHardMax);
390 if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
393 fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
395 if (fFragmentation) {
396 fPythia->SetMSTP(111,1);
398 fPythia->SetMSTP(111,0);
402 // initial state radiation
403 fPythia->SetMSTP(61,fGinit);
404 // final state radiation
405 fPythia->SetMSTP(71,fGfinal);
408 fPythia->SetMSTP(91,1);
409 fPythia->SetPARP(91,fPtKick);
410 fPythia->SetPARP(93, 4. * fPtKick);
412 fPythia->SetMSTP(91,0);
415 if (fReadLHEF) fPythia->OpenFortranFile(97, const_cast<char*>(fkNameLHEF));
418 fRL = AliRunLoader::Open(fkFileName, "Partons");
419 fRL->LoadKinematics();
425 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
426 // Forward Paramters to the AliPythia object
427 fDecayer->SetForceDecay(fForceDecay);
428 // Switch off Heavy Flavors on request
430 // Maximum number of quark flavours used in pdf
431 fPythia->SetMSTP(58, 3);
432 // Maximum number of flavors that can be used in showers
433 fPythia->SetMSTJ(45, 3);
434 // Switch off g->QQbar splitting in decay table
435 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
441 // Parent and Children Selection
444 case kPyOldUEQ2ordered:
445 case kPyOldUEQ2ordered2:
449 case kPyCharmUnforced:
450 case kPyCharmPbPbMNR:
453 case kPyCharmppMNRwmi:
455 fParentSelect[0] = 411;
456 fParentSelect[1] = 421;
457 fParentSelect[2] = 431;
458 fParentSelect[3] = 4122;
459 fParentSelect[4] = 4232;
460 fParentSelect[5] = 4132;
461 fParentSelect[6] = 4332;
467 fParentSelect[0] = 421;
470 case kPyDPlusPbPbMNR:
473 fParentSelect[0] = 411;
476 case kPyDPlusStrangePbPbMNR:
477 case kPyDPlusStrangepPbMNR:
478 case kPyDPlusStrangeppMNR:
479 fParentSelect[0] = 431;
482 case kPyLambdacppMNR:
483 fParentSelect[0] = 4122;
488 case kPyBeautyPbPbMNR:
489 case kPyBeautypPbMNR:
491 case kPyBeautyppMNRwmi:
493 fParentSelect[0]= 511;
494 fParentSelect[1]= 521;
495 fParentSelect[2]= 531;
496 fParentSelect[3]= 5122;
497 fParentSelect[4]= 5132;
498 fParentSelect[5]= 5232;
499 fParentSelect[6]= 5332;
502 case kPyBeautyUnforced:
503 fParentSelect[0] = 511;
504 fParentSelect[1] = 521;
505 fParentSelect[2] = 531;
506 fParentSelect[3] = 5122;
507 fParentSelect[4] = 5132;
508 fParentSelect[5] = 5232;
509 fParentSelect[6] = 5332;
514 fParentSelect[0] = 443;
517 case kPyMbAtlasTuneMC09:
519 case kPyMbWithDirectPhoton:
527 case kPyWPWHG: /// !!!!!!!!!! Change done for W prod with POWHEG !!!!!!!!!! :)
530 case kPyMBRSingleDiffraction:
531 case kPyMBRDoubleDiffraction:
532 case kPyMBRCentralDiffraction:
537 // JetFinder for Trigger
539 // Configure detector (EMCAL like)
541 fPythia->SetPARU(51, fPycellEtaMax);
542 fPythia->SetMSTU(51, fPycellNEta);
543 fPythia->SetMSTU(52, fPycellNPhi);
545 // Configure Jet Finder
547 fPythia->SetPARU(58, fPycellThreshold);
548 fPythia->SetPARU(52, fPycellEtSeed);
549 fPythia->SetPARU(53, fPycellMinEtJet);
550 fPythia->SetPARU(54, fPycellMaxRadius);
551 fPythia->SetMSTU(54, 2);
553 // This counts the total number of calls to Pyevnt() per run.
564 // Reset Lorentz boost if demanded
565 if(!fUseLorentzBoost) {
567 Warning("Init","Demand to discard Lorentz boost.\n");
574 Warning("Init","Deprecated function SetNuclei() used (nPDFs + no boost). Use SetProjectile + SetTarget + SetUseNuclearPDF + SetUseLorentzBoost instead.\n");
576 fPythia->SetPARJ(200, 0.0);
577 fPythia->SetPARJ(199, 0.0);
578 fPythia->SetPARJ(198, 0.0);
579 fPythia->SetPARJ(197, 0.0);
582 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
585 if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
588 // Nestor's change of the splittings
589 fPythia->SetPARJ(200, 0.8);
590 fPythia->SetMSTJ(41, 1); // QCD radiation only
591 fPythia->SetMSTJ(42, 2); // angular ordering
592 fPythia->SetMSTJ(44, 2); // option to run alpha_s
593 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
594 fPythia->SetMSTJ(50, 0); // No coherence in first branching
595 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
596 } else if (fQuench == 4) {
597 // Armesto-Cunqueiro-Salgado change of the splittings.
598 AliFastGlauber* glauber = AliFastGlauber::Instance();
600 //read and store transverse almonds corresponding to differnt
602 glauber->SetCentralityClass(0.,0.1);
603 fPythia->SetPARJ(200, 1.);
604 fPythia->SetPARJ(198, fQhat);
605 fPythia->SetPARJ(199, fLength);
606 fPythia->SetMSTJ(42, 2); // angular ordering
607 fPythia->SetMSTJ(44, 2); // option to run alpha_s
608 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
612 void AliGenPythia::Generate()
614 // Generate one event
615 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
616 fDecayer->ForceDecay();
618 Double_t polar[3] = {0,0,0};
619 Double_t origin[3] = {0,0,0};
621 // converts from mm/c to s
622 const Float_t kconv=0.001/2.999792458e8;
632 // Set collision vertex position
633 if (fVertexSmear == kPerEvent) Vertex();
642 // Switch hadronisation off
644 fPythia->SetMSTJ(1, 0);
648 // Quenching comes through medium-modified splitting functions.
649 AliFastGlauber::Instance()->GetRandomBHard(bimp);
650 fPythia->SetPARJ(197, bimp);
655 // Either produce new event or read partons from file
657 if (!fReadFromFile) {
663 fNpartons = fPythia->GetN();
665 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
666 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
668 LoadEvent(fRL->Stack(), 0 , 1);
673 // Run quenching routine
677 } else if (fQuench == 2){
678 fPythia->Pyquen(208., 0, 0.);
679 } else if (fQuench == 3) {
680 // Quenching is via multiplicative correction of the splittings
684 // Switch hadronisation on
686 if (fHadronisation) {
687 fPythia->SetMSTJ(1, 1);
689 // .. and perform hadronisation
690 // printf("Calling hadronisation %d\n", fPythia->GetN());
692 if (fPatchOmegaDalitz) {
693 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
695 fPythia->DalitzDecays();
696 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
701 fPythia->ImportParticles(&fParticles,"All");
703 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
704 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
712 Int_t np = fParticles.GetEntriesFast();
714 if (np == 0) continue;
718 Int_t* pParent = new Int_t[np];
719 Int_t* pSelected = new Int_t[np];
720 Int_t* trackIt = new Int_t[np];
721 for (i = 0; i < np; i++) {
727 Int_t nc = 0; // Total n. of selected particles
728 Int_t nParents = 0; // Selected parents
729 Int_t nTkbles = 0; // Trackable particles
730 if (fProcess != kPyMbDefault &&
732 fProcess != kPyMbAtlasTuneMC09 &&
733 fProcess != kPyMbWithDirectPhoton &&
734 fProcess != kPyJets &&
735 fProcess != kPyDirectGamma &&
736 fProcess != kPyMbNonDiffr &&
737 fProcess != kPyMbMSEL1 &&
740 fProcess != kPyCharmppMNRwmi &&
741 fProcess != kPyBeautyppMNRwmi &&
742 fProcess != kPyBeautyJets &&
743 fProcess != kPyWPWHG && /// !!!!!!!!!!!!!!!!! Change done for W with POHWEG !!!!!!!!!!!!!!!!!!! :)
744 fProcess != kPyJetsPWHG &&
745 fProcess != kPyCharmPWHG &&
746 fProcess != kPyBeautyPWHG) {
748 for (i = 0; i < np; i++) {
749 TParticle* iparticle = (TParticle *) fParticles.At(i);
750 Int_t ks = iparticle->GetStatusCode();
751 kf = CheckPDGCode(iparticle->GetPdgCode());
752 // No initial state partons
753 if (ks==21) continue;
755 // Heavy Flavor Selection
762 if (kfl > 100000) kfl %= 100000;
763 if (kfl > 10000) kfl %= 10000;
765 if (kfl > 10) kfl/=100;
767 if (kfl > 10) kfl/=10;
768 Int_t ipa = iparticle->GetFirstMother()-1;
771 // Establish mother daughter relation between heavy quarks and mesons
773 if (kf >= fFlavorSelect && kf <= 6) {
774 Int_t idau = iparticle->GetFirstDaughter() - 1;
776 TParticle* daughter = (TParticle *) fParticles.At(idau);
777 Int_t pdgD = daughter->GetPdgCode();
778 if (pdgD == 91 || pdgD == 92) {
779 Int_t jmin = daughter->GetFirstDaughter() - 1;
780 Int_t jmax = daughter->GetLastDaughter() - 1;
781 for (Int_t jp = jmin; jp <= jmax; jp++)
782 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
783 } // is string or cluster
789 TParticle * mother = (TParticle *) fParticles.At(ipa);
790 kfMo = TMath::Abs(mother->GetPdgCode());
793 // What to keep in Stack?
794 Bool_t flavorOK = kFALSE;
795 Bool_t selectOK = kFALSE;
797 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
799 if (kfl > fFlavorSelect) {
803 if (kfl == fFlavorSelect) flavorOK = kTRUE;
805 switch (fStackFillOpt) {
807 case kFlavorSelection:
810 case kParentSelection:
811 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
814 if (flavorOK && selectOK) {
816 // Heavy flavor hadron or quark
818 // Kinematic seletion on final state heavy flavor mesons
819 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
824 if (ParentSelected(kf)) ++nParents; // Update parent count
825 // printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
827 // Kinematic seletion on decay products
828 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
829 && !KinematicSelection(iparticle, 1))
835 // Select if mother was selected and is not tracked
837 if (pSelected[ipa] &&
838 !trackIt[ipa] && // mother will be tracked ?
839 kfMo != 5 && // mother is b-quark, don't store fragments
840 kfMo != 4 && // mother is c-quark, don't store fragments
841 kf != 92) // don't store string
844 // Semi-stable or de-selected: diselect decay products:
847 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
849 Int_t ipF = iparticle->GetFirstDaughter();
850 Int_t ipL = iparticle->GetLastDaughter();
851 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
853 // printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
854 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
857 if (pSelected[i] == -1) pSelected[i] = 0;
858 if (!pSelected[i]) continue;
859 // Count quarks only if you did not include fragmentation
860 if (fFragmentation && kf <= 10) continue;
863 // Decision on tracking
866 // Track final state particle
867 if (ks == 1) trackIt[i] = 1;
868 // Track semi-stable particles
869 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
870 // Track particles selected by process if undecayed.
871 if (fForceDecay == kNoDecay) {
872 if (ParentSelected(kf)) trackIt[i] = 1;
874 if (ParentSelected(kf)) trackIt[i] = 0;
876 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
880 } // particle selection loop
882 for (i = 0; i<np; i++) {
883 if (!pSelected[i]) continue;
884 TParticle * iparticle = (TParticle *) fParticles.At(i);
885 kf = CheckPDGCode(iparticle->GetPdgCode());
886 Int_t ks = iparticle->GetStatusCode();
887 p[0] = iparticle->Px();
888 p[1] = iparticle->Py();
889 p[2] = iparticle->Pz();
890 p[3] = iparticle->Energy();
892 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
893 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
894 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
896 Float_t tof = fTime + kconv*iparticle->T();
897 Int_t ipa = iparticle->GetFirstMother()-1;
898 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
900 PushTrack(fTrackIt*trackIt[i], iparent, kf,
901 p[0], p[1], p[2], p[3],
902 origin[0], origin[1], origin[2], tof,
903 polar[0], polar[1], polar[2],
904 kPPrimary, nt, 1., ks);
921 switch (fCountMode) {
923 // printf(" Count all \n");
927 // printf(" Count parents \n");
930 case kCountTrackables:
931 // printf(" Count trackable \n");
935 if (jev >= fNpart || fNpart == -1) {
936 fKineBias=Float_t(fNpart)/Float_t(fTrials);
938 fQ += fPythia->GetVINT(51);
939 fX1 += fPythia->GetVINT(41);
940 fX2 += fPythia->GetVINT(42);
941 fTrialsRun += fTrials;
948 SetHighWaterMark(nt);
949 // adjust weight due to kinematic selection
952 fXsection=fPythia->GetPARI(1);
955 Int_t AliGenPythia::GenerateMB()
958 // Min Bias selection and other global selections
960 Int_t i, kf, nt, iparent;
963 Double_t polar[3] = {0,0,0};
964 Double_t origin[3] = {0,0,0};
965 // converts from mm/c to s
966 const Float_t kconv=0.001/2.999792458e8;
969 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
973 Int_t* pParent = new Int_t[np];
974 for (i=0; i< np; i++) pParent[i] = -1;
975 if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
977 TParticle* jet1 = (TParticle *) fParticles.At(6);
978 TParticle* jet2 = (TParticle *) fParticles.At(7);
980 if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
987 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
988 // implemented primaryly for kPyJets, but extended to any kind of process.
989 if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
990 (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
991 Bool_t ok = TriggerOnSelectedParticles(np);
999 // Check for diffraction
1001 if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
1002 if(!CheckDiffraction()) {
1009 // Check for minimum multiplicity
1010 if (fTriggerMultiplicity > 0) {
1011 Int_t multiplicity = 0;
1012 for (i = 0; i < np; i++) {
1013 TParticle * iparticle = (TParticle *) fParticles.At(i);
1015 Int_t statusCode = iparticle->GetStatusCode();
1017 // Initial state particle
1018 if (statusCode != 1)
1021 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1024 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1027 TParticlePDG* pdgPart = iparticle->GetPDG();
1028 if (pdgPart && pdgPart->Charge() == 0)
1034 if (multiplicity < fTriggerMultiplicity) {
1038 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
1042 if (fTriggerParticle) {
1043 Bool_t triggered = kFALSE;
1044 for (i = 0; i < np; i++) {
1045 TParticle * iparticle = (TParticle *) fParticles.At(i);
1046 kf = CheckPDGCode(iparticle->GetPdgCode());
1047 if (kf != fTriggerParticle) continue;
1048 if (iparticle->Pt() == 0.) continue;
1049 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1050 if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
1061 // Check if there is a ccbar or bbbar pair with at least one of the two
1062 // in fYMin < y < fYMax
1064 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1065 TParticle *partCheck;
1067 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1068 Bool_t theChild=kFALSE;
1069 Bool_t triggered=kFALSE;
1071 Int_t pdg,mpdg,mpdgUpperFamily;
1072 for(i=0; i<np; i++) {
1073 partCheck = (TParticle*)fParticles.At(i);
1074 pdg = partCheck->GetPdgCode();
1075 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1076 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1077 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1078 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1079 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1081 if(fTriggerParticle) {
1082 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1084 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1085 Int_t mi = partCheck->GetFirstMother() - 1;
1087 mother = (TParticle*)fParticles.At(mi);
1088 mpdg=TMath::Abs(mother->GetPdgCode());
1089 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1090 if ( ParentSelected(mpdg) ||
1091 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1092 if (KinematicSelection(partCheck,1)) {
1098 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1102 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1106 if(fTriggerParticle && !triggered) { // particle requested is not produced
1113 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1115 fProcess == kPyWPWHG || /// !!!!!!!!!!!!!!!! Added for W with POWHEG !!!!!!!!!! :)
1118 fProcess == kPyMbDefault ||
1119 fProcess == kPyMb ||
1120 fProcess == kPyMbAtlasTuneMC09 ||
1121 fProcess == kPyMbWithDirectPhoton ||
1122 fProcess == kPyMbNonDiffr)
1123 && (fCutOnChild == 1) ) {
1124 if ( !CheckKinematicsOnChild() ) {
1131 for (i = 0; i < np; i++) {
1133 TParticle * iparticle = (TParticle *) fParticles.At(i);
1134 kf = CheckPDGCode(iparticle->GetPdgCode());
1135 Int_t ks = iparticle->GetStatusCode();
1136 Int_t km = iparticle->GetFirstMother();
1138 (((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
1139 ((fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) && ks == 21 && km == 0 && i>1)
1143 if (ks == 1) trackIt = 1;
1144 Int_t ipa = iparticle->GetFirstMother()-1;
1146 iparent = (ipa > -1) ? pParent[ipa] : -1;
1149 // store track information
1150 p[0] = iparticle->Px();
1151 p[1] = iparticle->Py();
1152 p[2] = iparticle->Pz();
1153 p[3] = iparticle->Energy();
1156 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1157 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1158 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1160 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
1162 PushTrack(fTrackIt*trackIt, iparent, kf,
1163 p[0], p[1], p[2], p[3],
1164 origin[0], origin[1], origin[2], tof,
1165 polar[0], polar[1], polar[2],
1166 kPPrimary, nt, 1., ks);
1170 SetHighWaterMark(nt);
1172 } // select particle
1181 void AliGenPythia::FinishRun()
1183 // Print x-section summary
1192 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1193 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1196 void AliGenPythia::AdjustWeights() const
1198 // Adjust the weights after generation of all events
1202 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1203 for (Int_t i=0; i<ntrack; i++) {
1204 part= gAlice->GetMCApp()->Particle(i);
1205 part->SetWeight(part->GetWeight()*fKineBias);
1210 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1212 // Treat protons as inside nuclei with mass numbers a1 and a2
1216 fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
1217 fUseNuclearPDF = kTRUE;
1222 void AliGenPythia::MakeHeader()
1225 // Make header for the simulated event
1228 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1229 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1232 // Builds the event header, to be called after each event
1233 if (fHeader) delete fHeader;
1234 fHeader = new AliGenPythiaEventHeader("Pythia");
1238 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1239 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1240 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1243 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1246 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1249 fHeader->SetPrimaryVertex(fVertex);
1250 fHeader->SetInteractionTime(fTime+fEventTime);
1252 // Number of primaries
1253 fHeader->SetNProduced(fNprimaries);
1256 fHeader->SetEventWeight(fPythia->GetVINT(97));
1258 // Jets that have triggered
1260 //Need to store jets for b-jet studies too!
1261 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
1264 Float_t jets[4][10];
1265 GetJets(njet, ntrig, jets);
1268 for (Int_t i = 0; i < ntrig; i++) {
1269 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1274 // Copy relevant information from external header, if present.
1279 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1280 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1282 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1285 exHeader->TriggerJet(i, uqJet);
1286 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1290 // Store quenching parameters
1293 Double_t z[4] = {0.};
1299 fPythia->GetQuenchingParameters(xp, yp, z);
1300 } else if (fQuench == 2){
1302 Double_t r1 = PARIMP.rb1;
1303 Double_t r2 = PARIMP.rb2;
1304 Double_t b = PARIMP.b1;
1305 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1306 Double_t phi = PARIMP.psib1;
1307 xp = r * TMath::Cos(phi);
1308 yp = r * TMath::Sin(phi);
1310 } else if (fQuench == 4) {
1314 AliFastGlauber::Instance()->GetSavedXY(xy);
1315 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1318 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1321 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1322 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1325 // Store pt^hard and cross-section
1326 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1327 ((AliGenPythiaEventHeader*) fHeader)->SetXsection(fPythia->GetPARI(1));
1330 // Store Event Weight
1331 ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7));
1340 Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1342 // Check the kinematic trigger condition
1345 eta[0] = jet1->Eta();
1346 eta[1] = jet2->Eta();
1348 phi[0] = jet1->Phi();
1349 phi[1] = jet2->Phi();
1351 pdg[0] = jet1->GetPdgCode();
1352 pdg[1] = jet2->GetPdgCode();
1353 Bool_t triggered = kFALSE;
1355 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) {
1358 Float_t jets[4][10];
1360 // Use Pythia clustering on parton level to determine jet axis
1362 GetJets(njets, ntrig, jets);
1364 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1369 if (pdg[0] == kGamma) {
1373 //Check eta range first...
1374 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1375 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1377 //Eta is okay, now check phi range
1378 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1379 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1390 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1392 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1394 Bool_t checking = kFALSE;
1395 Int_t j, kcode, ks, km;
1396 Int_t nPartAcc = 0; //number of particles in the acceptance range
1397 Int_t numberOfAcceptedParticles = 1;
1398 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1399 Int_t npart = fParticles.GetEntriesFast();
1401 for (j = 0; j<npart; j++) {
1402 TParticle * jparticle = (TParticle *) fParticles.At(j);
1403 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1404 ks = jparticle->GetStatusCode();
1405 km = jparticle->GetFirstMother();
1407 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1410 if( numberOfAcceptedParticles <= nPartAcc){
1419 void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1422 // Load event into Pythia Common Block
1425 Int_t npart = stack -> GetNprimary();
1429 (fPythia->GetPyjets())->N = npart;
1431 n0 = (fPythia->GetPyjets())->N;
1432 (fPythia->GetPyjets())->N = n0 + npart;
1436 for (Int_t part = 0; part < npart; part++) {
1437 TParticle *mPart = stack->Particle(part);
1439 Int_t kf = mPart->GetPdgCode();
1440 Int_t ks = mPart->GetStatusCode();
1441 Int_t idf = mPart->GetFirstDaughter();
1442 Int_t idl = mPart->GetLastDaughter();
1445 if (ks == 11 || ks == 12) {
1452 Float_t px = mPart->Px();
1453 Float_t py = mPart->Py();
1454 Float_t pz = mPart->Pz();
1455 Float_t e = mPart->Energy();
1456 Float_t m = mPart->GetCalcMass();
1459 (fPythia->GetPyjets())->P[0][part+n0] = px;
1460 (fPythia->GetPyjets())->P[1][part+n0] = py;
1461 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1462 (fPythia->GetPyjets())->P[3][part+n0] = e;
1463 (fPythia->GetPyjets())->P[4][part+n0] = m;
1465 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1466 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1467 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1468 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1469 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1473 void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1476 // Load event into Pythia Common Block
1479 Int_t npart = stack -> GetEntries();
1483 (fPythia->GetPyjets())->N = npart;
1485 n0 = (fPythia->GetPyjets())->N;
1486 (fPythia->GetPyjets())->N = n0 + npart;
1490 for (Int_t part = 0; part < npart; part++) {
1491 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1492 if (!mPart) continue;
1494 Int_t kf = mPart->GetPdgCode();
1495 Int_t ks = mPart->GetStatusCode();
1496 Int_t idf = mPart->GetFirstDaughter();
1497 Int_t idl = mPart->GetLastDaughter();
1500 if (ks == 11 || ks == 12) {
1507 Float_t px = mPart->Px();
1508 Float_t py = mPart->Py();
1509 Float_t pz = mPart->Pz();
1510 Float_t e = mPart->Energy();
1511 Float_t m = mPart->GetCalcMass();
1514 (fPythia->GetPyjets())->P[0][part+n0] = px;
1515 (fPythia->GetPyjets())->P[1][part+n0] = py;
1516 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1517 (fPythia->GetPyjets())->P[3][part+n0] = e;
1518 (fPythia->GetPyjets())->P[4][part+n0] = m;
1520 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1521 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1522 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1523 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1524 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1529 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1532 // Calls the Pythia jet finding algorithm to find jets in the current event
1537 Int_t n = fPythia->GetN();
1541 fPythia->Pycell(njets);
1543 for (i = 0; i < njets; i++) {
1544 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1545 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1546 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1547 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1558 void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1561 // Calls the Pythia clustering algorithm to find jets in the current event
1563 Int_t n = fPythia->GetN();
1566 if (fJetReconstruction == kCluster) {
1568 // Configure cluster algorithm
1570 fPythia->SetPARU(43, 2.);
1571 fPythia->SetMSTU(41, 1);
1573 // Call cluster algorithm
1575 fPythia->Pyclus(nJets);
1577 // Loading jets from common block
1583 fPythia->Pycell(nJets);
1587 for (i = 0; i < nJets; i++) {
1588 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1589 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1590 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1591 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1592 Float_t pt = TMath::Sqrt(px * px + py * py);
1593 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
1594 Float_t theta = TMath::ATan2(pt,pz);
1595 Float_t et = e * TMath::Sin(theta);
1596 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
1598 eta > fEtaMinJet && eta < fEtaMaxJet &&
1599 phi > fPhiMinJet && phi < fPhiMaxJet &&
1600 et > fEtMinJet && et < fEtMaxJet
1603 jets[0][nJetsTrig] = px;
1604 jets[1][nJetsTrig] = py;
1605 jets[2][nJetsTrig] = pz;
1606 jets[3][nJetsTrig] = e;
1608 // printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1610 // printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1615 void AliGenPythia::GetSubEventTime()
1617 // Calculates time of the next subevent
1620 TArrayF &array = *fEventsTime;
1621 fEventTime = array[fCurSubEvent++];
1623 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1627 Bool_t AliGenPythia::IsInBarrel(Float_t eta) const
1629 // Is particle in Central Barrel acceptance?
1631 if( eta < fTriggerEta )
1637 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
1639 // Is particle in EMCAL acceptance?
1640 // phi in degrees, etamin=-etamax
1641 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1648 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
1650 // Is particle in PHOS acceptance?
1651 // Acceptance slightly larger considered.
1652 // phi in degrees, etamin=-etamax
1653 // iparticle is the index of the particle to be checked, for PHOS rotation case
1655 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1660 if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1666 void AliGenPythia::RotatePhi(Bool_t& okdd)
1668 //Rotate event in phi to enhance events in PHOS acceptance
1670 if(fPHOSRotateCandidate < 0) return ;
1672 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1673 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1674 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1675 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
1677 //calculate deltaphi
1678 TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
1679 Double_t phphi = ph->Phi();
1680 Double_t deltaphi = phiPHOS - phphi;
1684 //loop for all particles and produce the phi rotation
1685 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1686 Double_t oldphi, newphi;
1687 Double_t newVx, newVy, r, vZ, time;
1688 Double_t newPx, newPy, pt, pz, e;
1689 for(Int_t i=0; i< np; i++) {
1690 TParticle* iparticle = (TParticle *) fParticles.At(i);
1691 oldphi = iparticle->Phi();
1692 newphi = oldphi + deltaphi;
1693 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1694 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1697 newVx = r * TMath::Cos(newphi);
1698 newVy = r * TMath::Sin(newphi);
1699 vZ = iparticle->Vz(); // don't transform
1700 time = iparticle->T(); // don't transform
1702 pt = iparticle->Pt();
1703 newPx = pt * TMath::Cos(newphi);
1704 newPy = pt * TMath::Sin(newphi);
1705 pz = iparticle->Pz(); // don't transform
1706 e = iparticle->Energy(); // don't transform
1709 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1710 iparticle->SetMomentum(newPx, newPy, pz, e);
1712 } //end particle loop
1714 // now let's check that we put correctly the candidate photon in PHOS
1715 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1716 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1717 if(IsInPHOS(phi,eta,-1))
1720 // reset the value for next event
1721 fPHOSRotateCandidate = -1;
1726 Bool_t AliGenPythia::CheckDiffraction()
1728 // use this method only with Perugia-0 tune!
1732 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1738 Double_t y2 = -1e10;
1740 const Int_t kNstable=20;
1741 const Int_t pdgStable[20] = {
1744 12, // Electron Neutrino
1746 14, // Muon Neutrino
1757 3112, // Sigma Minus
1764 for (Int_t i = 0; i < np; i++) {
1765 TParticle * part = (TParticle *) fParticles.At(i);
1767 Int_t statusCode = part->GetStatusCode();
1769 // Initial state particle
1770 if (statusCode != 1)
1773 Int_t pdg = TMath::Abs(part->GetPdgCode());
1774 Bool_t isStable = kFALSE;
1775 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1776 if (pdg == pdgStable[i1]) {
1784 Double_t y = part->Y();
1798 if(iPart1<0 || iPart2<0) return kFALSE;
1803 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1804 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1806 Int_t pdg1 = part1->GetPdgCode();
1807 Int_t pdg2 = part2->GetPdgCode();
1811 if (pdg1 == 2212 && pdg2 == 2212)
1819 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1822 else if (pdg1 == 2212)
1824 else if (pdg2 == 2212)
1833 TParticle * part = (TParticle *) fParticles.At(iPart);
1834 Double_t E= part->Energy();
1835 Double_t P= part->P();
1836 Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
1837 if(M2<0) return kFALSE;
1841 Double_t Mmin, Mmax, wSD, wDD, wND;
1842 if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1844 if(M>-1 && M<Mmin) return kFALSE;
1847 Int_t procType=fPythia->GetMSTI(1);
1849 if(procType== 94) proc0=1;
1850 if(procType== 92 || procType== 93) proc0=0;
1854 else if(proc0==1) proc=1;
1856 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1857 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
1860 // if(proc==1 || proc==2) return kFALSE;
1863 if(proc0!=0) fProcDiff = procType;
1864 else fProcDiff = 95;
1868 if(wSD<0) AliError("wSD<0 ! \n");
1870 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1872 // printf("iPart = %d\n", iPart);
1874 if(iPart==iPart1) fProcDiff=93;
1875 else if(iPart==iPart2) fProcDiff=92;
1877 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1886 Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1887 Double_t &wSD, Double_t &wDD, Double_t &wND)
1891 if(TMath::Abs(fEnergyCMS-900)<1 ){
1893 const Int_t nbin=400;
1895 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
1896 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
1897 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
1898 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
1899 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
1900 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
1901 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
1902 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
1903 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
1904 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
1905 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
1906 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
1907 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
1908 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
1909 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
1910 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
1911 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
1912 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
1913 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
1914 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
1915 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
1916 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
1917 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
1918 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
1919 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
1920 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
1921 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
1922 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
1923 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
1924 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
1925 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
1926 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
1927 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
1928 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1929 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1930 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1931 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1932 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1933 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1934 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1935 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1936 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1937 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
1938 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
1939 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
1940 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
1941 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
1942 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
1943 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
1944 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
1945 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
1946 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
1947 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
1948 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
1949 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
1950 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
1951 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
1952 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
1953 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
1954 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
1955 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
1956 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
1957 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
1958 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
1959 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
1960 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
1961 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1963 1.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
1964 0.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
1965 0.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
1966 0.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
1967 0.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
1968 0.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
1969 0.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
1970 0.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
1971 0.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
1972 0.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
1973 0.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
1974 0.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
1975 0.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
1976 0.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
1977 0.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
1978 0.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
1979 0.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
1980 0.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
1981 0.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
1982 0.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
1983 0.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
1984 0.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
1985 0.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
1986 0.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
1987 0.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
1988 0.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
1989 0.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
1990 0.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
1991 0.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
1992 0.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
1993 0.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
1994 0.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
1995 0.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
1996 0.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
1997 0.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
1998 0.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
1999 0.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
2000 0.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
2001 0.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
2002 0.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
2003 0.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
2004 0.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
2005 0.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
2006 0.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
2007 0.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
2008 0.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
2009 0.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
2010 0.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
2011 0.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
2012 0.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
2013 0.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
2014 0.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
2015 0.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
2016 0.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
2017 0.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
2018 0.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
2019 0.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
2020 0.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
2021 0.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
2022 0.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
2023 0.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
2024 0.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
2025 0.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
2026 0.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
2027 0.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
2028 0.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
2029 0.386112, 0.364314, 0.434375, 0.334629};
2036 if(M<Mmin || M>Mmax) return kTRUE;
2039 for(Int_t i=1; i<=nbin; i++)
2042 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2048 else if(TMath::Abs(fEnergyCMS-2760)<1 ){
2050 const Int_t nbin=400;
2052 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2053 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2054 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2055 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2056 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2057 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2058 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2059 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2060 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2061 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2062 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2063 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2064 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2065 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2066 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2067 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2068 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2069 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2070 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2071 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2072 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2073 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2074 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2075 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2076 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2077 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2078 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2079 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2080 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2081 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2082 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2083 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2084 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2085 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2086 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2087 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2088 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2089 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2090 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2091 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2092 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2093 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2094 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2095 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2096 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2097 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2098 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2099 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2100 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2101 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2102 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2103 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2104 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2105 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2106 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2107 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2108 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2109 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2110 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2111 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2112 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2113 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2114 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2115 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2116 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2117 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2118 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2120 1.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
2121 0.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
2122 0.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
2123 0.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
2124 0.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
2125 0.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
2126 0.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
2127 0.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
2128 0.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
2129 0.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
2130 0.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
2131 0.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
2132 0.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
2133 0.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
2134 0.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
2135 0.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
2136 0.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
2137 0.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
2138 0.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
2139 0.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
2140 0.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
2141 0.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
2142 0.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
2143 0.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
2144 0.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
2145 0.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
2146 0.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
2147 0.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
2148 0.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
2149 0.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
2150 0.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
2151 0.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
2152 0.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
2153 0.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
2154 0.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
2155 0.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
2156 0.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
2157 0.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
2158 0.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
2159 0.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
2160 0.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
2161 0.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
2162 0.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
2163 0.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
2164 0.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
2165 0.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
2166 0.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
2167 0.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
2168 0.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
2169 0.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
2170 0.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
2171 0.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
2172 0.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
2173 0.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
2174 0.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
2175 0.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
2176 0.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
2177 0.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
2178 0.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
2179 0.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
2180 0.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
2181 0.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
2182 0.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
2183 0.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
2184 0.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
2185 0.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
2186 0.201712, 0.242225, 0.255565, 0.258738};
2193 if(M<Mmin || M>Mmax) return kTRUE;
2196 for(Int_t i=1; i<=nbin; i++)
2199 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2207 else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2208 const Int_t nbin=400;
2210 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2211 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2212 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2213 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2214 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2215 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2216 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2217 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2218 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2219 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2220 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2221 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2222 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2223 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2224 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2225 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2226 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2227 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2228 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2229 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2230 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2231 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2232 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2233 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2234 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2235 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2236 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2237 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2238 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2239 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2240 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2241 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2242 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2243 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2244 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2245 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2246 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2247 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2248 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2249 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2250 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2251 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2252 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2253 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2254 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2255 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2256 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2257 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2258 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2259 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2260 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2261 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2262 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2263 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2264 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2265 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2266 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2267 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2268 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2269 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2270 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2271 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2272 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2273 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2274 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2275 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2276 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2278 1.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
2279 0.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
2280 0.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
2281 0.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
2282 0.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
2283 0.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
2284 0.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
2285 0.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
2286 0.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
2287 0.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
2288 0.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
2289 0.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
2290 0.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
2291 0.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
2292 0.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
2293 0.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
2294 0.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
2295 0.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
2296 0.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
2297 0.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
2298 0.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
2299 0.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
2300 0.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
2301 0.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
2302 0.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
2303 0.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
2304 0.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
2305 0.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
2306 0.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
2307 0.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
2308 0.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
2309 0.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
2310 0.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
2311 0.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
2312 0.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
2313 0.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
2314 0.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
2315 0.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
2316 0.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
2317 0.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
2318 0.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
2319 0.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
2320 0.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
2321 0.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
2322 0.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
2323 0.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
2324 0.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
2325 0.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
2326 0.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
2327 0.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
2328 0.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
2329 0.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
2330 0.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
2331 0.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
2332 0.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
2333 0.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
2334 0.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
2335 0.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
2336 0.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
2337 0.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
2338 0.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
2339 0.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
2340 0.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
2341 0.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
2342 0.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
2343 0.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
2344 0.175006, 0.223482, 0.202706, 0.218108};
2352 if(M<Mmin || M>Mmax) return kTRUE;
2355 for(Int_t i=1; i<=nbin; i++)
2358 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2368 Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2370 // Check if this is a heavy flavor decay product
2371 TParticle * part = (TParticle *) fParticles.At(ipart);
2372 Int_t mpdg = TMath::Abs(part->GetPdgCode());
2373 Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2376 if (mfl >= 4 && mfl < 6) return kTRUE;
2378 Int_t imo = part->GetFirstMother()-1;
2379 TParticle* pm = part;
2381 // Heavy flavor hadron produced by generator
2383 pm = (TParticle*)fParticles.At(imo);
2384 mpdg = TMath::Abs(pm->GetPdgCode());
2385 mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2386 if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2387 imo = pm->GetFirstMother()-1;
2392 Bool_t AliGenPythia::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
2394 // check the eta/phi correspond to the detectors acceptance
2395 // iparticle is the index of the particle to be checked, for PHOS rotation case
2396 if (fCheckPHOS && IsInPHOS (phi,eta,iparticle)) return kTRUE;
2397 else if(fCheckEMCAL && IsInEMCAL (phi,eta)) return kTRUE;
2398 else if(fCheckBarrel && IsInBarrel( eta)) return kTRUE;
2402 Bool_t AliGenPythia::TriggerOnSelectedParticles(Int_t np)
2404 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2405 // implemented primaryly for kPyJets, but extended to any kind of process.
2407 //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2408 // fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
2411 for (Int_t i=0; i< np; i++) {
2413 TParticle* iparticle = (TParticle *) fParticles.At(i);
2415 Int_t pdg = iparticle->GetPdgCode();
2416 Int_t status = iparticle->GetStatusCode();
2417 Int_t imother = iparticle->GetFirstMother() - 1;
2419 TParticle* pmother = 0x0;
2420 Int_t momStatus = -1;
2423 pmother = (TParticle *) fParticles.At(imother);
2424 momStatus = pmother->GetStatusCode();
2425 momPdg = pmother->GetPdgCode();
2431 // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2434 if (fHadronInCalo && status == 1)
2436 if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
2437 // (in case neutral mesons were declared stable)
2441 else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2445 //Fragmentation photon
2446 else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2448 if(momStatus != 11) ok = kTRUE ; // No photon from hadron decay
2451 else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2453 if( momStatus == 11)
2455 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2456 // pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2457 ok = kTRUE ; // photon from hadron decay
2459 //In case only decays from pi0 or eta requested
2460 if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2461 if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2465 // Pi0 or Eta particle
2466 else if ((fPi0InCalo || fEtaInCalo))
2468 if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2470 if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2472 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2475 else if (fEtaInCalo && pdg == 221)
2477 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2484 // Check that the selected particle is in the calorimeter acceptance
2486 if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2488 //Just check if the selected particle falls in the acceptance
2489 if(!fForceNeutralMeson2PhotonDecay )
2491 //printf("\t Check acceptance! \n");
2492 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2493 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2495 if(CheckDetectorAcceptance(phi,eta,i))
2498 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));
2499 //printf("\t Accept \n");
2504 //Mesons have several decay modes, select only those decaying into 2 photons
2505 else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2507 // In case we want the pi0/eta trigger,
2508 // check the decay mode (2 photons)
2510 //printf("\t Force decay 2 gamma\n");
2512 Int_t ndaughters = iparticle->GetNDaughters();
2513 if(ndaughters != 2){
2518 TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2519 TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2525 //iparticle->Print();
2529 Int_t pdgD1 = d1->GetPdgCode();
2530 Int_t pdgD2 = d2->GetPdgCode();
2531 //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2532 //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2534 if(pdgD1 != 22 || pdgD2 != 22){
2539 //printf("\t accept decay\n");
2541 //Trigger on the meson, not on the daughter
2542 if(!fDecayPhotonInCalo){
2544 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2545 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2547 if(CheckDetectorAcceptance(phi,eta,i))
2549 //printf("\t Accept meson pdg %d\n",pdg);
2551 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));
2559 //printf("Check daughters acceptance\n");
2561 //Trigger on the meson daughters
2563 Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2564 Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
2565 if(d1->Pt() > fTriggerParticleMinPt)
2567 //printf("\t Check acceptance photon 1! \n");
2568 if(CheckDetectorAcceptance(phi,eta,i))
2570 //printf("\t Accept Photon 1\n");
2572 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));
2580 phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2581 eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2583 if(d2->Pt() > fTriggerParticleMinPt)
2585 //printf("\t Check acceptance photon 2! \n");
2586 if(CheckDetectorAcceptance(phi,eta,i))
2588 //printf("\t Accept Photon 2\n");
2590 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));
2596 } // force 2 photon daughters in pi0/eta decays
2598 } else ok = kFALSE; // check acceptance
2602 // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2603 // A particle passing all trigger conditions except phi position in PHOS, is used as reference