From: hristov Date: Fri, 10 Jun 2011 13:56:37 +0000 (+0000) Subject: Tuned crossections according to the ALICE measurements (Martin) X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=800be8b7a0f11e86a9b82a1f6ff7faf13cadec59 Tuned crossections according to the ALICE measurements (Martin) --- diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx index 5908d6c16d5..5b61d49d897 100644 --- a/PYTHIA6/AliGenPythia.cxx +++ b/PYTHIA6/AliGenPythia.cxx @@ -46,6 +46,7 @@ #include "AliMC.h" #include "AliLog.h" #include "PyquenCommon.h" +#include "AliLog.h" ClassImp(AliGenPythia) @@ -138,8 +139,9 @@ AliGenPythia::AliGenPythia(): fPHOSEta(0.13), fEMCALMinPhi(79.), fEMCALMaxPhi(191.), - fEMCALEta(0.71) - + fEMCALEta(0.71), + fkTuneForDiff(0), + fProcDiff(0) { // Default Constructor fEnergyCMS = 5500.; @@ -235,7 +237,9 @@ AliGenPythia::AliGenPythia(Int_t npart) fPHOSEta(0.13), fEMCALMinPhi(79.), fEMCALMaxPhi(191.), - fEMCALEta(0.71) + fEMCALEta(0.71), + fkTuneForDiff(0), + fProcDiff(0) { // default charm production at 5. 5 TeV // semimuonic decay @@ -979,6 +983,16 @@ Int_t AliGenPythia::GenerateMB() AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi)); } + // Check for diffraction + if(fkTuneForDiff) { + if(fItune==320) { + if(!CheckDiffraction()) { + delete [] pParent; + return 0; + } + } + } + // Check for minimum multiplicity if (fTriggerMultiplicity > 0) { Int_t multiplicity = 0; @@ -1236,6 +1250,12 @@ void AliGenPythia::MakeHeader() fHeader = new AliGenPythiaEventHeader("Pythia"); // // Event type + if(fProcDiff>0){ + // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n"); + // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff); + ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff); + } + else ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1)); // // Number of trials @@ -1639,7 +1659,7 @@ void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd) //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180; Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180; - Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax); + Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax); //calculate deltaphi TParticle* ph = (TParticle *) fParticles.At(iphcand); @@ -1687,3 +1707,217 @@ void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd) +Bool_t AliGenPythia::CheckDiffraction() +{ + // use this method only with Perugia-0 tune! + + // printf("AAA\n"); + + Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons; + + Int_t iPart1=-1; + Int_t iPart2=-1; + + Double_t y1 = 1e10; + Double_t y2 = -1e10; + + const Int_t kNstable=20; + const Int_t pdgStable[20] = { + 22, // Photon + 11, // Electron + 12, // Electron Neutrino + 13, // Muon + 14, // Muon Neutrino + 15, // Tau + 16, // Tau Neutrino + 211, // Pion + 321, // Kaon + 311, // K0 + 130, // K0s + 310, // K0l + 2212, // Proton + 2112, // Neutron + 3122, // Lambda_0 + 3112, // Sigma Minus + 3222, // Sigma Plus + 3312, // Xsi Minus + 3322, // Xsi0 + 3334 // Omega + }; + + for (Int_t i = 0; i < np; i++) { + TParticle * part = (TParticle *) fParticles.At(i); + + Int_t statusCode = part->GetStatusCode(); + + // Initial state particle + if (statusCode != 1) + continue; + + Int_t pdg = TMath::Abs(part->GetPdgCode()); + Bool_t isStable = kFALSE; + for (Int_t i1 = 0; i1 < kNstable; i1++) { + if (pdg == pdgStable[i1]) { + isStable = kTRUE; + break; + } + } + if(!isStable) + continue; + + Double_t y = part->Y(); + + if (y < y1) + { + y1 = y; + iPart1 = i; + } + if (y > y2) + { + y2 = y; + iPart2 = i; + } + } + + if(iPart1<0 || iPart2<0) return kFALSE; + + y1=TMath::Abs(y1); + y2=TMath::Abs(y2); + + TParticle * part1 = (TParticle *) fParticles.At(iPart1); + TParticle * part2 = (TParticle *) fParticles.At(iPart2); + + Int_t pdg1 = part1->GetPdgCode(); + Int_t pdg2 = part2->GetPdgCode(); + + + Int_t iPart = -1; + if (pdg1 == 2212 && pdg2 == 2212) + { + if(y1 > y2) + iPart = iPart1; + else if(y1 < y2) + iPart = iPart2; + else { + iPart = iPart1; + if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2; + } + } + else if (pdg1 == 2212) + iPart = iPart1; + else if (pdg2 == 2212) + iPart = iPart2; + + + + + + Double_t M=-1.; + if(iPart>0) { + TParticle * part = (TParticle *) fParticles.At(iPart); + Double_t E= part->Energy(); + Double_t P= part->P(); + M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P)); + } + +Int_t nbin=120; +Double_t bin[]={ +1.080000, 1.274258, 1.468516, 1.662773, 1.857031, 2.051289, +2.245547, 2.439805, 2.634062, 2.828320, 3.022578, 3.216836, +3.411094, 3.605352, 3.799609, 3.993867, 4.188125, 4.382383, +4.576641, 4.770898, 4.965156, 5.547930, 6.130703, 6.713477, +7.296250, 7.879023, 8.461797, 9.044570, 9.627344, 10.210117, +10.792891, 11.375664, 11.958437, 12.541211, 13.123984, 13.706758, +14.289531, 14.872305, 15.455078, 16.037852, 16.620625, 17.203398, +17.786172, 18.368945, 18.951719, 19.534492, 20.117266, 20.700039, +21.282812, 21.865586, 22.448359, 23.031133, 23.613906, 24.196680, +24.779453, 25.362227, 25.945000, 26.527773, 27.110547, 27.693320, +28.276094, 28.858867, 29.441641, 30.024414, 30.607187, 31.189961, +31.772734, 32.355508, 32.938281, 33.521055, 34.103828, 34.686602, +35.269375, 35.852148, 36.434922, 37.017695, 37.600469, 38.183242, +38.766016, 39.348789, 39.931562, 40.514336, 41.097109, 41.679883, +42.262656, 42.845430, 43.428203, 44.010977, 44.593750, 45.176523, +45.759297, 46.342070, 46.924844, 47.507617, 48.090391, 48.673164, +49.255937, 49.838711, 50.421484, 57.220508, 64.019531, 70.818555, +77.617578, 84.416602, 91.215625, 98.014648, 104.813672, 111.612695, +118.411719, 125.210742, 132.009766, 138.808789, 145.607812, 152.406836, +159.205859, 166.004883, 172.803906, 179.602930, 186.401953, 193.200977, +200.000000}; +Double_t w[]={ +1.000000, 0.275812, 0.267308, 0.263268, 0.263828, 0.263039, +0.260546, 0.259344, 0.255477, 0.253782, 0.253562, 0.252492, +0.251076, 0.247862, 0.248933, 0.243599, 0.244255, 0.238506, +0.239546, 0.237845, 0.237977, 0.229379, 0.224771, 0.217581, +0.208860, 0.207241, 0.198955, 0.196449, 0.193827, 0.190220, +0.184715, 0.183067, 0.178325, 0.175887, 0.171299, 0.168718, +0.167453, 0.165153, 0.163457, 0.159637, 0.156855, 0.155488, +0.154545, 0.155968, 0.150488, 0.148797, 0.145358, 0.146196, +0.145891, 0.142752, 0.145072, 0.141265, 0.141857, 0.138462, +0.142992, 0.141357, 0.139391, 0.139629, 0.135197, 0.135731, +0.133194, 0.137190, 0.135745, 0.134522, 0.136094, 0.130405, +0.128371, 0.131680, 0.130591, 0.133539, 0.129370, 0.128254, +0.128262, 0.131088, 0.128294, 0.130070, 0.124553, 0.131489, +0.128038, 0.122997, 0.130699, 0.125630, 0.124746, 0.123679, +0.127864, 0.125776, 0.126272, 0.121492, 0.124929, 0.122221, +0.127017, 0.123706, 0.122701, 0.123731, 0.122219, 0.120783, +0.120324, 0.120228, 0.123615, 0.120589, 0.119549, 0.118836, +0.118146, 0.120175, 0.122031, 0.122076, 0.122849, 0.122627, +0.126281, 0.127696, 0.129084, 0.130001, 0.134062, 0.134786, +0.137286, 0.136625, 0.139091, 0.143692, 0.144781, 0.146768}; + + Double_t wSD=1.; + Double_t wDD=0.178783; + //Double_t wND=0.220200; + Double_t wND=0.220200+0.001; + + if(M>-1 && Mbin[nbin]) M=-1; + + Int_t procType=fPythia->GetMSTI(1); + Int_t proc0=2; + if(procType== 94) proc0=1; + if(procType== 92 || procType== 93) proc0=0; + + + // printf("M = %f bin[nbin] = %f\n",M, bin[nbin]); + + Int_t proc=2; + if(M>0) proc=0; + else if(proc0==1) proc=1; + + if(proc==0 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wSD) return kFALSE; + if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE; + if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE; + + + // if(proc==1 || proc==2) return kFALSE; + + if(proc!=0) { + if(proc0!=0) fProcDiff = procType; + else fProcDiff = 95; + return kTRUE; + } + + Int_t ibin=-1; + for(Int_t i=0; ibin[i] && M<=bin[i+1]) { + ibin=i; + // printf("Mi> %f && Mi< %f\n", bin[i], bin[i+1]); + break; + } + + // printf("w[ibin] = %f\n", w[ibin]); + + if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> w[ibin]) return kFALSE; + + // printf("iPart = %d\n", iPart); + + if(iPart==iPart1) fProcDiff=93; + else if(iPart==iPart2) fProcDiff=92; + else { + printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n"); + + } + + return kTRUE; +} diff --git a/PYTHIA6/AliGenPythia.h b/PYTHIA6/AliGenPythia.h index 3d711fdfff4..f03376f9ab2 100644 --- a/PYTHIA6/AliGenPythia.h +++ b/PYTHIA6/AliGenPythia.h @@ -192,6 +192,8 @@ class AliGenPythia : public AliGenMC Bool_t CheckKinematicsOnChild(); void GetSubEventTime(); + void SetTuneForDiff(Bool_t a=kTRUE) {fkTuneForDiff=a;} + protected: // adjust the weight from kinematic cuts void AdjustWeights() const; @@ -299,11 +301,17 @@ class AliGenPythia : public AliGenMC Float_t fEMCALMinPhi; // Minimum phi EMCAL Float_t fEMCALMaxPhi; // Maximum phi EMCAL Float_t fEMCALEta; // Maximum eta EMCAL + + Bool_t fkTuneForDiff; // Pythia tune + Int_t fProcDiff; private: AliGenPythia(const AliGenPythia &Pythia); AliGenPythia & operator=(const AliGenPythia & rhs); - ClassDef(AliGenPythia, 10) // AliGenerator interface to Pythia + + Bool_t CheckDiffraction(); + + ClassDef(AliGenPythia, 11) // AliGenerator interface to Pythia }; #endif