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():
90 fDecayer(new AliDecayerPythia()),
98 fPhiMaxJet(2.* TMath::Pi()),
99 fJetReconstruction(kCell),
103 fPhiMaxGamma(2. * TMath::Pi()),
107 fPycellThreshold(0.),
109 fPycellMinEtJet(10.),
110 fPycellMaxRadius(1.),
111 fStackFillOpt(kFlavorSelection),
113 fFragmentation(kTRUE),
120 fTriggerMultiplicity(0),
121 fTriggerMultiplicityEta(0),
122 fTriggerMultiplicityPtMin(0),
123 fCountMode(kCountAll),
127 fFragPhotonInCalo(kFALSE),
129 fPhotonInCalo(kFALSE),
133 fCheckPHOSeta(kFALSE),
134 fFragPhotonOrPi0MinPt(0),
146 // Default Constructor
148 if (!AliPythiaRndm::GetPythiaRandom())
149 AliPythiaRndm::SetPythiaRandom(GetRandom());
152 AliGenPythia::AliGenPythia(Int_t npart)
164 fInteractionRate(0.),
178 fHadronisation(kTRUE),
179 fPatchOmegaDalitz(0),
181 fReadFromFile(kFALSE),
188 fDecayer(new AliDecayerPythia()),
189 fDebugEventFirst(-1),
196 fPhiMaxJet(2.* TMath::Pi()),
197 fJetReconstruction(kCell),
201 fPhiMaxGamma(2. * TMath::Pi()),
205 fPycellThreshold(0.),
207 fPycellMinEtJet(10.),
208 fPycellMaxRadius(1.),
209 fStackFillOpt(kFlavorSelection),
211 fFragmentation(kTRUE),
218 fTriggerMultiplicity(0),
219 fTriggerMultiplicityEta(0),
220 fTriggerMultiplicityPtMin(0),
221 fCountMode(kCountAll),
225 fFragPhotonInCalo(kFALSE),
227 fPhotonInCalo(kFALSE),
231 fCheckPHOSeta(kFALSE),
232 fFragPhotonOrPi0MinPt(0),
244 // default charm production at 5. 5 TeV
246 // structure function GRVHO
250 fTitle= "Particle Generator using PYTHIA";
252 // Set random number generator
253 if (!AliPythiaRndm::GetPythiaRandom())
254 AliPythiaRndm::SetPythiaRandom(GetRandom());
257 AliGenPythia::~AliGenPythia()
260 if(fEventsTime) delete fEventsTime;
263 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
265 // Generate pileup using user specified rate
266 fInteractionRate = rate;
267 fTimeWindow = timewindow;
271 void AliGenPythia::GeneratePileup()
273 // Generate sub events time for pileup
275 if(fInteractionRate == 0.) {
276 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
280 Int_t npart = NumberParticles();
282 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
286 if(fEventsTime) delete fEventsTime;
287 fEventsTime = new TArrayF(npart);
288 TArrayF &array = *fEventsTime;
289 for(Int_t ipart = 0; ipart < npart; ipart++)
292 Float_t eventtime = 0.;
295 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
296 if(eventtime > fTimeWindow) break;
297 array.Set(array.GetSize()+1);
298 array[array.GetSize()-1] = eventtime;
304 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
305 if(TMath::Abs(eventtime) > fTimeWindow) break;
306 array.Set(array.GetSize()+1);
307 array[array.GetSize()-1] = eventtime;
310 SetNumberParticles(fEventsTime->GetSize());
313 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
314 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
316 // Set pycell parameters
317 fPycellEtaMax = etamax;
320 fPycellThreshold = thresh;
321 fPycellEtSeed = etseed;
322 fPycellMinEtJet = minet;
323 fPycellMaxRadius = r;
328 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
330 // Set a range of event numbers, for which a table
331 // of generated particle will be printed
332 fDebugEventFirst = eventFirst;
333 fDebugEventLast = eventLast;
334 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
337 void AliGenPythia::Init()
341 SetMC(AliPythia::Instance());
342 fPythia=(AliPythia*) fMCEvGen;
345 fParentWeight=1./Float_t(fNpart);
349 fPythia->SetCKIN(3,fPtHardMin);
350 fPythia->SetCKIN(4,fPtHardMax);
351 fPythia->SetCKIN(7,fYHardMin);
352 fPythia->SetCKIN(8,fYHardMax);
354 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
356 if (fFragmentation) {
357 fPythia->SetMSTP(111,1);
359 fPythia->SetMSTP(111,0);
363 // initial state radiation
364 fPythia->SetMSTP(61,fGinit);
365 // final state radiation
366 fPythia->SetMSTP(71,fGfinal);
369 fPythia->SetMSTP(91,1);
370 fPythia->SetPARP(91,fPtKick);
371 fPythia->SetPARP(93, 4. * fPtKick);
373 fPythia->SetMSTP(91,0);
378 fRL = AliRunLoader::Open(fkFileName, "Partons");
379 fRL->LoadKinematics();
385 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
386 // Forward Paramters to the AliPythia object
387 fDecayer->SetForceDecay(fForceDecay);
388 // Switch off Heavy Flavors on request
390 // Maximum number of quark flavours used in pdf
391 fPythia->SetMSTP(58, 3);
392 // Maximum number of flavors that can be used in showers
393 fPythia->SetMSTJ(45, 3);
394 // Switch off g->QQbar splitting in decay table
395 ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
401 // Parent and Children Selection
404 case kPyOldUEQ2ordered:
405 case kPyOldUEQ2ordered2:
409 case kPyCharmUnforced:
410 case kPyCharmPbPbMNR:
413 case kPyCharmppMNRwmi:
414 fParentSelect[0] = 411;
415 fParentSelect[1] = 421;
416 fParentSelect[2] = 431;
417 fParentSelect[3] = 4122;
418 fParentSelect[4] = 4232;
419 fParentSelect[5] = 4132;
420 fParentSelect[6] = 4332;
426 fParentSelect[0] = 421;
429 case kPyDPlusPbPbMNR:
432 fParentSelect[0] = 411;
435 case kPyDPlusStrangePbPbMNR:
436 case kPyDPlusStrangepPbMNR:
437 case kPyDPlusStrangeppMNR:
438 fParentSelect[0] = 431;
441 case kPyLambdacppMNR:
442 fParentSelect[0] = 4122;
447 case kPyBeautyPbPbMNR:
448 case kPyBeautypPbMNR:
450 case kPyBeautyppMNRwmi:
451 fParentSelect[0]= 511;
452 fParentSelect[1]= 521;
453 fParentSelect[2]= 531;
454 fParentSelect[3]= 5122;
455 fParentSelect[4]= 5132;
456 fParentSelect[5]= 5232;
457 fParentSelect[6]= 5332;
460 case kPyBeautyUnforced:
461 fParentSelect[0] = 511;
462 fParentSelect[1] = 521;
463 fParentSelect[2] = 531;
464 fParentSelect[3] = 5122;
465 fParentSelect[4] = 5132;
466 fParentSelect[5] = 5232;
467 fParentSelect[6] = 5332;
472 fParentSelect[0] = 443;
475 case kPyMbAtlasTuneMC09:
477 case kPyMbWithDirectPhoton:
490 // JetFinder for Trigger
492 // Configure detector (EMCAL like)
494 fPythia->SetPARU(51, fPycellEtaMax);
495 fPythia->SetMSTU(51, fPycellNEta);
496 fPythia->SetMSTU(52, fPycellNPhi);
498 // Configure Jet Finder
500 fPythia->SetPARU(58, fPycellThreshold);
501 fPythia->SetPARU(52, fPycellEtSeed);
502 fPythia->SetPARU(53, fPycellMinEtJet);
503 fPythia->SetPARU(54, fPycellMaxRadius);
504 fPythia->SetMSTU(54, 2);
506 // This counts the total number of calls to Pyevnt() per run.
521 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
524 fPythia->SetPARJ(200, 0.0);
525 fPythia->SetPARJ(199, 0.0);
526 fPythia->SetPARJ(198, 0.0);
527 fPythia->SetPARJ(197, 0.0);
530 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
534 // Nestor's change of the splittings
535 fPythia->SetPARJ(200, 0.8);
536 fPythia->SetMSTJ(41, 1); // QCD radiation only
537 fPythia->SetMSTJ(42, 2); // angular ordering
538 fPythia->SetMSTJ(44, 2); // option to run alpha_s
539 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
540 fPythia->SetMSTJ(50, 0); // No coherence in first branching
541 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
542 } else if (fQuench == 4) {
543 // Armesto-Cunqueiro-Salgado change of the splittings.
544 AliFastGlauber* glauber = AliFastGlauber::Instance();
546 //read and store transverse almonds corresponding to differnt
548 glauber->SetCentralityClass(0.,0.1);
549 fPythia->SetPARJ(200, 1.);
550 fPythia->SetPARJ(198, fQhat);
551 fPythia->SetPARJ(199, fLength);
552 fPythia->SetMSTJ(42, 2); // angular ordering
553 fPythia->SetMSTJ(44, 2); // option to run alpha_s
554 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
558 void AliGenPythia::Generate()
560 // Generate one event
561 if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
562 fDecayer->ForceDecay();
564 Double_t polar[3] = {0,0,0};
565 Double_t origin[3] = {0,0,0};
567 // converts from mm/c to s
568 const Float_t kconv=0.001/2.999792458e8;
578 // Set collision vertex position
579 if (fVertexSmear == kPerEvent) Vertex();
588 // Switch hadronisation off
590 fPythia->SetMSTJ(1, 0);
594 // Quenching comes through medium-modified splitting functions.
595 AliFastGlauber::Instance()->GetRandomBHard(bimp);
596 fPythia->SetPARJ(197, bimp);
601 // Either produce new event or read partons from file
603 if (!fReadFromFile) {
609 fNpartons = fPythia->GetN();
611 printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
612 fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
614 LoadEvent(fRL->Stack(), 0 , 1);
619 // Run quenching routine
623 } else if (fQuench == 2){
624 fPythia->Pyquen(208., 0, 0.);
625 } else if (fQuench == 3) {
626 // Quenching is via multiplicative correction of the splittings
630 // Switch hadronisation on
632 if (fHadronisation) {
633 fPythia->SetMSTJ(1, 1);
635 // .. and perform hadronisation
636 // printf("Calling hadronisation %d\n", fPythia->GetN());
638 if (fPatchOmegaDalitz) {
639 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
641 fPythia->DalitzDecays();
642 fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
647 fPythia->ImportParticles(&fParticles,"All");
649 if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
650 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
658 Int_t np = fParticles.GetEntriesFast();
660 if (np == 0) continue;
664 Int_t* pParent = new Int_t[np];
665 Int_t* pSelected = new Int_t[np];
666 Int_t* trackIt = new Int_t[np];
667 for (i = 0; i < np; i++) {
673 Int_t nc = 0; // Total n. of selected particles
674 Int_t nParents = 0; // Selected parents
675 Int_t nTkbles = 0; // Trackable particles
676 if (fProcess != kPyMbDefault &&
678 fProcess != kPyMbAtlasTuneMC09 &&
679 fProcess != kPyMbWithDirectPhoton &&
680 fProcess != kPyJets &&
681 fProcess != kPyDirectGamma &&
682 fProcess != kPyMbNonDiffr &&
683 fProcess != kPyMbMSEL1 &&
686 fProcess != kPyCharmppMNRwmi &&
687 fProcess != kPyBeautyppMNRwmi &&
688 fProcess != kPyBeautyJets) {
690 for (i = 0; i < np; i++) {
691 TParticle* iparticle = (TParticle *) fParticles.At(i);
692 Int_t ks = iparticle->GetStatusCode();
693 kf = CheckPDGCode(iparticle->GetPdgCode());
694 // No initial state partons
695 if (ks==21) continue;
697 // Heavy Flavor Selection
704 if (kfl > 100000) kfl %= 100000;
705 if (kfl > 10000) kfl %= 10000;
707 if (kfl > 10) kfl/=100;
709 if (kfl > 10) kfl/=10;
710 Int_t ipa = iparticle->GetFirstMother()-1;
713 // Establish mother daughter relation between heavy quarks and mesons
715 if (kf >= fFlavorSelect && kf <= 6) {
716 Int_t idau = iparticle->GetFirstDaughter() - 1;
718 TParticle* daughter = (TParticle *) fParticles.At(idau);
719 Int_t pdgD = daughter->GetPdgCode();
720 if (pdgD == 91 || pdgD == 92) {
721 Int_t jmin = daughter->GetFirstDaughter() - 1;
722 Int_t jmax = daughter->GetLastDaughter() - 1;
723 for (Int_t jp = jmin; jp <= jmax; jp++)
724 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
725 } // is string or cluster
731 TParticle * mother = (TParticle *) fParticles.At(ipa);
732 kfMo = TMath::Abs(mother->GetPdgCode());
735 // What to keep in Stack?
736 Bool_t flavorOK = kFALSE;
737 Bool_t selectOK = kFALSE;
739 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
741 if (kfl > fFlavorSelect) {
745 if (kfl == fFlavorSelect) flavorOK = kTRUE;
747 switch (fStackFillOpt) {
748 case kFlavorSelection:
751 case kParentSelection:
752 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
755 if (flavorOK && selectOK) {
757 // Heavy flavor hadron or quark
759 // Kinematic seletion on final state heavy flavor mesons
760 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
765 if (ParentSelected(kf)) ++nParents; // Update parent count
766 // printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
768 // Kinematic seletion on decay products
769 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
770 && !KinematicSelection(iparticle, 1))
776 // Select if mother was selected and is not tracked
778 if (pSelected[ipa] &&
779 !trackIt[ipa] && // mother will be tracked ?
780 kfMo != 5 && // mother is b-quark, don't store fragments
781 kfMo != 4 && // mother is c-quark, don't store fragments
782 kf != 92) // don't store string
785 // Semi-stable or de-selected: diselect decay products:
788 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
790 Int_t ipF = iparticle->GetFirstDaughter();
791 Int_t ipL = iparticle->GetLastDaughter();
792 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
794 // printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
795 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
798 if (pSelected[i] == -1) pSelected[i] = 0;
799 if (!pSelected[i]) continue;
800 // Count quarks only if you did not include fragmentation
801 if (fFragmentation && kf <= 10) continue;
804 // Decision on tracking
807 // Track final state particle
808 if (ks == 1) trackIt[i] = 1;
809 // Track semi-stable particles
810 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
811 // Track particles selected by process if undecayed.
812 if (fForceDecay == kNoDecay) {
813 if (ParentSelected(kf)) trackIt[i] = 1;
815 if (ParentSelected(kf)) trackIt[i] = 0;
817 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
821 } // particle selection loop
823 for (i = 0; i<np; i++) {
824 if (!pSelected[i]) continue;
825 TParticle * iparticle = (TParticle *) fParticles.At(i);
826 kf = CheckPDGCode(iparticle->GetPdgCode());
827 Int_t ks = iparticle->GetStatusCode();
828 p[0] = iparticle->Px();
829 p[1] = iparticle->Py();
830 p[2] = iparticle->Pz();
831 p[3] = iparticle->Energy();
833 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
834 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
835 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
837 Float_t tof = kconv*iparticle->T();
838 Int_t ipa = iparticle->GetFirstMother()-1;
839 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
841 PushTrack(fTrackIt*trackIt[i], iparent, kf,
842 p[0], p[1], p[2], p[3],
843 origin[0], origin[1], origin[2], tof,
844 polar[0], polar[1], polar[2],
845 kPPrimary, nt, 1., ks);
862 switch (fCountMode) {
864 // printf(" Count all \n");
868 // printf(" Count parents \n");
871 case kCountTrackables:
872 // printf(" Count trackable \n");
876 if (jev >= fNpart || fNpart == -1) {
877 fKineBias=Float_t(fNpart)/Float_t(fTrials);
879 fQ += fPythia->GetVINT(51);
880 fX1 += fPythia->GetVINT(41);
881 fX2 += fPythia->GetVINT(42);
882 fTrialsRun += fTrials;
889 SetHighWaterMark(nt);
890 // adjust weight due to kinematic selection
893 fXsection=fPythia->GetPARI(1);
896 Int_t AliGenPythia::GenerateMB()
899 // Min Bias selection and other global selections
901 Int_t i, kf, nt, iparent;
904 Double_t polar[3] = {0,0,0};
905 Double_t origin[3] = {0,0,0};
906 // converts from mm/c to s
907 const Float_t kconv=0.001/2.999792458e8;
911 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
915 Int_t* pParent = new Int_t[np];
916 for (i=0; i< np; i++) pParent[i] = -1;
917 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
918 TParticle* jet1 = (TParticle *) fParticles.At(6);
919 TParticle* jet2 = (TParticle *) fParticles.At(7);
920 if (!CheckTrigger(jet1, jet2)) {
926 // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
927 if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
932 if (fFragPhotonInCalo) pdg = 22 ; // Photon
933 else if (fPi0InCalo) pdg = 111 ; // Pi0
935 for (i=0; i< np; i++) {
936 TParticle* iparticle = (TParticle *) fParticles.At(i);
937 if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
938 iparticle->Pt() > fFragPhotonOrPi0MinPt){
939 Int_t imother = iparticle->GetFirstMother() - 1;
940 TParticle* pmother = (TParticle *) fParticles.At(imother);
942 (pdg == 22 && pmother->GetStatusCode() != 11)) //No photon from hadron decay
944 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
945 Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
946 if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
947 (fCheckPHOS && IsInPHOS(phi,eta)) )
958 // Select beauty jets with electron in EMCAL
959 if (fProcess == kPyBeautyJets && fEleInEMCAL) {
963 Int_t pdg = 11; //electron
968 for (i=0; i< np; i++) {
969 TParticle* iparticle = (TParticle *) fParticles.At(i);
970 if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
971 iparticle->Pt() > fElectronMinPt){
972 pt = iparticle->Pt();
973 phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
974 eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
975 if(IsInEMCAL(phi,eta))
984 AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
986 // Check for diffraction
989 if(!CheckDiffraction()) {
996 // Check for minimum multiplicity
997 if (fTriggerMultiplicity > 0) {
998 Int_t multiplicity = 0;
999 for (i = 0; i < np; i++) {
1000 TParticle * iparticle = (TParticle *) fParticles.At(i);
1002 Int_t statusCode = iparticle->GetStatusCode();
1004 // Initial state particle
1005 if (statusCode != 1)
1008 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1011 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1014 TParticlePDG* pdgPart = iparticle->GetPDG();
1015 if (pdgPart && pdgPart->Charge() == 0)
1021 if (multiplicity < fTriggerMultiplicity) {
1025 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
1028 // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
1029 if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
1031 Bool_t okd = kFALSE;
1035 for (i=0; i< np; i++) {
1036 TParticle* iparticle = (TParticle *) fParticles.At(i);
1037 Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1038 Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
1040 if(iparticle->GetStatusCode() == 1
1041 && iparticle->GetPdgCode() == pdg
1042 && iparticle->Pt() > fPhotonMinPt
1045 // first check if the photon is in PHOS phi
1046 if(IsInPHOS(phi,eta)){
1050 if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
1055 if(!okd && iphcand != -1) // execute rotation in phi
1056 RotatePhi(iphcand,okd);
1064 if (fTriggerParticle) {
1065 Bool_t triggered = kFALSE;
1066 for (i = 0; i < np; i++) {
1067 TParticle * iparticle = (TParticle *) fParticles.At(i);
1068 kf = CheckPDGCode(iparticle->GetPdgCode());
1069 if (kf != fTriggerParticle) continue;
1070 if (iparticle->Pt() == 0.) continue;
1071 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1082 // Check if there is a ccbar or bbbar pair with at least one of the two
1083 // in fYMin < y < fYMax
1085 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1086 TParticle *partCheck;
1088 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1089 Bool_t theChild=kFALSE;
1090 Bool_t triggered=kFALSE;
1092 Int_t pdg,mpdg,mpdgUpperFamily;
1093 for(i=0; i<np; i++) {
1094 partCheck = (TParticle*)fParticles.At(i);
1095 pdg = partCheck->GetPdgCode();
1096 if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1097 if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1098 y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1099 (partCheck->Energy()-partCheck->Pz()+1.e-13));
1100 if(y>fYMin && y<fYMax) inYcut=kTRUE;
1102 if(fTriggerParticle) {
1103 if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1105 if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1106 Int_t mi = partCheck->GetFirstMother() - 1;
1108 mother = (TParticle*)fParticles.At(mi);
1109 mpdg=TMath::Abs(mother->GetPdgCode());
1110 mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1111 if ( ParentSelected(mpdg) ||
1112 (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1113 if (KinematicSelection(partCheck,1)) {
1119 if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1123 if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1127 if(fTriggerParticle && !triggered) { // particle requested is not produced
1134 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1135 if ( (fProcess == kPyW ||
1137 fProcess == kPyMbDefault ||
1138 fProcess == kPyMb ||
1139 fProcess == kPyMbAtlasTuneMC09 ||
1140 fProcess == kPyMbWithDirectPhoton ||
1141 fProcess == kPyMbNonDiffr)
1142 && (fCutOnChild == 1) ) {
1143 if ( !CheckKinematicsOnChild() ) {
1150 for (i = 0; i < np; i++) {
1152 TParticle * iparticle = (TParticle *) fParticles.At(i);
1153 kf = CheckPDGCode(iparticle->GetPdgCode());
1154 Int_t ks = iparticle->GetStatusCode();
1155 Int_t km = iparticle->GetFirstMother();
1156 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
1158 ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
1160 if (ks == 1) trackIt = 1;
1161 Int_t ipa = iparticle->GetFirstMother()-1;
1163 iparent = (ipa > -1) ? pParent[ipa] : -1;
1166 // store track information
1167 p[0] = iparticle->Px();
1168 p[1] = iparticle->Py();
1169 p[2] = iparticle->Pz();
1170 p[3] = iparticle->Energy();
1173 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1174 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1175 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1177 Float_t tof = fEventTime + kconv * iparticle->T();
1179 PushTrack(fTrackIt*trackIt, iparent, kf,
1180 p[0], p[1], p[2], p[3],
1181 origin[0], origin[1], origin[2], tof,
1182 polar[0], polar[1], polar[2],
1183 kPPrimary, nt, 1., ks);
1187 SetHighWaterMark(nt);
1189 } // select particle
1198 void AliGenPythia::FinishRun()
1200 // Print x-section summary
1209 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1210 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1213 void AliGenPythia::AdjustWeights() const
1215 // Adjust the weights after generation of all events
1219 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1220 for (Int_t i=0; i<ntrack; i++) {
1221 part= gAlice->GetMCApp()->Particle(i);
1222 part->SetWeight(part->GetWeight()*fKineBias);
1227 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1229 // Treat protons as inside nuclei with mass numbers a1 and a2
1233 fNucPdf = pdfset; // 0 EKS98 1 EPS08
1238 void AliGenPythia::MakeHeader()
1241 // Make header for the simulated event
1244 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1245 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1248 // Builds the event header, to be called after each event
1249 if (fHeader) delete fHeader;
1250 fHeader = new AliGenPythiaEventHeader("Pythia");
1254 // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1255 // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1256 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1259 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1262 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1265 fHeader->SetPrimaryVertex(fVertex);
1266 fHeader->SetInteractionTime(fEventTime);
1268 // Number of primaries
1269 fHeader->SetNProduced(fNprimaries);
1271 // Jets that have triggered
1273 //Need to store jets for b-jet studies too!
1274 if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
1277 Float_t jets[4][10];
1278 GetJets(njet, ntrig, jets);
1281 for (Int_t i = 0; i < ntrig; i++) {
1282 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1287 // Copy relevant information from external header, if present.
1292 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1293 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1295 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1298 exHeader->TriggerJet(i, uqJet);
1299 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1303 // Store quenching parameters
1306 Double_t z[4] = {0.};
1312 fPythia->GetQuenchingParameters(xp, yp, z);
1313 } else if (fQuench == 2){
1315 Double_t r1 = PARIMP.rb1;
1316 Double_t r2 = PARIMP.rb2;
1317 Double_t b = PARIMP.b1;
1318 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1319 Double_t phi = PARIMP.psib1;
1320 xp = r * TMath::Cos(phi);
1321 yp = r * TMath::Sin(phi);
1323 } else if (fQuench == 4) {
1327 AliFastGlauber::Instance()->GetSavedXY(xy);
1328 AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1331 ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1334 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1335 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1339 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1347 Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1349 // Check the kinematic trigger condition
1352 eta[0] = jet1->Eta();
1353 eta[1] = jet2->Eta();
1355 phi[0] = jet1->Phi();
1356 phi[1] = jet2->Phi();
1358 pdg[0] = jet1->GetPdgCode();
1359 pdg[1] = jet2->GetPdgCode();
1360 Bool_t triggered = kFALSE;
1362 if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
1365 Float_t jets[4][10];
1367 // Use Pythia clustering on parton level to determine jet axis
1369 GetJets(njets, ntrig, jets);
1371 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1376 if (pdg[0] == kGamma) {
1380 //Check eta range first...
1381 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1382 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1384 //Eta is okay, now check phi range
1385 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1386 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1397 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1399 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1401 Bool_t checking = kFALSE;
1402 Int_t j, kcode, ks, km;
1403 Int_t nPartAcc = 0; //number of particles in the acceptance range
1404 Int_t numberOfAcceptedParticles = 1;
1405 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1406 Int_t npart = fParticles.GetEntriesFast();
1408 for (j = 0; j<npart; j++) {
1409 TParticle * jparticle = (TParticle *) fParticles.At(j);
1410 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1411 ks = jparticle->GetStatusCode();
1412 km = jparticle->GetFirstMother();
1414 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1417 if( numberOfAcceptedParticles <= nPartAcc){
1426 void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1429 // Load event into Pythia Common Block
1432 Int_t npart = stack -> GetNprimary();
1436 (fPythia->GetPyjets())->N = npart;
1438 n0 = (fPythia->GetPyjets())->N;
1439 (fPythia->GetPyjets())->N = n0 + npart;
1443 for (Int_t part = 0; part < npart; part++) {
1444 TParticle *mPart = stack->Particle(part);
1446 Int_t kf = mPart->GetPdgCode();
1447 Int_t ks = mPart->GetStatusCode();
1448 Int_t idf = mPart->GetFirstDaughter();
1449 Int_t idl = mPart->GetLastDaughter();
1452 if (ks == 11 || ks == 12) {
1459 Float_t px = mPart->Px();
1460 Float_t py = mPart->Py();
1461 Float_t pz = mPart->Pz();
1462 Float_t e = mPart->Energy();
1463 Float_t m = mPart->GetCalcMass();
1466 (fPythia->GetPyjets())->P[0][part+n0] = px;
1467 (fPythia->GetPyjets())->P[1][part+n0] = py;
1468 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1469 (fPythia->GetPyjets())->P[3][part+n0] = e;
1470 (fPythia->GetPyjets())->P[4][part+n0] = m;
1472 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1473 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1474 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1475 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1476 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1480 void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1483 // Load event into Pythia Common Block
1486 Int_t npart = stack -> GetEntries();
1490 (fPythia->GetPyjets())->N = npart;
1492 n0 = (fPythia->GetPyjets())->N;
1493 (fPythia->GetPyjets())->N = n0 + npart;
1497 for (Int_t part = 0; part < npart; part++) {
1498 TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1499 if (!mPart) continue;
1501 Int_t kf = mPart->GetPdgCode();
1502 Int_t ks = mPart->GetStatusCode();
1503 Int_t idf = mPart->GetFirstDaughter();
1504 Int_t idl = mPart->GetLastDaughter();
1507 if (ks == 11 || ks == 12) {
1514 Float_t px = mPart->Px();
1515 Float_t py = mPart->Py();
1516 Float_t pz = mPart->Pz();
1517 Float_t e = mPart->Energy();
1518 Float_t m = mPart->GetCalcMass();
1521 (fPythia->GetPyjets())->P[0][part+n0] = px;
1522 (fPythia->GetPyjets())->P[1][part+n0] = py;
1523 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1524 (fPythia->GetPyjets())->P[3][part+n0] = e;
1525 (fPythia->GetPyjets())->P[4][part+n0] = m;
1527 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1528 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1529 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1530 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1531 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1536 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1539 // Calls the Pythia jet finding algorithm to find jets in the current event
1544 Int_t n = fPythia->GetN();
1548 fPythia->Pycell(njets);
1550 for (i = 0; i < njets; i++) {
1551 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1552 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1553 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1554 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1565 void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1568 // Calls the Pythia clustering algorithm to find jets in the current event
1570 Int_t n = fPythia->GetN();
1573 if (fJetReconstruction == kCluster) {
1575 // Configure cluster algorithm
1577 fPythia->SetPARU(43, 2.);
1578 fPythia->SetMSTU(41, 1);
1580 // Call cluster algorithm
1582 fPythia->Pyclus(nJets);
1584 // Loading jets from common block
1590 fPythia->Pycell(nJets);
1594 for (i = 0; i < nJets; i++) {
1595 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1596 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1597 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1598 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1599 Float_t pt = TMath::Sqrt(px * px + py * py);
1600 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
1601 Float_t theta = TMath::ATan2(pt,pz);
1602 Float_t et = e * TMath::Sin(theta);
1603 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
1605 eta > fEtaMinJet && eta < fEtaMaxJet &&
1606 phi > fPhiMinJet && phi < fPhiMaxJet &&
1607 et > fEtMinJet && et < fEtMaxJet
1610 jets[0][nJetsTrig] = px;
1611 jets[1][nJetsTrig] = py;
1612 jets[2][nJetsTrig] = pz;
1613 jets[3][nJetsTrig] = e;
1615 // printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1617 // printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1622 void AliGenPythia::GetSubEventTime()
1624 // Calculates time of the next subevent
1627 TArrayF &array = *fEventsTime;
1628 fEventTime = array[fCurSubEvent++];
1630 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1634 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
1636 // Is particle in EMCAL acceptance?
1637 // phi in degrees, etamin=-etamax
1638 if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
1645 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta) const
1647 // Is particle in PHOS acceptance?
1648 // Acceptance slightly larger considered.
1649 // phi in degrees, etamin=-etamax
1650 if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
1657 void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1659 //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1660 Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1661 Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1662 Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
1664 //calculate deltaphi
1665 TParticle* ph = (TParticle *) fParticles.At(iphcand);
1666 Double_t phphi = ph->Phi();
1667 Double_t deltaphi = phiPHOS - phphi;
1671 //loop for all particles and produce the phi rotation
1672 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1673 Double_t oldphi, newphi;
1674 Double_t newVx, newVy, r, vZ, time;
1675 Double_t newPx, newPy, pt, pz, e;
1676 for(Int_t i=0; i< np; i++) {
1677 TParticle* iparticle = (TParticle *) fParticles.At(i);
1678 oldphi = iparticle->Phi();
1679 newphi = oldphi + deltaphi;
1680 if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1681 if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1684 newVx = r * TMath::Cos(newphi);
1685 newVy = r * TMath::Sin(newphi);
1686 vZ = iparticle->Vz(); // don't transform
1687 time = iparticle->T(); // don't transform
1689 pt = iparticle->Pt();
1690 newPx = pt * TMath::Cos(newphi);
1691 newPy = pt * TMath::Sin(newphi);
1692 pz = iparticle->Pz(); // don't transform
1693 e = iparticle->Energy(); // don't transform
1696 iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1697 iparticle->SetMomentum(newPx, newPy, pz, e);
1699 } //end particle loop
1701 // now let's check that we put correctly the candidate photon in PHOS
1702 Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1703 Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1704 if(IsInPHOS(phi,eta))
1710 Bool_t AliGenPythia::CheckDiffraction()
1712 // use this method only with Perugia-0 tune!
1716 Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1722 Double_t y2 = -1e10;
1724 const Int_t kNstable=20;
1725 const Int_t pdgStable[20] = {
1728 12, // Electron Neutrino
1730 14, // Muon Neutrino
1741 3112, // Sigma Minus
1748 for (Int_t i = 0; i < np; i++) {
1749 TParticle * part = (TParticle *) fParticles.At(i);
1751 Int_t statusCode = part->GetStatusCode();
1753 // Initial state particle
1754 if (statusCode != 1)
1757 Int_t pdg = TMath::Abs(part->GetPdgCode());
1758 Bool_t isStable = kFALSE;
1759 for (Int_t i1 = 0; i1 < kNstable; i1++) {
1760 if (pdg == pdgStable[i1]) {
1768 Double_t y = part->Y();
1782 if(iPart1<0 || iPart2<0) return kFALSE;
1787 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
1788 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
1790 Int_t pdg1 = part1->GetPdgCode();
1791 Int_t pdg2 = part2->GetPdgCode();
1795 if (pdg1 == 2212 && pdg2 == 2212)
1803 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1806 else if (pdg1 == 2212)
1808 else if (pdg2 == 2212)
1817 TParticle * part = (TParticle *) fParticles.At(iPart);
1818 Double_t E= part->Energy();
1819 Double_t P= part->P();
1820 M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
1825 1.080000, 1.274258, 1.468516, 1.662773, 1.857031, 2.051289,
1826 2.245547, 2.439805, 2.634062, 2.828320, 3.022578, 3.216836,
1827 3.411094, 3.605352, 3.799609, 3.993867, 4.188125, 4.382383,
1828 4.576641, 4.770898, 4.965156, 5.547930, 6.130703, 6.713477,
1829 7.296250, 7.879023, 8.461797, 9.044570, 9.627344, 10.210117,
1830 10.792891, 11.375664, 11.958437, 12.541211, 13.123984, 13.706758,
1831 14.289531, 14.872305, 15.455078, 16.037852, 16.620625, 17.203398,
1832 17.786172, 18.368945, 18.951719, 19.534492, 20.117266, 20.700039,
1833 21.282812, 21.865586, 22.448359, 23.031133, 23.613906, 24.196680,
1834 24.779453, 25.362227, 25.945000, 26.527773, 27.110547, 27.693320,
1835 28.276094, 28.858867, 29.441641, 30.024414, 30.607187, 31.189961,
1836 31.772734, 32.355508, 32.938281, 33.521055, 34.103828, 34.686602,
1837 35.269375, 35.852148, 36.434922, 37.017695, 37.600469, 38.183242,
1838 38.766016, 39.348789, 39.931562, 40.514336, 41.097109, 41.679883,
1839 42.262656, 42.845430, 43.428203, 44.010977, 44.593750, 45.176523,
1840 45.759297, 46.342070, 46.924844, 47.507617, 48.090391, 48.673164,
1841 49.255937, 49.838711, 50.421484, 57.220508, 64.019531, 70.818555,
1842 77.617578, 84.416602, 91.215625, 98.014648, 104.813672, 111.612695,
1843 118.411719, 125.210742, 132.009766, 138.808789, 145.607812, 152.406836,
1844 159.205859, 166.004883, 172.803906, 179.602930, 186.401953, 193.200977,
1847 1.000000, 0.275812, 0.267308, 0.263268, 0.263828, 0.263039,
1848 0.260546, 0.259344, 0.255477, 0.253782, 0.253562, 0.252492,
1849 0.251076, 0.247862, 0.248933, 0.243599, 0.244255, 0.238506,
1850 0.239546, 0.237845, 0.237977, 0.229379, 0.224771, 0.217581,
1851 0.208860, 0.207241, 0.198955, 0.196449, 0.193827, 0.190220,
1852 0.184715, 0.183067, 0.178325, 0.175887, 0.171299, 0.168718,
1853 0.167453, 0.165153, 0.163457, 0.159637, 0.156855, 0.155488,
1854 0.154545, 0.155968, 0.150488, 0.148797, 0.145358, 0.146196,
1855 0.145891, 0.142752, 0.145072, 0.141265, 0.141857, 0.138462,
1856 0.142992, 0.141357, 0.139391, 0.139629, 0.135197, 0.135731,
1857 0.133194, 0.137190, 0.135745, 0.134522, 0.136094, 0.130405,
1858 0.128371, 0.131680, 0.130591, 0.133539, 0.129370, 0.128254,
1859 0.128262, 0.131088, 0.128294, 0.130070, 0.124553, 0.131489,
1860 0.128038, 0.122997, 0.130699, 0.125630, 0.124746, 0.123679,
1861 0.127864, 0.125776, 0.126272, 0.121492, 0.124929, 0.122221,
1862 0.127017, 0.123706, 0.122701, 0.123731, 0.122219, 0.120783,
1863 0.120324, 0.120228, 0.123615, 0.120589, 0.119549, 0.118836,
1864 0.118146, 0.120175, 0.122031, 0.122076, 0.122849, 0.122627,
1865 0.126281, 0.127696, 0.129084, 0.130001, 0.134062, 0.134786,
1866 0.137286, 0.136625, 0.139091, 0.143692, 0.144781, 0.146768};
1869 Double_t wDD=0.178783;
1870 //Double_t wND=0.220200;
1871 Double_t wND=0.220200+0.001;
1873 if(M>-1 && M<bin[0]) return kFALSE;
1874 if(M>bin[nbin]) M=-1;
1876 Int_t procType=fPythia->GetMSTI(1);
1878 if(procType== 94) proc0=1;
1879 if(procType== 92 || procType== 93) proc0=0;
1882 // printf("M = %f bin[nbin] = %f\n",M, bin[nbin]);
1886 else if(proc0==1) proc=1;
1888 if(proc==0 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wSD) return kFALSE;
1889 if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1890 if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
1893 // if(proc==1 || proc==2) return kFALSE;
1896 if(proc0!=0) fProcDiff = procType;
1897 else fProcDiff = 95;
1902 for(Int_t i=0; i<nbin; i++)
1903 if(M>bin[i] && M<=bin[i+1]) {
1905 // printf("Mi> %f && Mi< %f\n", bin[i], bin[i+1]);
1909 // printf("w[ibin] = %f\n", w[ibin]);
1911 if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> w[ibin]) return kFALSE;
1913 // printf("iPart = %d\n", iPart);
1915 if(iPart==iPart1) fProcDiff=93;
1916 else if(iPart==iPart2) fProcDiff=92;
1918 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");