Tuned crossections according to the ALICE measurements (Martin)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 10 Jun 2011 13:56:37 +0000 (13:56 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 10 Jun 2011 13:56:37 +0000 (13:56 +0000)
PYTHIA6/AliGenPythia.cxx
PYTHIA6/AliGenPythia.h

index 5908d6c..5b61d49 100644 (file)
@@ -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 && M<bin[0]) return kFALSE;
+ if(M>bin[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; i<nbin; i++) 
+      if(M>bin[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;
+}
index 3d711fd..f03376f 100644 (file)
@@ -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