]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliPythia6.cxx
Using detector quality flag (taken from ALICE logbook) to decide whether to rpodcue...
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia6.cxx
index b0b126d0e730ae00d2f4d4cd77da0b52c2fee603..6c9ff6f8a8374755a335701fa5b2cc37fa159bcd 100644 (file)
@@ -1,4 +1,3 @@
-
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
@@ -32,25 +31,38 @@ ClassImp(AliPythia6)
 # define pyclus pyclus_
 # define pycell pycell_
 # define pyshow pyshow_
+# define pyshowq pyshowq_
 # define pyrobo pyrobo_
 # define pyquen pyquen_
 # define pyevnw pyevnw_
+# define pyjoin pyjoin_
+# define qpygin0 qpygin0_
+# define setpowwght setpowwght_
 # define type_of_call
 #else
 # define pyclus PYCLUS
 # define pycell PYCELL
+# define pyshow PYSHOW
+# define pyshowq PYSHOWQ
 # define pyrobo PYROBO
 # define pyquen PYQUEN
 # define pyevnw PYEVNW
+# define pyjoin PYJOIN
+# define qpygin0 QPYGIN0
+# define setpowwght SETPOWWGHT
 # define type_of_call _stdcall
 #endif
 
+extern "C" void type_of_call pyjoin(Int_t &, Int_t * );
 extern "C" void type_of_call pyclus(Int_t & );
 extern "C" void type_of_call pycell(Int_t & );
 extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
+extern "C" void type_of_call pyshowq(Int_t &, Int_t &, Double_t &);
+extern "C" void type_of_call qpygin0();
 extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
 extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
 extern "C" void type_of_call pyevnw();
+extern "C" void type_of_call setpowwght(Double_t &);
 
 
 //_____________________________________________________________________________
@@ -63,6 +75,8 @@ AliPythia6::AliPythia6():
     fProcess(kPyMb),
     fEcms(0.),
     fStrucFunc(kCTEQ5L),
+    fProjectile("p"),
+    fTarget("p"),
     fXJet(0.),
     fYJet(0.),
     fNGmax(30),
@@ -73,6 +87,11 @@ AliPythia6::AliPythia6():
 // Default Constructor
 //
 //  Set random number
+    Int_t i;
+    for (i = 0; i <  501; i++) fDefMDCY[i] = 0;
+    for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
+    for (i = 0; i <    4; i++) fZQuench[i] = 0;
+
     if (!AliPythiaRndm::GetPythiaRandom()) 
       AliPythiaRndm::SetPythiaRandom(GetRandom());
     fGlauber          = 0;
@@ -85,6 +104,8 @@ AliPythia6::AliPythia6(const AliPythia6& pythia):
     fProcess(kPyMb),
     fEcms(0.),
     fStrucFunc(kCTEQ5L),
+    fProjectile("p"),
+    fTarget("p"),
     fXJet(0.),
     fYJet(0.),
     fNGmax(30),
@@ -93,10 +114,14 @@ AliPythia6::AliPythia6(const AliPythia6& pythia):
     fQuenchingWeights(0)
 {
     // Copy Constructor
+    Int_t i;
+    for (i = 0; i <  501; i++) fDefMDCY[i] = 0;
+    for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
+    for (i = 0; i <    4; i++) fZQuench[i] = 0;
     pythia.Copy(*this);
 }
 
-void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
+void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc, Int_t /*tune*/)
 {
 // Initialise the process to generate 
     if (!AliPythiaRndm::GetPythiaRandom()) 
@@ -249,6 +274,19 @@ void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfu
 
        AtlasTuning();
        break;
+    case kPyMbAtlasTuneMC09:
+// Minimum Bias pp-Collisions
+//
+//   
+//      select Pythia min. bias model
+       SetMSEL(0);
+       SetMSUB(92,1);             // single diffraction AB-->XB
+       SetMSUB(93,1);             // single diffraction AB-->AX
+       SetMSUB(94,1);             // double diffraction
+       SetMSUB(95,1);             // low pt production
+
+       AtlasTuningMC09();
+       break;
 
     case kPyMbWithDirectPhoton:
 // Minimum Bias pp-Collisions with direct photon processes added 
@@ -392,6 +430,7 @@ void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfu
     case kPyD0ppMNR:
     case kPyDPlusppMNR:
     case kPyDPlusStrangeppMNR:
+    case kPyLambdacppMNR:
       // Tuning of Pythia parameters aimed to get a resonable agreement
       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
       // c-cbar single inclusive and double differential distributions.
@@ -494,6 +533,7 @@ void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfu
        // Set b-quark mass
        SetPMAS(5,1,4.75);
       break;
+     case kPyBeautyJets: 
      case kPyBeautyppMNRwmi:
       // Tuning of Pythia parameters aimed to get a resonable agreement
       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
@@ -571,15 +611,50 @@ void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfu
       // With parton showers on we are generating "Z inclusive process"
       SetMSTP(61,1); //Initial QCD & QED showers on
       SetMSTP(71,1); //Final QCD & QED showers on
-      
+      break;
+    case kPyZgamma:
+      //Inclusive production of Z
+      SetMSEL(0);
+      //f fbar -> Z/gamma
+      SetMSUB(1,1);
+      // Initial/final parton shower on (Pythia default)
+      // With parton showers on we are generating "Z inclusive process"
+      SetMSTP(61,1); //Initial QCD & QED showers on
+      SetMSTP(71,1); //Final QCD & QED showers on
+      break;
+      case kPyMBRSingleDiffraction:
+      case kPyMBRDoubleDiffraction:
+      case kPyMBRCentralDiffraction:
       break;  
+      case kPyJetsPWHG:
+      //    N.B.
+      //    ====
+      //    For the case of jet production the following parameter setting
+      //    limits the transverse momentum of secondary scatterings, due
+      //    to multiple parton interactions, to be less than that of the
+      //    primary interaction (see POWHEG Dijet paper arXiv:1012.3380
+      //    [hep-ph] sec. 4.1 and also the PYTHIA Manual).
+      SetMSTP(86,1);
+      //    maximum number of errors before pythia aborts (def=10)
+      SetMSTU(22,10);
+      //    number of warnings printed on the shell
+      SetMSTU(26,20);
+      break;
+      case kPyCharmPWHG:
+      case kPyBeautyPWHG:
+      case kPyWPWHG:
+      //    number of warnings printed on the shell
+      SetMSTU(26,20);
 
     }
 //
 //  Initialize PYTHIA
     SetMSTP(41,1);   // all resonance decays switched on
-    Initialize("CMS","p","p",fEcms);
-    
+    if (process == kPyJetsPWHG || process == kPyCharmPWHG || process == kPyBeautyPWHG || process == kPyWPWHG) {
+      Initialize("USER","","",0.);
+    } else {   
+      Initialize("CMS",fProjectile,fTarget,fEcms);
+    }
 }
 
 Int_t AliPythia6::CheckedLuComp(Int_t kf)
@@ -659,6 +734,28 @@ void  AliPythia6::SetDecayTable()
     for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
 }
 
+void  AliPythia6::Pyjoin(Int_t& npart, Int_t *ipart)
+{
+//  Call Pythia join alogorithm to set up a string between
+//  npart partons, given by indices in array ipart[npart]
+//
+    pyjoin(npart, ipart);
+}
+
+void  AliPythia6::Pyshowq(Int_t ip1, Int_t ip2, Double_t qmax)
+{
+//  Call qPythia showering
+//
+    pyshowq(ip1, ip2, qmax);
+}
+
+void AliPythia6::Qpygin0()
+{
+    //position of the hard scattering in the nuclear overlapping area.
+    //just for qpythia.
+    qpygin0();
+}
+
 void  AliPythia6::Pyclus(Int_t& njet)
 {
 //  Call Pythia clustering algorithm
@@ -704,7 +801,7 @@ void AliPythia6::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECM
 // (2) The nuclear geometry using the Glauber Model
 //     
     
-    fGlauber = new AliFastGlauber();
+    fGlauber = AliFastGlauber::Instance();
     fGlauber->Init(2);
     fGlauber->SetCentralityClass(cMin, cMax); 
 
@@ -749,7 +846,7 @@ void  AliPythia6::Quench()
     Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
     Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
     Bool_t  quenched[4];
-    Double_t wjtKick[4];
+    Double_t wjtKick[4] = {0., 0., 0., 0.};
     Int_t nGluon[4];
     Int_t qPdg[4];
     Int_t   imo, kst, pdg;
@@ -813,7 +910,8 @@ void  AliPythia6::Quench()
 
 
        Int_t itype = (qPdg[j] == 21) ? 2 : 1;
-       Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
+       //      Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
+       Double_t eloss = fQuenchingWeights->GetELossRandomK(itype, int0[j], int1[j], eq[j]);
 
        if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
            fZQuench[j] = 0.;
@@ -1142,7 +1240,7 @@ void  AliPythia6::Quench()
                //
                // Isotropic decay ????
                Double_t cost = 2. * gRandom->Rndm() - 1.;
-               Double_t sint = TMath::Sqrt(1. - cost * cost);
+               Double_t sint = TMath::Sqrt((1.-cost)*(1.+cost));
                Double_t phis =  2. * TMath::Pi() * gRandom->Rndm();
                
                Double_t pz1 =   pst * cost;
@@ -1382,6 +1480,43 @@ void AliPythia6::AtlasTuning()
        SetPARP(67,1);             // Regulates Initial State Radiation
 }
 
+void AliPythia6::AtlasTuningMC09()
+{
+    //
+    // Configuration for the ATLAS tuning
+    printf("ATLAS New TUNE MC09\n");
+    SetMSTP(81,21);             // treatment for MI, ISR, FSR and beam remnants: MI on, new model
+    SetMSTP(82, 4);             // Double Gaussian Model
+    SetMSTP(52, 2);             // External PDF
+    SetMSTP(51, 20650);         // MRST LO*
+  
+    
+    SetMSTP(70, 0);             // (was 2: def manual 1, def code 0) virtuality scale for ISR 
+    SetMSTP(72, 1);             // (was 0: def 1) maximum scale for FSR
+    SetMSTP(88, 1);             // (was 0: def 1) strategy for qq junction to di-quark or baryon in beam remnant
+    SetMSTP(90, 0);             // (was 1: def 0) strategy of compensate the primordial kT
+
+    SetPARP(78, 0.3);           // the amount of color reconnection in the final state
+    SetPARP(80, 0.1);           // probability of color partons kicked out from beam remnant
+    SetPARP(82, 2.3);           // [GeV]    PT_min at Ref. energy    
+    SetPARP(83, 0.8);           // Core density in proton matter distribution (def.value)    
+    SetPARP(84, 0.7);           // Core radius
+    SetPARP(90, 0.25);          //  2*epsilon (exponent in power law)
+    SetPARJ(81, 0.29);          // (was 0.14: def 0.29) Labmda value in running alpha_s for parton showers
+
+    SetMSTP(95, 6);
+    SetPARJ(41, 0.3);           // a and b parameters of the symmm. Lund FF
+    SetPARJ(42, 0.58);
+    SetPARJ(46, 0.75);          // mod. of the Lund FF for heavy end-point quarks
+    SetPARP(89,1800.);         // [GeV]   Ref. energy
+}
+
+void AliPythia6::SetWeightPower(Double_t pow)
+{
+    setpowwght(pow);
+    SetMSTP(142, 1); // Tell Pythia to use pyevwt to calculate event wghts
+}
+
 void AliPythia6::SetPtHardRange(Float_t ptmin, Float_t ptmax)
 {
     // Set the pt hard range