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():
95 fDecayer(new AliDecayerPythia()),
103 fPhiMaxJet(2.* TMath::Pi()),
104 fJetReconstruction(kCell),
108 fPhiMaxGamma(2. * TMath::Pi()),
112 fPycellThreshold(0.),
114 fPycellMinEtJet(10.),
115 fPycellMaxRadius(1.),
116 fStackFillOpt(kFlavorSelection),
118 fFragmentation(kTRUE),
127 fTriggerMultiplicity(0),
128 fTriggerMultiplicityEta(0),
129 fTriggerMultiplicityPtMin(0),
130 fCountMode(kCountAll),
134 fFragPhotonInCalo(kFALSE),
135 fHadronInCalo(kFALSE) ,
138 fPhotonInCalo(kFALSE), // not in use
139 fDecayPhotonInCalo(kFALSE),
140 fForceNeutralMeson2PhotonDecay(kFALSE),
142 fEleInEMCAL(kFALSE), // not in use
143 fCheckBarrel(kFALSE),
146 fCheckPHOSeta(kFALSE),
147 fPHOSRotateCandidate(-1),
148 fTriggerParticleMinPt(0),
149 fPhotonMinPt(0), // not in use
150 fElectronMinPt(0), // not in use
160 // Default Constructor
162 if (!AliPythiaRndm::GetPythiaRandom())
163 AliPythiaRndm::SetPythiaRandom(GetRandom());
166 AliGenPythia::AliGenPythia(Int_t npart)
178 fInteractionRate(0.),
192 fHadronisation(kTRUE),
193 fPatchOmegaDalitz(0),
195 fReadFromFile(kFALSE),
207 fDecayer(new AliDecayerPythia()),
208 fDebugEventFirst(-1),
215 fPhiMaxJet(2.* TMath::Pi()),
216 fJetReconstruction(kCell),
220 fPhiMaxGamma(2. * TMath::Pi()),
224 fPycellThreshold(0.),
226 fPycellMinEtJet(10.),
227 fPycellMaxRadius(1.),
228 fStackFillOpt(kFlavorSelection),
230 fFragmentation(kTRUE),
239 fTriggerMultiplicity(0),
240 fTriggerMultiplicityEta(0),
241 fTriggerMultiplicityPtMin(0),
242 fCountMode(kCountAll),
246 fFragPhotonInCalo(kFALSE),
247 fHadronInCalo(kFALSE) ,
250 fPhotonInCalo(kFALSE), // not in use
251 fDecayPhotonInCalo(kFALSE),
252 fForceNeutralMeson2PhotonDecay(kFALSE),
254 fEleInEMCAL(kFALSE), // not in use
255 fCheckBarrel(kFALSE),
258 fCheckPHOSeta(kFALSE),
259 fPHOSRotateCandidate(-1),
260 fTriggerParticleMinPt(0),
261 fPhotonMinPt(0), // not in use
262 fElectronMinPt(0), // not in use
272 // default charm production at 5. 5 TeV
274 // structure function GRVHO
278 fTitle= "Particle Generator using PYTHIA";
280 // Set random number generator
281 if (!AliPythiaRndm::GetPythiaRandom())
282 AliPythiaRndm::SetPythiaRandom(GetRandom());
285 AliGenPythia::~AliGenPythia()
288 if(fEventsTime) delete fEventsTime;
291 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
293 // Generate pileup using user specified rate
294 fInteractionRate = rate;
295 fTimeWindow = timewindow;
299 void AliGenPythia::GeneratePileup()
301 // Generate sub events time for pileup
303 if(fInteractionRate == 0.) {
304 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
308 Int_t npart = NumberParticles();
310 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
314 if(fEventsTime) delete fEventsTime;
315 fEventsTime = new TArrayF(npart);
316 TArrayF &array = *fEventsTime;
317 for(Int_t ipart = 0; ipart < npart; ipart++)
320 Float_t eventtime = 0.;
323 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
324 if(eventtime > fTimeWindow) break;
325 array.Set(array.GetSize()+1);
326 array[array.GetSize()-1] = eventtime;
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;
338 SetNumberParticles(fEventsTime->GetSize());
341 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
342 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
344 // Set pycell parameters
345 fPycellEtaMax = etamax;
348 fPycellThreshold = thresh;
349 fPycellEtSeed = etseed;
350 fPycellMinEtJet = minet;
351 fPycellMaxRadius = r;
356 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
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;
365 void AliGenPythia::Init()
369 SetMC(AliPythia::Instance());
370 fPythia=(AliPythia*) fMCEvGen;
373 fParentWeight=1./Float_t(fNpart);
377 fPythia->SetCKIN(3,fPtHardMin);
378 fPythia->SetCKIN(4,fPtHardMax);
379 fPythia->SetCKIN(7,fYHardMin);
380 fPythia->SetCKIN(8,fYHardMax);
382 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
384 if (fFragmentation) {
385 fPythia->SetMSTP(111,1);
387 fPythia->SetMSTP(111,0);
391 // initial state radiation
392 fPythia->SetMSTP(61,fGinit);
393 // final state radiation
394 fPythia->SetMSTP(71,fGfinal);
397 fPythia->SetMSTP(91,1);
398 fPythia->SetPARP(91,fPtKick);
399 fPythia->SetPARP(93, 4. * fPtKick);
401 fPythia->SetMSTP(91,0);
406 fRL = AliRunLoader::Open(fkFileName, "Partons");
407 fRL->LoadKinematics();
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
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();
429 // Parent and Children Selection
432 case kPyOldUEQ2ordered:
433 case kPyOldUEQ2ordered2:
437 case kPyCharmUnforced:
438 case kPyCharmPbPbMNR:
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;
454 fParentSelect[0] = 421;
457 case kPyDPlusPbPbMNR:
460 fParentSelect[0] = 411;
463 case kPyDPlusStrangePbPbMNR:
464 case kPyDPlusStrangepPbMNR:
465 case kPyDPlusStrangeppMNR:
466 fParentSelect[0] = 431;
469 case kPyLambdacppMNR:
470 fParentSelect[0] = 4122;
475 case kPyBeautyPbPbMNR:
476 case kPyBeautypPbMNR:
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;
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;
500 fParentSelect[0] = 443;
503 case kPyMbAtlasTuneMC09:
505 case kPyMbWithDirectPhoton:
518 // JetFinder for Trigger
520 // Configure detector (EMCAL like)
522 fPythia->SetPARU(51, fPycellEtaMax);
523 fPythia->SetMSTU(51, fPycellNEta);
524 fPythia->SetMSTU(52, fPycellNPhi);
526 // Configure Jet Finder
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);
534 // This counts the total number of calls to Pyevnt() per run.
549 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
552 fPythia->SetPARJ(200, 0.0);
553 fPythia->SetPARJ(199, 0.0);
554 fPythia->SetPARJ(198, 0.0);
555 fPythia->SetPARJ(197, 0.0);
558 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
561 if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
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();
576 //read and store transverse almonds corresponding to differnt
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
588 void AliGenPythia::Generate()
590 // Generate one event
591 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
592 fDecayer->ForceDecay();
594 Double_t polar[3] = {0,0,0};
595 Double_t origin[3] = {0,0,0};
597 // converts from mm/c to s
598 const Float_t kconv=0.001/2.999792458e8;
608 // Set collision vertex position
609 if (fVertexSmear == kPerEvent) Vertex();
618 // Switch hadronisation off
620 fPythia->SetMSTJ(1, 0);
624 // Quenching comes through medium-modified splitting functions.
625 AliFastGlauber::Instance()->GetRandomBHard(bimp);
626 fPythia->SetPARJ(197, bimp);
631 // Either produce new event or read partons from file
633 if (!fReadFromFile) {
639 fNpartons = fPythia->GetN();
641 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
642 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
644 LoadEvent(fRL->Stack(), 0 , 1);
649 // Run quenching routine
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
660 // Switch hadronisation on
662 if (fHadronisation) {
663 fPythia->SetMSTJ(1, 1);
665 // .. and perform hadronisation
666 // printf("Calling hadronisation %d\n", fPythia->GetN());
668 if (fPatchOmegaDalitz) {
669 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
671 fPythia->DalitzDecays();
672 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
677 fPythia->ImportParticles(&fParticles,"All");
679 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
680 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
688 Int_t np = fParticles.GetEntriesFast();
690 if (np == 0) continue;
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++) {
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 &&
708 fProcess != kPyMbAtlasTuneMC09 &&
709 fProcess != kPyMbWithDirectPhoton &&
710 fProcess != kPyJets &&
711 fProcess != kPyDirectGamma &&
712 fProcess != kPyMbNonDiffr &&
713 fProcess != kPyMbMSEL1 &&
716 fProcess != kPyCharmppMNRwmi &&
717 fProcess != kPyBeautyppMNRwmi &&
718 fProcess != kPyBeautyJets) {
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;
727 // Heavy Flavor Selection
734 if (kfl > 100000) kfl %= 100000;
735 if (kfl > 10000) kfl %= 10000;
737 if (kfl > 10) kfl/=100;
739 if (kfl > 10) kfl/=10;
740 Int_t ipa = iparticle->GetFirstMother()-1;
743 // Establish mother daughter relation between heavy quarks and mesons
745 if (kf >= fFlavorSelect && kf <= 6) {
746 Int_t idau = iparticle->GetFirstDaughter() - 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
761 TParticle * mother = (TParticle *) fParticles.At(ipa);
762 kfMo = TMath::Abs(mother->GetPdgCode());
765 // What to keep in Stack?
766 Bool_t flavorOK = kFALSE;
767 Bool_t selectOK = kFALSE;
769 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
771 if (kfl > fFlavorSelect) {
775 if (kfl == fFlavorSelect) flavorOK = kTRUE;
777 switch (fStackFillOpt) {
779 case kFlavorSelection:
782 case kParentSelection:
783 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
786 if (flavorOK && selectOK) {
788 // Heavy flavor hadron or quark
790 // Kinematic seletion on final state heavy flavor mesons
791 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
796 if (ParentSelected(kf)) ++nParents; // Update parent count
797 // printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
799 // Kinematic seletion on decay products
800 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
801 && !KinematicSelection(iparticle, 1))
807 // Select if mother was selected and is not tracked
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
816 // Semi-stable or de-selected: diselect decay products:
819 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
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;
825 // printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
826 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
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;
835 // Decision on tracking
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;
846 if (ParentSelected(kf)) trackIt[i] = 0;
848 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
852 } // particle selection loop
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();
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]
868 Float_t tof = fTime + kconv*iparticle->T();
869 Int_t ipa = iparticle->GetFirstMother()-1;
870 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
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);
893 switch (fCountMode) {
895 // printf(" Count all \n");
899 // printf(" Count parents \n");
902 case kCountTrackables:
903 // printf(" Count trackable \n");
907 if (jev >= fNpart || fNpart == -1) {
908 fKineBias=Float_t(fNpart)/Float_t(fTrials);
910 fQ += fPythia->GetVINT(51);
911 fX1 += fPythia->GetVINT(41);
912 fX2 += fPythia->GetVINT(42);
913 fTrialsRun += fTrials;
920 SetHighWaterMark(nt);
921 // adjust weight due to kinematic selection
924 fXsection=fPythia->GetPARI(1);
927 Int_t AliGenPythia::GenerateMB()
930 // Min Bias selection and other global selections
932 Int_t i, kf, nt, iparent;
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;
941 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
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)) {
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);
969 // Check for diffraction
971 if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
972 if(!CheckDiffraction()) {
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);
985 Int_t statusCode = iparticle->GetStatusCode();
987 // Initial state particle
991 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
994 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
997 TParticlePDG* pdgPart = iparticle->GetPDG();
998 if (pdgPart && pdgPart->Charge() == 0)
1004 if (multiplicity < fTriggerMultiplicity) {
1008 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
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;
1031 // Check if there is a ccbar or bbbar pair with at least one of the two
1032 // in fYMin < y < fYMax
1034 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1035 TParticle *partCheck;
1037 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1038 Bool_t theChild=kFALSE;
1039 Bool_t triggered=kFALSE;
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;
1051 if(fTriggerParticle) {
1052 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1054 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1055 Int_t mi = partCheck->GetFirstMother() - 1;
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)) {
1068 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1072 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1076 if(fTriggerParticle && !triggered) { // particle requested is not produced
1083 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1084 if ( (fProcess == kPyW ||
1086 fProcess == kPyMbDefault ||
1087 fProcess == kPyMb ||
1088 fProcess == kPyMbAtlasTuneMC09 ||
1089 fProcess == kPyMbWithDirectPhoton ||
1090 fProcess == kPyMbNonDiffr)
1091 && (fCutOnChild == 1) ) {
1092 if ( !CheckKinematicsOnChild() ) {
1099 for (i = 0; i < np; i++) {
1101 TParticle * iparticle = (TParticle *) fParticles.At(i);
1102 kf = CheckPDGCode(iparticle->GetPdgCode());
1103 Int_t ks = iparticle->GetStatusCode();
1104 Int_t km = iparticle->GetFirstMother();
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)
1111 if (ks == 1) trackIt = 1;
1112 Int_t ipa = iparticle->GetFirstMother()-1;
1114 iparent = (ipa > -1) ? pParent[ipa] : -1;
1117 // store track information
1118 p[0] = iparticle->Px();
1119 p[1] = iparticle->Py();
1120 p[2] = iparticle->Pz();
1121 p[3] = iparticle->Energy();
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]
1128 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
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);
1138 SetHighWaterMark(nt);
1140 } // select particle
1149 void AliGenPythia::FinishRun()
1151 // Print x-section summary
1160 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1161 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1164 void AliGenPythia::AdjustWeights() const
1166 // Adjust the weights after generation of all events
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);
1178 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1180 // Treat protons as inside nuclei with mass numbers a1 and a2
1184 fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
1189 void AliGenPythia::MakeHeader()
1192 // Make header for the simulated event
1195 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1196 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1199 // Builds the event header, to be called after each event
1200 if (fHeader) delete fHeader;
1201 fHeader = new AliGenPythiaEventHeader("Pythia");
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);
1210 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1213 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1216 fHeader->SetPrimaryVertex(fVertex);
1217 fHeader->SetInteractionTime(fTime+fEventTime);
1219 // Number of primaries
1220 fHeader->SetNProduced(fNprimaries);
1222 // Jets that have triggered
1224 //Need to store jets for b-jet studies too!
1225 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
1228 Float_t jets[4][10];
1229 GetJets(njet, ntrig, jets);
1232 for (Int_t i = 0; i < ntrig; i++) {
1233 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1238 // Copy relevant information from external header, if present.
1243 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1244 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1246 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1249 exHeader->TriggerJet(i, uqJet);
1250 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1254 // Store quenching parameters
1257 Double_t z[4] = {0.};
1263 fPythia->GetQuenchingParameters(xp, yp, z);
1264 } else if (fQuench == 2){
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);
1274 } else if (fQuench == 4) {
1278 AliFastGlauber::Instance()->GetSavedXY(xy);
1279 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1282 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1285 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1286 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1290 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1298 Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1300 // Check the kinematic trigger condition
1303 eta[0] = jet1->Eta();
1304 eta[1] = jet2->Eta();
1306 phi[0] = jet1->Phi();
1307 phi[1] = jet2->Phi();
1309 pdg[0] = jet1->GetPdgCode();
1310 pdg[1] = jet2->GetPdgCode();
1311 Bool_t triggered = kFALSE;
1313 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
1316 Float_t jets[4][10];
1318 // Use Pythia clustering on parton level to determine jet axis
1320 GetJets(njets, ntrig, jets);
1322 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1327 if (pdg[0] == kGamma) {
1331 //Check eta range first...
1332 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1333 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1335 //Eta is okay, now check phi range
1336 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1337 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1348 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1350 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
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();
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();
1365 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1368 if( numberOfAcceptedParticles <= nPartAcc){
1377 void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1380 // Load event into Pythia Common Block
1383 Int_t npart = stack -> GetNprimary();
1387 (fPythia->GetPyjets())->N = npart;
1389 n0 = (fPythia->GetPyjets())->N;
1390 (fPythia->GetPyjets())->N = n0 + npart;
1394 for (Int_t part = 0; part < npart; part++) {
1395 TParticle *mPart = stack->Particle(part);
1397 Int_t kf = mPart->GetPdgCode();
1398 Int_t ks = mPart->GetStatusCode();
1399 Int_t idf = mPart->GetFirstDaughter();
1400 Int_t idl = mPart->GetLastDaughter();
1403 if (ks == 11 || ks == 12) {
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();
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;
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;
1431 void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1434 // Load event into Pythia Common Block
1437 Int_t npart = stack -> GetEntries();
1441 (fPythia->GetPyjets())->N = npart;
1443 n0 = (fPythia->GetPyjets())->N;
1444 (fPythia->GetPyjets())->N = n0 + npart;
1448 for (Int_t part = 0; part < npart; part++) {
1449 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1450 if (!mPart) continue;
1452 Int_t kf = mPart->GetPdgCode();
1453 Int_t ks = mPart->GetStatusCode();
1454 Int_t idf = mPart->GetFirstDaughter();
1455 Int_t idl = mPart->GetLastDaughter();
1458 if (ks == 11 || ks == 12) {
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();
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;
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;
1487 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1490 // Calls the Pythia jet finding algorithm to find jets in the current event
1495 Int_t n = fPythia->GetN();
1499 fPythia->Pycell(njets);
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];
1516 void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1519 // Calls the Pythia clustering algorithm to find jets in the current event
1521 Int_t n = fPythia->GetN();
1524 if (fJetReconstruction == kCluster) {
1526 // Configure cluster algorithm
1528 fPythia->SetPARU(43, 2.);
1529 fPythia->SetMSTU(41, 1);
1531 // Call cluster algorithm
1533 fPythia->Pyclus(nJets);
1535 // Loading jets from common block
1541 fPythia->Pycell(nJets);
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.));
1556 eta > fEtaMinJet && eta < fEtaMaxJet &&
1557 phi > fPhiMinJet && phi < fPhiMaxJet &&
1558 et > fEtMinJet && et < fEtMaxJet
1561 jets[0][nJetsTrig] = px;
1562 jets[1][nJetsTrig] = py;
1563 jets[2][nJetsTrig] = pz;
1564 jets[3][nJetsTrig] = e;
1566 // printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1568 // printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1573 void AliGenPythia::GetSubEventTime()
1575 // Calculates time of the next subevent
1578 TArrayF &array = *fEventsTime;
1579 fEventTime = array[fCurSubEvent++];
1581 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1585 Bool_t AliGenPythia::IsInBarrel(const Float_t eta) const
1587 // Is particle in Central Barrel acceptance?
1589 if( eta < fTriggerEta )
1595 Bool_t AliGenPythia::IsInEMCAL(const Float_t phi, const Float_t eta) const
1597 // Is particle in EMCAL acceptance?
1598 // phi in degrees, etamin=-etamax
1599 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1606 Bool_t AliGenPythia::IsInPHOS(const Float_t phi, const Float_t eta, const Int_t iparticle)
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
1613 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1618 if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1624 void AliGenPythia::RotatePhi(Bool_t& okdd)
1626 //Rotate event in phi to enhance events in PHOS acceptance
1628 if(fPHOSRotateCandidate < 0) return ;
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);
1635 //calculate deltaphi
1636 TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
1637 Double_t phphi = ph->Phi();
1638 Double_t deltaphi = phiPHOS - phphi;
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
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
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
1667 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1668 iparticle->SetMomentum(newPx, newPy, pz, e);
1670 } //end particle loop
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))
1678 // reset the value for next event
1679 fPHOSRotateCandidate = -1;
1684 Bool_t AliGenPythia::CheckDiffraction()
1686 // use this method only with Perugia-0 tune!
1690 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1696 Double_t y2 = -1e10;
1698 const Int_t kNstable=20;
1699 const Int_t pdgStable[20] = {
1702 12, // Electron Neutrino
1704 14, // Muon Neutrino
1715 3112, // Sigma Minus
1722 for (Int_t i = 0; i < np; i++) {
1723 TParticle * part = (TParticle *) fParticles.At(i);
1725 Int_t statusCode = part->GetStatusCode();
1727 // Initial state particle
1728 if (statusCode != 1)
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]) {
1742 Double_t y = part->Y();
1756 if(iPart1<0 || iPart2<0) return kFALSE;
1761 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1762 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1764 Int_t pdg1 = part1->GetPdgCode();
1765 Int_t pdg2 = part2->GetPdgCode();
1769 if (pdg1 == 2212 && pdg2 == 2212)
1777 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1780 else if (pdg1 == 2212)
1782 else if (pdg2 == 2212)
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));
1797 Double_t Mmin, Mmax, wSD, wDD, wND;
1798 if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1800 if(M>-1 && M<Mmin) return kFALSE;
1803 Int_t procType=fPythia->GetMSTI(1);
1805 if(procType== 94) proc0=1;
1806 if(procType== 92 || procType== 93) proc0=0;
1810 else if(proc0==1) proc=1;
1812 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1813 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
1816 // if(proc==1 || proc==2) return kFALSE;
1819 if(proc0!=0) fProcDiff = procType;
1820 else fProcDiff = 95;
1824 if(wSD<0) AliError("wSD<0 ! \n");
1826 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1828 // printf("iPart = %d\n", iPart);
1830 if(iPart==iPart1) fProcDiff=93;
1831 else if(iPart==iPart2) fProcDiff=92;
1833 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1842 Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1843 Double_t &wSD, Double_t &wDD, Double_t &wND)
1847 if(TMath::Abs(fEnergyCMS-900)<1 ){
1849 const Int_t nbin=400;
1851 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
1852 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
1853 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
1854 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
1855 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
1856 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
1857 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
1858 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
1859 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
1860 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
1861 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
1862 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
1863 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
1864 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
1865 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
1866 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
1867 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
1868 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
1869 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
1870 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
1871 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
1872 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
1873 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
1874 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
1875 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
1876 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
1877 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
1878 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
1879 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
1880 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
1881 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
1882 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
1883 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
1884 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1885 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1886 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1887 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1888 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1889 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1890 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1891 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1892 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1893 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
1894 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
1895 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
1896 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
1897 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
1898 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
1899 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
1900 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
1901 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
1902 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
1903 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
1904 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
1905 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
1906 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
1907 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
1908 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
1909 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
1910 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
1911 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
1912 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
1913 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
1914 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
1915 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
1916 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
1917 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1919 1.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
1920 0.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
1921 0.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
1922 0.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
1923 0.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
1924 0.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
1925 0.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
1926 0.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
1927 0.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
1928 0.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
1929 0.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
1930 0.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
1931 0.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
1932 0.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
1933 0.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
1934 0.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
1935 0.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
1936 0.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
1937 0.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
1938 0.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
1939 0.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
1940 0.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
1941 0.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
1942 0.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
1943 0.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
1944 0.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
1945 0.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
1946 0.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
1947 0.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
1948 0.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
1949 0.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
1950 0.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
1951 0.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
1952 0.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
1953 0.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
1954 0.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
1955 0.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
1956 0.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
1957 0.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
1958 0.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
1959 0.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
1960 0.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
1961 0.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
1962 0.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
1963 0.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
1964 0.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
1965 0.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
1966 0.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
1967 0.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
1968 0.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
1969 0.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
1970 0.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
1971 0.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
1972 0.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
1973 0.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
1974 0.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
1975 0.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
1976 0.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
1977 0.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
1978 0.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
1979 0.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
1980 0.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
1981 0.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
1982 0.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
1983 0.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
1984 0.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
1985 0.386112, 0.364314, 0.434375, 0.334629};
1992 if(M<Mmin || M>Mmax) return kTRUE;
1995 for(Int_t i=1; i<=nbin; i++)
1998 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2004 else if(TMath::Abs(fEnergyCMS-2760)<1 ){
2006 const Int_t nbin=400;
2008 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2009 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2010 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2011 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2012 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2013 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2014 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2015 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2016 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2017 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2018 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2019 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2020 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2021 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2022 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2023 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2024 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2025 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2026 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2027 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2028 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2029 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2030 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2031 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2032 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2033 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2034 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2035 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2036 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2037 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2038 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2039 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2040 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2041 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2042 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2043 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2044 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2045 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2046 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2047 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2048 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2049 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2050 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2051 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2052 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2053 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2054 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2055 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2056 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2057 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2058 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2059 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2060 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2061 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2062 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2063 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2064 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2065 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2066 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2067 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2068 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2069 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2070 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2071 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2072 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2073 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2074 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2076 1.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
2077 0.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
2078 0.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
2079 0.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
2080 0.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
2081 0.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
2082 0.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
2083 0.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
2084 0.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
2085 0.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
2086 0.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
2087 0.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
2088 0.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
2089 0.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
2090 0.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
2091 0.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
2092 0.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
2093 0.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
2094 0.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
2095 0.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
2096 0.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
2097 0.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
2098 0.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
2099 0.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
2100 0.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
2101 0.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
2102 0.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
2103 0.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
2104 0.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
2105 0.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
2106 0.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
2107 0.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
2108 0.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
2109 0.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
2110 0.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
2111 0.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
2112 0.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
2113 0.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
2114 0.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
2115 0.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
2116 0.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
2117 0.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
2118 0.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
2119 0.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
2120 0.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
2121 0.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
2122 0.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
2123 0.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
2124 0.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
2125 0.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
2126 0.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
2127 0.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
2128 0.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
2129 0.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
2130 0.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
2131 0.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
2132 0.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
2133 0.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
2134 0.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
2135 0.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
2136 0.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
2137 0.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
2138 0.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
2139 0.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
2140 0.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
2141 0.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
2142 0.201712, 0.242225, 0.255565, 0.258738};
2149 if(M<Mmin || M>Mmax) return kTRUE;
2152 for(Int_t i=1; i<=nbin; i++)
2155 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2163 else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2164 const Int_t nbin=400;
2166 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2167 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2168 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2169 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2170 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2171 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2172 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2173 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2174 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2175 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2176 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2177 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2178 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2179 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2180 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2181 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2182 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2183 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2184 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2185 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2186 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2187 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2188 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2189 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2190 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2191 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2192 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2193 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2194 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2195 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2196 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2197 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2198 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2199 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2200 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2201 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2202 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2203 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2204 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2205 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2206 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2207 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2208 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2209 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2210 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2211 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2212 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2213 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2214 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2215 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2216 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2217 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2218 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2219 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2220 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2221 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2222 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2223 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2224 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2225 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2226 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2227 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2228 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2229 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2230 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2231 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2232 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2234 1.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
2235 0.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
2236 0.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
2237 0.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
2238 0.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
2239 0.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
2240 0.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
2241 0.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
2242 0.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
2243 0.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
2244 0.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
2245 0.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
2246 0.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
2247 0.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
2248 0.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
2249 0.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
2250 0.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
2251 0.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
2252 0.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
2253 0.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
2254 0.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
2255 0.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
2256 0.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
2257 0.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
2258 0.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
2259 0.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
2260 0.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
2261 0.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
2262 0.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
2263 0.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
2264 0.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
2265 0.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
2266 0.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
2267 0.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
2268 0.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
2269 0.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
2270 0.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
2271 0.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
2272 0.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
2273 0.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
2274 0.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
2275 0.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
2276 0.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
2277 0.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
2278 0.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
2279 0.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
2280 0.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
2281 0.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
2282 0.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
2283 0.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
2284 0.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
2285 0.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
2286 0.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
2287 0.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
2288 0.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
2289 0.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
2290 0.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
2291 0.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
2292 0.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
2293 0.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
2294 0.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
2295 0.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
2296 0.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
2297 0.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
2298 0.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
2299 0.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
2300 0.175006, 0.223482, 0.202706, 0.218108};
2308 if(M<Mmin || M>Mmax) return kTRUE;
2311 for(Int_t i=1; i<=nbin; i++)
2314 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2324 Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
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))));
2332 if (mfl >= 4 && mfl < 6) return kTRUE;
2334 Int_t imo = part->GetFirstMother()-1;
2335 TParticle* pm = part;
2337 // Heavy flavor hadron produced by generator
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;
2348 Bool_t AliGenPythia::CheckDetectorAcceptance(const Float_t phi, const Float_t eta, const Int_t iparticle)
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;
2358 Bool_t AliGenPythia::TriggerOnSelectedParticles(const Int_t np)
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.
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);
2367 for (Int_t i=0; i< np; i++) {
2369 TParticle* iparticle = (TParticle *) fParticles.At(i);
2371 Int_t pdg = iparticle->GetPdgCode();
2372 Int_t status = iparticle->GetStatusCode();
2373 Int_t imother = iparticle->GetFirstMother() - 1;
2375 TParticle* pmother = 0x0;
2376 Int_t momStatus = -1;
2379 pmother = (TParticle *) fParticles.At(imother);
2380 momStatus = pmother->GetStatusCode();
2381 momPdg = pmother->GetPdgCode();
2387 // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2390 if (fHadronInCalo && status == 1)
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)
2397 else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2401 //Fragmentation photon
2402 else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2404 if(momStatus != 11) ok = kTRUE ; // No photon from hadron decay
2407 else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2409 if( momStatus == 11)
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
2415 //In case only decays from pi0 or eta requested
2416 if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2417 if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2421 // Pi0 or Eta particle
2422 else if ((fPi0InCalo || fEtaInCalo))
2424 if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2426 if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2428 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2431 else if (fEtaInCalo && pdg == 221)
2433 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2440 // Check that the selected particle is in the calorimeter acceptance
2442 if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2444 //Just check if the selected particle falls in the acceptance
2445 if(!fForceNeutralMeson2PhotonDecay )
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
2451 if(CheckDetectorAcceptance(phi,eta,i))
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");
2460 //Mesons have several decay modes, select only those decaying into 2 photons
2461 else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2463 // In case we want the pi0/eta trigger,
2464 // check the decay mode (2 photons)
2466 //printf("\t Force decay 2 gamma\n");
2468 Int_t ndaughters = iparticle->GetNDaughters();
2469 if(ndaughters != 2){
2474 TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2475 TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2481 //iparticle->Print();
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());
2490 if(pdgD1 != 22 || pdgD2 != 22){
2495 //printf("\t accept decay\n");
2497 //Trigger on the meson, not on the daughter
2498 if(!fDecayPhotonInCalo){
2500 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2501 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2503 if(CheckDetectorAcceptance(phi,eta,i))
2505 //printf("\t Accept meson pdg %d\n",pdg);
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));
2515 //printf("Check daughters acceptance\n");
2517 //Trigger on the meson daughters
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)
2523 //printf("\t Check acceptance photon 1! \n");
2524 if(CheckDetectorAcceptance(phi,eta,i))
2526 //printf("\t Accept Photon 1\n");
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));
2536 phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2537 eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2539 if(d2->Pt() > fTriggerParticleMinPt)
2541 //printf("\t Check acceptance photon 2! \n");
2542 if(CheckDetectorAcceptance(phi,eta,i))
2544 //printf("\t Accept Photon 2\n");
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));
2552 } // force 2 photon daughters in pi0/eta decays
2554 } else ok = kFALSE; // check acceptance
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