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);
383 if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
386 if (fFragmentation) {
387 fPythia->SetMSTP(111,1);
389 fPythia->SetMSTP(111,0);
393 // initial state radiation
394 fPythia->SetMSTP(61,fGinit);
395 // final state radiation
396 fPythia->SetMSTP(71,fGfinal);
399 fPythia->SetMSTP(91,1);
400 fPythia->SetPARP(91,fPtKick);
401 fPythia->SetPARP(93, 4. * fPtKick);
403 fPythia->SetMSTP(91,0);
408 fRL = AliRunLoader::Open(fkFileName, "Partons");
409 fRL->LoadKinematics();
415 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
416 // Forward Paramters to the AliPythia object
417 fDecayer->SetForceDecay(fForceDecay);
418 // Switch off Heavy Flavors on request
420 // Maximum number of quark flavours used in pdf
421 fPythia->SetMSTP(58, 3);
422 // Maximum number of flavors that can be used in showers
423 fPythia->SetMSTJ(45, 3);
424 // Switch off g->QQbar splitting in decay table
425 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
431 // Parent and Children Selection
434 case kPyOldUEQ2ordered:
435 case kPyOldUEQ2ordered2:
439 case kPyCharmUnforced:
440 case kPyCharmPbPbMNR:
443 case kPyCharmppMNRwmi:
444 fParentSelect[0] = 411;
445 fParentSelect[1] = 421;
446 fParentSelect[2] = 431;
447 fParentSelect[3] = 4122;
448 fParentSelect[4] = 4232;
449 fParentSelect[5] = 4132;
450 fParentSelect[6] = 4332;
456 fParentSelect[0] = 421;
459 case kPyDPlusPbPbMNR:
462 fParentSelect[0] = 411;
465 case kPyDPlusStrangePbPbMNR:
466 case kPyDPlusStrangepPbMNR:
467 case kPyDPlusStrangeppMNR:
468 fParentSelect[0] = 431;
471 case kPyLambdacppMNR:
472 fParentSelect[0] = 4122;
477 case kPyBeautyPbPbMNR:
478 case kPyBeautypPbMNR:
480 case kPyBeautyppMNRwmi:
481 fParentSelect[0]= 511;
482 fParentSelect[1]= 521;
483 fParentSelect[2]= 531;
484 fParentSelect[3]= 5122;
485 fParentSelect[4]= 5132;
486 fParentSelect[5]= 5232;
487 fParentSelect[6]= 5332;
490 case kPyBeautyUnforced:
491 fParentSelect[0] = 511;
492 fParentSelect[1] = 521;
493 fParentSelect[2] = 531;
494 fParentSelect[3] = 5122;
495 fParentSelect[4] = 5132;
496 fParentSelect[5] = 5232;
497 fParentSelect[6] = 5332;
502 fParentSelect[0] = 443;
505 case kPyMbAtlasTuneMC09:
507 case kPyMbWithDirectPhoton:
516 case kPyMBRSingleDiffraction:
517 case kPyMBRDoubleDiffraction:
518 case kPyMBRCentralDiffraction:
523 // JetFinder for Trigger
525 // Configure detector (EMCAL like)
527 fPythia->SetPARU(51, fPycellEtaMax);
528 fPythia->SetMSTU(51, fPycellNEta);
529 fPythia->SetMSTU(52, fPycellNPhi);
531 // Configure Jet Finder
533 fPythia->SetPARU(58, fPycellThreshold);
534 fPythia->SetPARU(52, fPycellEtSeed);
535 fPythia->SetPARU(53, fPycellMinEtJet);
536 fPythia->SetPARU(54, fPycellMaxRadius);
537 fPythia->SetMSTU(54, 2);
539 // This counts the total number of calls to Pyevnt() per run.
554 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
557 fPythia->SetPARJ(200, 0.0);
558 fPythia->SetPARJ(199, 0.0);
559 fPythia->SetPARJ(198, 0.0);
560 fPythia->SetPARJ(197, 0.0);
563 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
566 if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
569 // Nestor's change of the splittings
570 fPythia->SetPARJ(200, 0.8);
571 fPythia->SetMSTJ(41, 1); // QCD radiation only
572 fPythia->SetMSTJ(42, 2); // angular ordering
573 fPythia->SetMSTJ(44, 2); // option to run alpha_s
574 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
575 fPythia->SetMSTJ(50, 0); // No coherence in first branching
576 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
577 } else if (fQuench == 4) {
578 // Armesto-Cunqueiro-Salgado change of the splittings.
579 AliFastGlauber* glauber = AliFastGlauber::Instance();
581 //read and store transverse almonds corresponding to differnt
583 glauber->SetCentralityClass(0.,0.1);
584 fPythia->SetPARJ(200, 1.);
585 fPythia->SetPARJ(198, fQhat);
586 fPythia->SetPARJ(199, fLength);
587 fPythia->SetMSTJ(42, 2); // angular ordering
588 fPythia->SetMSTJ(44, 2); // option to run alpha_s
589 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
593 void AliGenPythia::Generate()
595 // Generate one event
596 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
597 fDecayer->ForceDecay();
599 Double_t polar[3] = {0,0,0};
600 Double_t origin[3] = {0,0,0};
602 // converts from mm/c to s
603 const Float_t kconv=0.001/2.999792458e8;
613 // Set collision vertex position
614 if (fVertexSmear == kPerEvent) Vertex();
623 // Switch hadronisation off
625 fPythia->SetMSTJ(1, 0);
629 // Quenching comes through medium-modified splitting functions.
630 AliFastGlauber::Instance()->GetRandomBHard(bimp);
631 fPythia->SetPARJ(197, bimp);
636 // Either produce new event or read partons from file
638 if (!fReadFromFile) {
644 fNpartons = fPythia->GetN();
646 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
647 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
649 LoadEvent(fRL->Stack(), 0 , 1);
654 // Run quenching routine
658 } else if (fQuench == 2){
659 fPythia->Pyquen(208., 0, 0.);
660 } else if (fQuench == 3) {
661 // Quenching is via multiplicative correction of the splittings
665 // Switch hadronisation on
667 if (fHadronisation) {
668 fPythia->SetMSTJ(1, 1);
670 // .. and perform hadronisation
671 // printf("Calling hadronisation %d\n", fPythia->GetN());
673 if (fPatchOmegaDalitz) {
674 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
676 fPythia->DalitzDecays();
677 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
682 fPythia->ImportParticles(&fParticles,"All");
684 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
685 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
693 Int_t np = fParticles.GetEntriesFast();
695 if (np == 0) continue;
699 Int_t* pParent = new Int_t[np];
700 Int_t* pSelected = new Int_t[np];
701 Int_t* trackIt = new Int_t[np];
702 for (i = 0; i < np; i++) {
708 Int_t nc = 0; // Total n. of selected particles
709 Int_t nParents = 0; // Selected parents
710 Int_t nTkbles = 0; // Trackable particles
711 if (fProcess != kPyMbDefault &&
713 fProcess != kPyMbAtlasTuneMC09 &&
714 fProcess != kPyMbWithDirectPhoton &&
715 fProcess != kPyJets &&
716 fProcess != kPyDirectGamma &&
717 fProcess != kPyMbNonDiffr &&
718 fProcess != kPyMbMSEL1 &&
721 fProcess != kPyCharmppMNRwmi &&
722 fProcess != kPyBeautyppMNRwmi &&
723 fProcess != kPyBeautyJets) {
725 for (i = 0; i < np; i++) {
726 TParticle* iparticle = (TParticle *) fParticles.At(i);
727 Int_t ks = iparticle->GetStatusCode();
728 kf = CheckPDGCode(iparticle->GetPdgCode());
729 // No initial state partons
730 if (ks==21) continue;
732 // Heavy Flavor Selection
739 if (kfl > 100000) kfl %= 100000;
740 if (kfl > 10000) kfl %= 10000;
742 if (kfl > 10) kfl/=100;
744 if (kfl > 10) kfl/=10;
745 Int_t ipa = iparticle->GetFirstMother()-1;
748 // Establish mother daughter relation between heavy quarks and mesons
750 if (kf >= fFlavorSelect && kf <= 6) {
751 Int_t idau = iparticle->GetFirstDaughter() - 1;
753 TParticle* daughter = (TParticle *) fParticles.At(idau);
754 Int_t pdgD = daughter->GetPdgCode();
755 if (pdgD == 91 || pdgD == 92) {
756 Int_t jmin = daughter->GetFirstDaughter() - 1;
757 Int_t jmax = daughter->GetLastDaughter() - 1;
758 for (Int_t jp = jmin; jp <= jmax; jp++)
759 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
760 } // is string or cluster
766 TParticle * mother = (TParticle *) fParticles.At(ipa);
767 kfMo = TMath::Abs(mother->GetPdgCode());
770 // What to keep in Stack?
771 Bool_t flavorOK = kFALSE;
772 Bool_t selectOK = kFALSE;
774 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
776 if (kfl > fFlavorSelect) {
780 if (kfl == fFlavorSelect) flavorOK = kTRUE;
782 switch (fStackFillOpt) {
784 case kFlavorSelection:
787 case kParentSelection:
788 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
791 if (flavorOK && selectOK) {
793 // Heavy flavor hadron or quark
795 // Kinematic seletion on final state heavy flavor mesons
796 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
801 if (ParentSelected(kf)) ++nParents; // Update parent count
802 // printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
804 // Kinematic seletion on decay products
805 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
806 && !KinematicSelection(iparticle, 1))
812 // Select if mother was selected and is not tracked
814 if (pSelected[ipa] &&
815 !trackIt[ipa] && // mother will be tracked ?
816 kfMo != 5 && // mother is b-quark, don't store fragments
817 kfMo != 4 && // mother is c-quark, don't store fragments
818 kf != 92) // don't store string
821 // Semi-stable or de-selected: diselect decay products:
824 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
826 Int_t ipF = iparticle->GetFirstDaughter();
827 Int_t ipL = iparticle->GetLastDaughter();
828 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
830 // printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
831 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
834 if (pSelected[i] == -1) pSelected[i] = 0;
835 if (!pSelected[i]) continue;
836 // Count quarks only if you did not include fragmentation
837 if (fFragmentation && kf <= 10) continue;
840 // Decision on tracking
843 // Track final state particle
844 if (ks == 1) trackIt[i] = 1;
845 // Track semi-stable particles
846 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
847 // Track particles selected by process if undecayed.
848 if (fForceDecay == kNoDecay) {
849 if (ParentSelected(kf)) trackIt[i] = 1;
851 if (ParentSelected(kf)) trackIt[i] = 0;
853 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
857 } // particle selection loop
859 for (i = 0; i<np; i++) {
860 if (!pSelected[i]) continue;
861 TParticle * iparticle = (TParticle *) fParticles.At(i);
862 kf = CheckPDGCode(iparticle->GetPdgCode());
863 Int_t ks = iparticle->GetStatusCode();
864 p[0] = iparticle->Px();
865 p[1] = iparticle->Py();
866 p[2] = iparticle->Pz();
867 p[3] = iparticle->Energy();
869 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
870 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
871 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
873 Float_t tof = fTime + kconv*iparticle->T();
874 Int_t ipa = iparticle->GetFirstMother()-1;
875 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
877 PushTrack(fTrackIt*trackIt[i], iparent, kf,
878 p[0], p[1], p[2], p[3],
879 origin[0], origin[1], origin[2], tof,
880 polar[0], polar[1], polar[2],
881 kPPrimary, nt, 1., ks);
898 switch (fCountMode) {
900 // printf(" Count all \n");
904 // printf(" Count parents \n");
907 case kCountTrackables:
908 // printf(" Count trackable \n");
912 if (jev >= fNpart || fNpart == -1) {
913 fKineBias=Float_t(fNpart)/Float_t(fTrials);
915 fQ += fPythia->GetVINT(51);
916 fX1 += fPythia->GetVINT(41);
917 fX2 += fPythia->GetVINT(42);
918 fTrialsRun += fTrials;
925 SetHighWaterMark(nt);
926 // adjust weight due to kinematic selection
929 fXsection=fPythia->GetPARI(1);
932 Int_t AliGenPythia::GenerateMB()
935 // Min Bias selection and other global selections
937 Int_t i, kf, nt, iparent;
940 Double_t polar[3] = {0,0,0};
941 Double_t origin[3] = {0,0,0};
942 // converts from mm/c to s
943 const Float_t kconv=0.001/2.999792458e8;
946 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
950 Int_t* pParent = new Int_t[np];
951 for (i=0; i< np; i++) pParent[i] = -1;
952 if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
954 TParticle* jet1 = (TParticle *) fParticles.At(6);
955 TParticle* jet2 = (TParticle *) fParticles.At(7);
957 if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
964 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
965 // implemented primaryly for kPyJets, but extended to any kind of process.
966 if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
967 (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
968 Bool_t ok = TriggerOnSelectedParticles(np);
976 // Check for diffraction
978 if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
979 if(!CheckDiffraction()) {
986 // Check for minimum multiplicity
987 if (fTriggerMultiplicity > 0) {
988 Int_t multiplicity = 0;
989 for (i = 0; i < np; i++) {
990 TParticle * iparticle = (TParticle *) fParticles.At(i);
992 Int_t statusCode = iparticle->GetStatusCode();
994 // Initial state particle
998 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1001 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1004 TParticlePDG* pdgPart = iparticle->GetPDG();
1005 if (pdgPart && pdgPart->Charge() == 0)
1011 if (multiplicity < fTriggerMultiplicity) {
1015 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
1019 if (fTriggerParticle) {
1020 Bool_t triggered = kFALSE;
1021 for (i = 0; i < np; i++) {
1022 TParticle * iparticle = (TParticle *) fParticles.At(i);
1023 kf = CheckPDGCode(iparticle->GetPdgCode());
1024 if (kf != fTriggerParticle) continue;
1025 if (iparticle->Pt() == 0.) continue;
1026 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1027 if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
1038 // Check if there is a ccbar or bbbar pair with at least one of the two
1039 // in fYMin < y < fYMax
1041 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1042 TParticle *partCheck;
1044 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1045 Bool_t theChild=kFALSE;
1046 Bool_t triggered=kFALSE;
1048 Int_t pdg,mpdg,mpdgUpperFamily;
1049 for(i=0; i<np; i++) {
1050 partCheck = (TParticle*)fParticles.At(i);
1051 pdg = partCheck->GetPdgCode();
1052 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1053 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1054 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1055 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1056 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1058 if(fTriggerParticle) {
1059 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1061 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1062 Int_t mi = partCheck->GetFirstMother() - 1;
1064 mother = (TParticle*)fParticles.At(mi);
1065 mpdg=TMath::Abs(mother->GetPdgCode());
1066 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1067 if ( ParentSelected(mpdg) ||
1068 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1069 if (KinematicSelection(partCheck,1)) {
1075 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1079 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1083 if(fTriggerParticle && !triggered) { // particle requested is not produced
1090 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1091 if ( (fProcess == kPyW ||
1093 fProcess == kPyMbDefault ||
1094 fProcess == kPyMb ||
1095 fProcess == kPyMbAtlasTuneMC09 ||
1096 fProcess == kPyMbWithDirectPhoton ||
1097 fProcess == kPyMbNonDiffr)
1098 && (fCutOnChild == 1) ) {
1099 if ( !CheckKinematicsOnChild() ) {
1106 for (i = 0; i < np; i++) {
1108 TParticle * iparticle = (TParticle *) fParticles.At(i);
1109 kf = CheckPDGCode(iparticle->GetPdgCode());
1110 Int_t ks = iparticle->GetStatusCode();
1111 Int_t km = iparticle->GetFirstMother();
1113 (((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
1114 ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
1118 if (ks == 1) trackIt = 1;
1119 Int_t ipa = iparticle->GetFirstMother()-1;
1121 iparent = (ipa > -1) ? pParent[ipa] : -1;
1124 // store track information
1125 p[0] = iparticle->Px();
1126 p[1] = iparticle->Py();
1127 p[2] = iparticle->Pz();
1128 p[3] = iparticle->Energy();
1131 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1132 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1133 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1135 Float_t tof = fTime + fEventTime + kconv * iparticle->T();
1137 PushTrack(fTrackIt*trackIt, iparent, kf,
1138 p[0], p[1], p[2], p[3],
1139 origin[0], origin[1], origin[2], tof,
1140 polar[0], polar[1], polar[2],
1141 kPPrimary, nt, 1., ks);
1145 SetHighWaterMark(nt);
1147 } // select particle
1156 void AliGenPythia::FinishRun()
1158 // Print x-section summary
1167 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1168 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1171 void AliGenPythia::AdjustWeights() const
1173 // Adjust the weights after generation of all events
1177 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1178 for (Int_t i=0; i<ntrack; i++) {
1179 part= gAlice->GetMCApp()->Particle(i);
1180 part->SetWeight(part->GetWeight()*fKineBias);
1185 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1187 // Treat protons as inside nuclei with mass numbers a1 and a2
1191 fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
1196 void AliGenPythia::MakeHeader()
1199 // Make header for the simulated event
1202 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1203 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1206 // Builds the event header, to be called after each event
1207 if (fHeader) delete fHeader;
1208 fHeader = new AliGenPythiaEventHeader("Pythia");
1212 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1213 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1214 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1217 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1220 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1223 fHeader->SetPrimaryVertex(fVertex);
1224 fHeader->SetInteractionTime(fTime+fEventTime);
1226 // Number of primaries
1227 fHeader->SetNProduced(fNprimaries);
1229 // Jets that have triggered
1231 //Need to store jets for b-jet studies too!
1232 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
1235 Float_t jets[4][10];
1236 GetJets(njet, ntrig, jets);
1239 for (Int_t i = 0; i < ntrig; i++) {
1240 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1245 // Copy relevant information from external header, if present.
1250 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1251 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1253 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1256 exHeader->TriggerJet(i, uqJet);
1257 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1261 // Store quenching parameters
1264 Double_t z[4] = {0.};
1270 fPythia->GetQuenchingParameters(xp, yp, z);
1271 } else if (fQuench == 2){
1273 Double_t r1 = PARIMP.rb1;
1274 Double_t r2 = PARIMP.rb2;
1275 Double_t b = PARIMP.b1;
1276 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1277 Double_t phi = PARIMP.psib1;
1278 xp = r * TMath::Cos(phi);
1279 yp = r * TMath::Sin(phi);
1281 } else if (fQuench == 4) {
1285 AliFastGlauber::Instance()->GetSavedXY(xy);
1286 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1289 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1292 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1293 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1297 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1305 Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1307 // Check the kinematic trigger condition
1310 eta[0] = jet1->Eta();
1311 eta[1] = jet2->Eta();
1313 phi[0] = jet1->Phi();
1314 phi[1] = jet2->Phi();
1316 pdg[0] = jet1->GetPdgCode();
1317 pdg[1] = jet2->GetPdgCode();
1318 Bool_t triggered = kFALSE;
1320 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
1323 Float_t jets[4][10];
1325 // Use Pythia clustering on parton level to determine jet axis
1327 GetJets(njets, ntrig, jets);
1329 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1334 if (pdg[0] == kGamma) {
1338 //Check eta range first...
1339 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1340 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1342 //Eta is okay, now check phi range
1343 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1344 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1355 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1357 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1359 Bool_t checking = kFALSE;
1360 Int_t j, kcode, ks, km;
1361 Int_t nPartAcc = 0; //number of particles in the acceptance range
1362 Int_t numberOfAcceptedParticles = 1;
1363 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1364 Int_t npart = fParticles.GetEntriesFast();
1366 for (j = 0; j<npart; j++) {
1367 TParticle * jparticle = (TParticle *) fParticles.At(j);
1368 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1369 ks = jparticle->GetStatusCode();
1370 km = jparticle->GetFirstMother();
1372 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1375 if( numberOfAcceptedParticles <= nPartAcc){
1384 void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1387 // Load event into Pythia Common Block
1390 Int_t npart = stack -> GetNprimary();
1394 (fPythia->GetPyjets())->N = npart;
1396 n0 = (fPythia->GetPyjets())->N;
1397 (fPythia->GetPyjets())->N = n0 + npart;
1401 for (Int_t part = 0; part < npart; part++) {
1402 TParticle *mPart = stack->Particle(part);
1404 Int_t kf = mPart->GetPdgCode();
1405 Int_t ks = mPart->GetStatusCode();
1406 Int_t idf = mPart->GetFirstDaughter();
1407 Int_t idl = mPart->GetLastDaughter();
1410 if (ks == 11 || ks == 12) {
1417 Float_t px = mPart->Px();
1418 Float_t py = mPart->Py();
1419 Float_t pz = mPart->Pz();
1420 Float_t e = mPart->Energy();
1421 Float_t m = mPart->GetCalcMass();
1424 (fPythia->GetPyjets())->P[0][part+n0] = px;
1425 (fPythia->GetPyjets())->P[1][part+n0] = py;
1426 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1427 (fPythia->GetPyjets())->P[3][part+n0] = e;
1428 (fPythia->GetPyjets())->P[4][part+n0] = m;
1430 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1431 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1432 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1433 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1434 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1438 void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1441 // Load event into Pythia Common Block
1444 Int_t npart = stack -> GetEntries();
1448 (fPythia->GetPyjets())->N = npart;
1450 n0 = (fPythia->GetPyjets())->N;
1451 (fPythia->GetPyjets())->N = n0 + npart;
1455 for (Int_t part = 0; part < npart; part++) {
1456 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1457 if (!mPart) continue;
1459 Int_t kf = mPart->GetPdgCode();
1460 Int_t ks = mPart->GetStatusCode();
1461 Int_t idf = mPart->GetFirstDaughter();
1462 Int_t idl = mPart->GetLastDaughter();
1465 if (ks == 11 || ks == 12) {
1472 Float_t px = mPart->Px();
1473 Float_t py = mPart->Py();
1474 Float_t pz = mPart->Pz();
1475 Float_t e = mPart->Energy();
1476 Float_t m = mPart->GetCalcMass();
1479 (fPythia->GetPyjets())->P[0][part+n0] = px;
1480 (fPythia->GetPyjets())->P[1][part+n0] = py;
1481 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1482 (fPythia->GetPyjets())->P[3][part+n0] = e;
1483 (fPythia->GetPyjets())->P[4][part+n0] = m;
1485 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1486 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1487 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1488 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1489 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1494 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1497 // Calls the Pythia jet finding algorithm to find jets in the current event
1502 Int_t n = fPythia->GetN();
1506 fPythia->Pycell(njets);
1508 for (i = 0; i < njets; i++) {
1509 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1510 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1511 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1512 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1523 void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1526 // Calls the Pythia clustering algorithm to find jets in the current event
1528 Int_t n = fPythia->GetN();
1531 if (fJetReconstruction == kCluster) {
1533 // Configure cluster algorithm
1535 fPythia->SetPARU(43, 2.);
1536 fPythia->SetMSTU(41, 1);
1538 // Call cluster algorithm
1540 fPythia->Pyclus(nJets);
1542 // Loading jets from common block
1548 fPythia->Pycell(nJets);
1552 for (i = 0; i < nJets; i++) {
1553 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1554 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1555 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1556 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1557 Float_t pt = TMath::Sqrt(px * px + py * py);
1558 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
1559 Float_t theta = TMath::ATan2(pt,pz);
1560 Float_t et = e * TMath::Sin(theta);
1561 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
1563 eta > fEtaMinJet && eta < fEtaMaxJet &&
1564 phi > fPhiMinJet && phi < fPhiMaxJet &&
1565 et > fEtMinJet && et < fEtMaxJet
1568 jets[0][nJetsTrig] = px;
1569 jets[1][nJetsTrig] = py;
1570 jets[2][nJetsTrig] = pz;
1571 jets[3][nJetsTrig] = e;
1573 // printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1575 // printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1580 void AliGenPythia::GetSubEventTime()
1582 // Calculates time of the next subevent
1585 TArrayF &array = *fEventsTime;
1586 fEventTime = array[fCurSubEvent++];
1588 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1592 Bool_t AliGenPythia::IsInBarrel(Float_t eta) const
1594 // Is particle in Central Barrel acceptance?
1596 if( eta < fTriggerEta )
1602 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
1604 // Is particle in EMCAL acceptance?
1605 // phi in degrees, etamin=-etamax
1606 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1613 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
1615 // Is particle in PHOS acceptance?
1616 // Acceptance slightly larger considered.
1617 // phi in degrees, etamin=-etamax
1618 // iparticle is the index of the particle to be checked, for PHOS rotation case
1620 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1625 if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1631 void AliGenPythia::RotatePhi(Bool_t& okdd)
1633 //Rotate event in phi to enhance events in PHOS acceptance
1635 if(fPHOSRotateCandidate < 0) return ;
1637 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1638 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1639 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1640 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
1642 //calculate deltaphi
1643 TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
1644 Double_t phphi = ph->Phi();
1645 Double_t deltaphi = phiPHOS - phphi;
1649 //loop for all particles and produce the phi rotation
1650 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1651 Double_t oldphi, newphi;
1652 Double_t newVx, newVy, r, vZ, time;
1653 Double_t newPx, newPy, pt, pz, e;
1654 for(Int_t i=0; i< np; i++) {
1655 TParticle* iparticle = (TParticle *) fParticles.At(i);
1656 oldphi = iparticle->Phi();
1657 newphi = oldphi + deltaphi;
1658 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1659 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1662 newVx = r * TMath::Cos(newphi);
1663 newVy = r * TMath::Sin(newphi);
1664 vZ = iparticle->Vz(); // don't transform
1665 time = iparticle->T(); // don't transform
1667 pt = iparticle->Pt();
1668 newPx = pt * TMath::Cos(newphi);
1669 newPy = pt * TMath::Sin(newphi);
1670 pz = iparticle->Pz(); // don't transform
1671 e = iparticle->Energy(); // don't transform
1674 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1675 iparticle->SetMomentum(newPx, newPy, pz, e);
1677 } //end particle loop
1679 // now let's check that we put correctly the candidate photon in PHOS
1680 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1681 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1682 if(IsInPHOS(phi,eta,-1))
1685 // reset the value for next event
1686 fPHOSRotateCandidate = -1;
1691 Bool_t AliGenPythia::CheckDiffraction()
1693 // use this method only with Perugia-0 tune!
1697 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1703 Double_t y2 = -1e10;
1705 const Int_t kNstable=20;
1706 const Int_t pdgStable[20] = {
1709 12, // Electron Neutrino
1711 14, // Muon Neutrino
1722 3112, // Sigma Minus
1729 for (Int_t i = 0; i < np; i++) {
1730 TParticle * part = (TParticle *) fParticles.At(i);
1732 Int_t statusCode = part->GetStatusCode();
1734 // Initial state particle
1735 if (statusCode != 1)
1738 Int_t pdg = TMath::Abs(part->GetPdgCode());
1739 Bool_t isStable = kFALSE;
1740 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1741 if (pdg == pdgStable[i1]) {
1749 Double_t y = part->Y();
1763 if(iPart1<0 || iPart2<0) return kFALSE;
1768 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1769 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1771 Int_t pdg1 = part1->GetPdgCode();
1772 Int_t pdg2 = part2->GetPdgCode();
1776 if (pdg1 == 2212 && pdg2 == 2212)
1784 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1787 else if (pdg1 == 2212)
1789 else if (pdg2 == 2212)
1798 TParticle * part = (TParticle *) fParticles.At(iPart);
1799 Double_t E= part->Energy();
1800 Double_t P= part->P();
1801 M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
1804 Double_t Mmin, Mmax, wSD, wDD, wND;
1805 if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1807 if(M>-1 && M<Mmin) return kFALSE;
1810 Int_t procType=fPythia->GetMSTI(1);
1812 if(procType== 94) proc0=1;
1813 if(procType== 92 || procType== 93) proc0=0;
1817 else if(proc0==1) proc=1;
1819 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1820 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
1823 // if(proc==1 || proc==2) return kFALSE;
1826 if(proc0!=0) fProcDiff = procType;
1827 else fProcDiff = 95;
1831 if(wSD<0) AliError("wSD<0 ! \n");
1833 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1835 // printf("iPart = %d\n", iPart);
1837 if(iPart==iPart1) fProcDiff=93;
1838 else if(iPart==iPart2) fProcDiff=92;
1840 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
1849 Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1850 Double_t &wSD, Double_t &wDD, Double_t &wND)
1854 if(TMath::Abs(fEnergyCMS-900)<1 ){
1856 const Int_t nbin=400;
1858 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
1859 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
1860 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
1861 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
1862 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
1863 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
1864 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
1865 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
1866 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
1867 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
1868 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
1869 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
1870 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
1871 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
1872 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
1873 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
1874 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
1875 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
1876 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
1877 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
1878 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
1879 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
1880 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
1881 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
1882 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
1883 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
1884 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
1885 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
1886 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
1887 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
1888 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
1889 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
1890 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
1891 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1892 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1893 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1894 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1895 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1896 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1897 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1898 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1899 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1900 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
1901 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
1902 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
1903 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
1904 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
1905 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
1906 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
1907 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
1908 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
1909 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
1910 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
1911 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
1912 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
1913 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
1914 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
1915 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
1916 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
1917 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
1918 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
1919 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
1920 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
1921 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
1922 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
1923 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
1924 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1926 1.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
1927 0.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
1928 0.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
1929 0.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
1930 0.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
1931 0.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
1932 0.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
1933 0.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
1934 0.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
1935 0.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
1936 0.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
1937 0.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
1938 0.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
1939 0.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
1940 0.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
1941 0.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
1942 0.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
1943 0.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
1944 0.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
1945 0.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
1946 0.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
1947 0.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
1948 0.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
1949 0.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
1950 0.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
1951 0.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
1952 0.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
1953 0.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
1954 0.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
1955 0.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
1956 0.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
1957 0.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
1958 0.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
1959 0.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
1960 0.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
1961 0.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
1962 0.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
1963 0.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
1964 0.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
1965 0.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
1966 0.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
1967 0.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
1968 0.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
1969 0.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
1970 0.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
1971 0.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
1972 0.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
1973 0.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
1974 0.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
1975 0.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
1976 0.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
1977 0.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
1978 0.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
1979 0.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
1980 0.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
1981 0.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
1982 0.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
1983 0.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
1984 0.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
1985 0.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
1986 0.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
1987 0.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
1988 0.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
1989 0.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
1990 0.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
1991 0.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
1992 0.386112, 0.364314, 0.434375, 0.334629};
1999 if(M<Mmin || M>Mmax) return kTRUE;
2002 for(Int_t i=1; i<=nbin; i++)
2005 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2011 else if(TMath::Abs(fEnergyCMS-2760)<1 ){
2013 const Int_t nbin=400;
2015 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2016 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2017 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2018 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2019 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2020 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2021 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2022 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2023 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2024 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2025 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2026 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2027 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2028 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2029 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2030 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2031 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2032 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2033 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2034 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2035 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2036 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2037 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2038 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2039 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2040 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2041 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2042 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2043 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2044 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2045 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2046 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2047 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2048 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2049 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2050 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2051 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2052 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2053 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2054 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2055 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2056 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2057 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2058 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2059 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2060 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2061 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2062 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2063 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2064 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2065 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2066 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2067 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2068 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2069 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2070 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2071 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2072 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2073 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2074 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2075 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2076 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2077 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2078 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2079 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2080 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2081 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2083 1.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
2084 0.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
2085 0.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
2086 0.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
2087 0.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
2088 0.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
2089 0.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
2090 0.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
2091 0.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
2092 0.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
2093 0.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
2094 0.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
2095 0.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
2096 0.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
2097 0.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
2098 0.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
2099 0.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
2100 0.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
2101 0.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
2102 0.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
2103 0.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
2104 0.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
2105 0.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
2106 0.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
2107 0.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
2108 0.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
2109 0.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
2110 0.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
2111 0.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
2112 0.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
2113 0.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
2114 0.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
2115 0.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
2116 0.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
2117 0.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
2118 0.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
2119 0.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
2120 0.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
2121 0.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
2122 0.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
2123 0.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
2124 0.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
2125 0.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
2126 0.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
2127 0.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
2128 0.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
2129 0.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
2130 0.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
2131 0.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
2132 0.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
2133 0.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
2134 0.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
2135 0.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
2136 0.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
2137 0.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
2138 0.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
2139 0.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
2140 0.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
2141 0.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
2142 0.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
2143 0.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
2144 0.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
2145 0.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
2146 0.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
2147 0.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
2148 0.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
2149 0.201712, 0.242225, 0.255565, 0.258738};
2156 if(M<Mmin || M>Mmax) return kTRUE;
2159 for(Int_t i=1; i<=nbin; i++)
2162 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2170 else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2171 const Int_t nbin=400;
2173 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2174 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2175 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2176 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2177 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2178 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2179 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2180 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2181 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2182 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2183 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2184 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2185 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2186 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2187 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2188 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2189 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2190 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2191 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2192 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2193 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2194 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2195 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2196 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2197 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2198 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2199 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2200 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2201 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2202 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2203 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2204 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2205 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2206 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2207 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2208 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2209 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2210 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2211 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2212 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2213 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2214 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2215 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2216 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2217 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2218 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2219 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2220 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2221 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2222 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2223 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2224 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2225 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2226 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2227 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2228 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2229 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2230 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2231 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2232 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2233 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2234 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2235 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2236 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2237 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2238 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2239 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2241 1.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
2242 0.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
2243 0.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
2244 0.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
2245 0.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
2246 0.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
2247 0.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
2248 0.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
2249 0.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
2250 0.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
2251 0.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
2252 0.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
2253 0.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
2254 0.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
2255 0.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
2256 0.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
2257 0.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
2258 0.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
2259 0.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
2260 0.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
2261 0.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
2262 0.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
2263 0.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
2264 0.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
2265 0.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
2266 0.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
2267 0.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
2268 0.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
2269 0.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
2270 0.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
2271 0.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
2272 0.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
2273 0.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
2274 0.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
2275 0.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
2276 0.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
2277 0.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
2278 0.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
2279 0.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
2280 0.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
2281 0.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
2282 0.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
2283 0.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
2284 0.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
2285 0.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
2286 0.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
2287 0.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
2288 0.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
2289 0.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
2290 0.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
2291 0.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
2292 0.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
2293 0.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
2294 0.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
2295 0.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
2296 0.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
2297 0.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
2298 0.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
2299 0.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
2300 0.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
2301 0.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
2302 0.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
2303 0.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
2304 0.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
2305 0.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
2306 0.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
2307 0.175006, 0.223482, 0.202706, 0.218108};
2315 if(M<Mmin || M>Mmax) return kTRUE;
2318 for(Int_t i=1; i<=nbin; i++)
2321 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2331 Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2333 // Check if this is a heavy flavor decay product
2334 TParticle * part = (TParticle *) fParticles.At(ipart);
2335 Int_t mpdg = TMath::Abs(part->GetPdgCode());
2336 Int_t mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2339 if (mfl >= 4 && mfl < 6) return kTRUE;
2341 Int_t imo = part->GetFirstMother()-1;
2342 TParticle* pm = part;
2344 // Heavy flavor hadron produced by generator
2346 pm = (TParticle*)fParticles.At(imo);
2347 mpdg = TMath::Abs(pm->GetPdgCode());
2348 mfl = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2349 if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2350 imo = pm->GetFirstMother()-1;
2355 Bool_t AliGenPythia::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
2357 // check the eta/phi correspond to the detectors acceptance
2358 // iparticle is the index of the particle to be checked, for PHOS rotation case
2359 if (fCheckPHOS && IsInPHOS (phi,eta,iparticle)) return kTRUE;
2360 else if(fCheckEMCAL && IsInEMCAL (phi,eta)) return kTRUE;
2361 else if(fCheckBarrel && IsInBarrel( eta)) return kTRUE;
2365 Bool_t AliGenPythia::TriggerOnSelectedParticles(Int_t np)
2367 // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2368 // implemented primaryly for kPyJets, but extended to any kind of process.
2370 //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2371 // fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
2374 for (Int_t i=0; i< np; i++) {
2376 TParticle* iparticle = (TParticle *) fParticles.At(i);
2378 Int_t pdg = iparticle->GetPdgCode();
2379 Int_t status = iparticle->GetStatusCode();
2380 Int_t imother = iparticle->GetFirstMother() - 1;
2382 TParticle* pmother = 0x0;
2383 Int_t momStatus = -1;
2386 pmother = (TParticle *) fParticles.At(imother);
2387 momStatus = pmother->GetStatusCode();
2388 momPdg = pmother->GetPdgCode();
2394 // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2397 if (fHadronInCalo && status == 1)
2399 if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
2400 // (in case neutral mesons were declared stable)
2404 else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2408 //Fragmentation photon
2409 else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2411 if(momStatus != 11) ok = kTRUE ; // No photon from hadron decay
2414 else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2416 if( momStatus == 11)
2418 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2419 // pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2420 ok = kTRUE ; // photon from hadron decay
2422 //In case only decays from pi0 or eta requested
2423 if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2424 if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2428 // Pi0 or Eta particle
2429 else if ((fPi0InCalo || fEtaInCalo))
2431 if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2433 if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2435 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2438 else if (fEtaInCalo && pdg == 221)
2440 //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2447 // Check that the selected particle is in the calorimeter acceptance
2449 if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2451 //Just check if the selected particle falls in the acceptance
2452 if(!fForceNeutralMeson2PhotonDecay )
2454 //printf("\t Check acceptance! \n");
2455 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2456 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2458 if(CheckDetectorAcceptance(phi,eta,i))
2461 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));
2462 //printf("\t Accept \n");
2467 //Mesons have several decay modes, select only those decaying into 2 photons
2468 else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2470 // In case we want the pi0/eta trigger,
2471 // check the decay mode (2 photons)
2473 //printf("\t Force decay 2 gamma\n");
2475 Int_t ndaughters = iparticle->GetNDaughters();
2476 if(ndaughters != 2){
2481 TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2482 TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2488 //iparticle->Print();
2492 Int_t pdgD1 = d1->GetPdgCode();
2493 Int_t pdgD2 = d2->GetPdgCode();
2494 //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2495 //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2497 if(pdgD1 != 22 || pdgD2 != 22){
2502 //printf("\t accept decay\n");
2504 //Trigger on the meson, not on the daughter
2505 if(!fDecayPhotonInCalo){
2507 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2508 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2510 if(CheckDetectorAcceptance(phi,eta,i))
2512 //printf("\t Accept meson pdg %d\n",pdg);
2514 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));
2522 //printf("Check daughters acceptance\n");
2524 //Trigger on the meson daughters
2526 Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2527 Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
2528 if(d1->Pt() > fTriggerParticleMinPt)
2530 //printf("\t Check acceptance photon 1! \n");
2531 if(CheckDetectorAcceptance(phi,eta,i))
2533 //printf("\t Accept Photon 1\n");
2535 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));
2543 phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2544 eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2546 if(d2->Pt() > fTriggerParticleMinPt)
2548 //printf("\t Check acceptance photon 2! \n");
2549 if(CheckDetectorAcceptance(phi,eta,i))
2551 //printf("\t Accept Photon 2\n");
2553 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 } // force 2 photon daughters in pi0/eta decays
2561 } else ok = kFALSE; // check acceptance
2565 // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2566 // A particle passing all trigger conditions except phi position in PHOS, is used as reference