]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliPythia.cxx
lhapdf.f used from the original LHAPDF distribution. Modifications to avoid data...
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
index cd1c2134fa223febaadeaf162d1cc6479a609ada..9173c056a66f04a4b56d8566ac17ec5194b70794 100644 (file)
@@ -1,3 +1,4 @@
+
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
 
 #include "AliPythia.h"
 #include "AliPythiaRndm.h"
-#include "../FASTSIM/AliFastGlauber.h"
-#include "../FASTSIM/AliQuenchingWeights.h"
+#include "AliFastGlauber.h"
+#include "AliQuenchingWeights.h"
 #include "TVector3.h"
+#include "PyquenCommon.h"
 
 ClassImp(AliPythia)
 
@@ -51,7 +53,16 @@ extern "C" void type_of_call pyevnw(){;}
 
 AliPythia* AliPythia::fgAliPythia=NULL;
 
-AliPythia::AliPythia()
+AliPythia::AliPythia():
+    fProcess(kPyMb),
+    fEcms(0.),
+    fStrucFunc(kCTEQ5L),
+    fXJet(0.),
+    fYJet(0.),
+    fNGmax(30),
+    fZmax(0.97),
+    fGlauber(0),
+    fQuenchingWeights(0)
 {
 // Default Constructor
 //
@@ -62,6 +73,23 @@ AliPythia::AliPythia()
     fQuenchingWeights = 0;
 }
 
+AliPythia::AliPythia(const AliPythia& pythia):
+    TPythia6(pythia), 
+    AliRndm(pythia),
+    fProcess(kPyMb),
+    fEcms(0.),
+    fStrucFunc(kCTEQ5L),
+    fXJet(0.),
+    fYJet(0.),
+    fNGmax(30),
+    fZmax(0.97),
+    fGlauber(0),
+    fQuenchingWeights(0)
+{
+    // Copy Constructor
+    pythia.Copy(*this);
+}
+
 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
 {
 // Initialise the process to generate 
@@ -81,9 +109,13 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     SetMDCY(Pycomp(3312),1,0);
     SetMDCY(Pycomp(3322),1,0);
     SetMDCY(Pycomp(3334),1,0);
-    //  select structure function 
+    // Select structure function 
     SetMSTP(52,2);
     SetMSTP(51,strucfunc);
+    // Particles produced in string fragmentation point directly to either of the two endpoints
+    // of the string (depending in the side they were generated from).
+    SetMSTU(16,2);
+
 //
 // Pythia initialisation for selected processes//
 //
@@ -146,24 +178,13 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        SetMSTP(88,2); 
 // (D=1)see can be used to form  baryons (BARYON JUNCTION)
        SetMSTJ(1,1);  
-       SetMSTP(51,kCTEQ5L);// CTEQ 5L        ! CTEQ5L pdf
-       SetMSTP(81,1);      // Multiple Interactions ON
-       SetMSTP(82,4);      // Double Gaussian Model         
-       SetPARP(82,1.8);    // [GeV]    PT_min at Ref. energy
-       SetPARP(89,1000.);  // [GeV]   Ref. energy
-       SetPARP(90,0.16);   // 2*epsilon (exponent in power law)
-       SetPARP(83,0.5);    // Core density in proton matter dist. (def.value)
-       SetPARP(84,0.5);    // Core radius
-       SetPARP(85,0.33);   // Regulates gluon prod. mechanism
-       SetPARP(86,0.66);   // Regulates gluon prod. mechanism
-       SetPARP(67,1);      // Regulate gluon prod. mechanism
+       AtlasTuning();
        break;
     case kPyCharm:
        SetMSEL(4);
 //  heavy quark masses
 
        SetPMAS(4,1,1.2);
-       SetMSTU(16,2);
 //
 //    primordial pT
        SetMSTP(91,1);
@@ -174,7 +195,6 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     case kPyBeauty:
        SetMSEL(5);
        SetPMAS(5,1,4.75);
-       SetMSTU(16,2);
        break;
     case kPyJpsi:
        SetMSEL(0);
@@ -221,22 +241,44 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        SetMSUB(94,1);             // double diffraction
        SetMSUB(95,1);             // low pt production
 
+       AtlasTuning();
+       break;
+    case kPyMbDefault:
+// Minimum Bias pp-Collisions
 //
-// ATLAS Tuning
-//
-       
-        SetMSTP(51, kCTEQ5L);      // CTEQ5L pdf
+//   
+//      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
+
+       break;
+    case kPyLhwgMb:
+// Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
+//  -> Pythia 6.3 or above is needed
+//   
+       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
+
+        SetMSTP(51,kCTEQ6ll);      // CTEQ6ll pdf
+       SetMSTP(52,2);
+       SetMSTP(68,1);
+       SetMSTP(70,2);
        SetMSTP(81,1);             // Multiple Interactions ON
        SetMSTP(82,4);             // Double Gaussian Model
+       SetMSTP(88,1);
 
-       SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
-       SetPARP(89,1000.);         // [GeV]   Ref. energy
-       SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
+       SetPARP(82,2.3);           // [GeV]    PT_min at Ref. energy
        SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
        SetPARP(84,0.5);           // Core radius
-       SetPARP(85,0.33);          // Regulates gluon prod. mechanism
-       SetPARP(86,0.66);          // Regulates gluon prod. mechanism
-       SetPARP(67,1);             // Regulates Initial State Radiation
+       SetPARP(85,0.9);           // Regulates gluon prod. mechanism
+       SetPARP(90,0.2);           // 2*epsilon (exponent in power law)
+
        break;
     case kPyMbNonDiffr:
 // Minimum Bias pp-Collisions
@@ -246,22 +288,20 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        SetMSEL(0);
        SetMSUB(95,1);             // low pt production
 
-//
-// ATLAS Tuning
-//
-       
-       SetMSTP(51,kCTEQ5L);       // CTEQ5L pdf
-       SetMSTP(81,1);             // Multiple Interactions ON
-       SetMSTP(82,4);             // Double Gaussian Model
-
-       SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
-       SetPARP(89,1000.);         // [GeV]   Ref. energy
-       SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
-       SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
-       SetPARP(84,0.5);           // Core radius
-       SetPARP(85,0.33);          // Regulates gluon prod. mechanism
-       SetPARP(86,0.66);          // Regulates gluon prod. mechanism
-       SetPARP(67,1);             // Regulates Initial State Radiation
+       AtlasTuning();
+       break;
+    case kPyMbMSEL1:
+       ConfigHeavyFlavor();
+// Intrinsic <kT^2>
+        SetMSTP(91,1);// Width (1=gaussian) primordial kT dist. inside hadrons
+        SetPARP(91,1.);     // <kT^2> = PARP(91,1.)^2
+        SetPARP(93,5.);     // Upper cut-off
+// Set Q-quark mass
+        SetPMAS(4,1,1.2);   // Charm quark mass
+        SetPMAS(5,1,4.78);  // Beauty quark mass
+       SetPARP(71,4.);     // Defaut value
+// Atlas Tuning
+       AtlasTuning();
        break;
     case kPyJets:
 //
@@ -270,7 +310,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        SetMSEL(1);
  // Pythia Tune A (CDF)
  //
-       SetPARP(67,4.);            // Regulates Initial State Radiation
+       SetPARP(67,2.5);           // Regulates Initial State Radiation (value from best fit to D0 dijet analysis)
        SetMSTP(82,4);             // Double Gaussian Model
        SetPARP(82,2.0);           // [GeV]    PT_min at Ref. energy
        SetPARP(84,0.4);           // Core radius
@@ -285,6 +325,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     case kPyCharmPbPbMNR:
     case kPyD0PbPbMNR:
     case kPyDPlusPbPbMNR:
+    case kPyDPlusStrangePbPbMNR:
       // 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.
@@ -303,6 +344,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     case kPyCharmpPbMNR:
     case kPyD0pPbMNR:
     case kPyDPluspPbMNR:
+    case kPyDPlusStrangepPbMNR:
       // 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.
@@ -322,6 +364,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     case kPyCharmppMNR:
     case kPyD0ppMNR:
     case kPyDPlusppMNR:
+    case kPyDPlusStrangeppMNR:
       // 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.
@@ -338,6 +381,33 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       // Set c-quark mass
        SetPMAS(4,1,1.2);
       break;
+    case kPyCharmppMNRwmi:
+      // 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.
+      // This parameter settings are meant to work with pp collisions
+      // and with kCTEQ5L PDFs.
+      // Added multiple interactions according to ATLAS tune settings.
+      // To get a "reasonable" agreement with MNR results, events have to be 
+      // generated with the minimum ptHard (AliGenPythia::SetPtHard)
+      // set to 2.76 GeV.
+      // To get a "perfect" agreement with MNR results, events have to be 
+      // generated in four ptHard bins with the following relative 
+      // normalizations:
+      // 2.76-3 GeV: 25%
+      //    3-4 GeV: 40%
+      //    4-8 GeV: 29%
+      //     >8 GeV:  6%
+       ConfigHeavyFlavor();
+      // Intrinsic <kT^2>
+       SetMSTP(91,1);
+       SetPARP(91,1.);
+       SetPARP(93,5.);
+
+      // Set c-quark mass
+       SetPMAS(4,1,1.2);
+       AtlasTuning();
+       break;
     case kPyBeautyPbPbMNR:
       // Tuning of Pythia parameters aimed to get a resonable agreement
       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
@@ -397,7 +467,38 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        // Set b-quark mass
        SetPMAS(5,1,4.75);
       break;
+     case kPyBeautyppMNRwmi:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // b-bbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with pp collisions
+      // and with kCTEQ5L PDFs.
+      // Added multiple interactions according to ATLAS tune settings.
+      // To get a "reasonable" agreement with MNR results, events have to be 
+      // generated with the minimum ptHard (AliGenPythia::SetPtHard)
+      // set to 2.76 GeV.
+      // To get a "perfect" agreement with MNR results, events have to be 
+      // generated in four ptHard bins with the following relative 
+      // normalizations:
+      // 2.76-4 GeV:  5% 
+      //    4-6 GeV: 31%
+      //    6-8 GeV: 28%
+      //     >8 GeV: 36%
+        ConfigHeavyFlavor();
+      // QCD scales
+        SetPARP(67,1.0);
+        SetPARP(71,1.0);
+        
+        // Intrinsic <kT>
+        SetMSTP(91,1);
+        SetPARP(91,1.);
+        SetPARP(93,5.);
+
+      // Set b-quark mass
+        SetPMAS(5,1,4.75);
+
+        AtlasTuning();
+        break; 
     case kPyW:
 
       //Inclusive production of W+/-
@@ -450,9 +551,8 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
 //
 //  Initialize PYTHIA
     SetMSTP(41,1);   // all resonance decays switched on
-
     Initialize("CMS","p","p",fEcms);
-
+    
 }
 
 Int_t AliPythia::CheckedLuComp(Int_t kf)
@@ -561,14 +661,13 @@ void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_
 
 
 
-void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod)
+void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax)
 {
 // Initializes 
 // (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
 // (2) The nuclear geometry using the Glauber Model
 //     
-
-
+    
     fGlauber = new AliFastGlauber();
     fGlauber->Init(2);
     fGlauber->SetCentralityClass(cMin, cMax); 
@@ -577,6 +676,9 @@ void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMe
     fQuenchingWeights->InitMult();
     fQuenchingWeights->SetK(k);
     fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
+    fNGmax = ngmax;
+    fZmax  = zmax;
+    
 }
 
 
@@ -694,7 +796,7 @@ void  AliPythia::Quench()
            //
            // Avoid complete loss
            //
-           if (fZQuench[j] == 1.) fZQuench[j] = 0.95;
+           if (fZQuench[j] == 1.) fZQuench[j] = fZmax;
            //
            // Some debug printing
 
@@ -723,7 +825,7 @@ void  AliPythia::Quench()
        if (!quenched[isys]) continue;
        
        nGluon[isys]   = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
-       if (nGluon[isys] > 6) nGluon[isys] = 6;
+       if (nGluon[isys] > fNGmax) nGluon[isys] = fNGmax;
        zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
        wjtKick[isys]  = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
 
@@ -1108,9 +1210,25 @@ void  AliPythia::Quench()
 void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
 {
     // Igor Lokthine's quenching routine
+    // http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt
+
     pyquen(a, ibf, b);
 }
 
+void AliPythia::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl)
+{
+    // Set the parameters for the PYQUEN package.
+    // See comments in PyquenCommon.h
+    
+    
+    PYQPAR.t0    = t0;
+    PYQPAR.tau0  = tau0;
+    PYQPAR.nf    = nf;
+    PYQPAR.iengl = iengl;
+    PYQPAR.iangl = iangl;
+}
+
+
 void AliPythia::Pyevnw()
 {
     // New multiple interaction scenario
@@ -1137,9 +1255,8 @@ void AliPythia::ConfigHeavyFlavor()
     
     // No multiple interactions
     SetMSTP(81,0);
-    SetPARP(81,0.0);
-    SetPARP(82,0.0);
-    
+    SetPARP(81, 0.);
+    SetPARP(82, 0.);    
     // Initial/final parton shower on (Pythia default)
     SetMSTP(61,1);
     SetMSTP(71,1);
@@ -1151,3 +1268,36 @@ void AliPythia::ConfigHeavyFlavor()
     SetMSTP(32,2);
     SetPARP(34,1.0);
 }
+
+void AliPythia::AtlasTuning()
+{
+    //
+    // Configuration for the ATLAS tuning
+        SetMSTP(51, kCTEQ5L);      // CTEQ5L pdf
+       SetMSTP(81,1);             // Multiple Interactions ON
+       SetMSTP(82,4);             // Double Gaussian Model
+       SetPARP(81,1.9);           // Min. pt for multiple interactions (default in 6.2-14) 
+       SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
+       SetPARP(89,1000.);         // [GeV]   Ref. energy
+       SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
+       SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
+       SetPARP(84,0.5);           // Core radius
+       SetPARP(85,0.33);          // Regulates gluon prod. mechanism
+       SetPARP(86,0.66);          // Regulates gluon prod. mechanism
+       SetPARP(67,1);             // Regulates Initial State Radiation
+}
+
+AliPythia& AliPythia::operator=(const  AliPythia& rhs)
+{
+// Assignment operator
+    rhs.Copy(*this);
+    return *this;
+}
+
+ void AliPythia::Copy(TObject&) const
+{
+    //
+    // Copy 
+    //
+    Fatal("Copy","Not implemented!\n");
+}