bugfix: boundery check for static hit array
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
index e1e055a..9173c05 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)
 
@@ -28,11 +30,15 @@ ClassImp(AliPythia)
 # define pycell pycell_
 # define pyshow pyshow_
 # define pyrobo pyrobo_
+# define pyquen pyquen_
+# define pyevnw pyevnw_
 # define type_of_call
 #else
 # define pyclus PYCLUS
 # define pycell PYCELL
 # define pyrobo PYROBO
+# define pyquen PYQUEN
+# define pyevnw PYEVNW
 # define type_of_call _stdcall
 #endif
 
@@ -40,12 +46,23 @@ 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 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(){;}
 
 //_____________________________________________________________________________
 
 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
 //
@@ -56,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 
@@ -65,11 +99,23 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     fProcess = process;
     fEcms = energy;
     fStrucFunc = strucfunc;
-//  don't decay p0
-    SetMDCY(Pycomp(111),1,0);
-//  select structure function 
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+    SetMDCY(Pycomp(111) ,1,0);
+    SetMDCY(Pycomp(310) ,1,0);
+    SetMDCY(Pycomp(3122),1,0);
+    SetMDCY(Pycomp(3112),1,0);
+    SetMDCY(Pycomp(3212),1,0);
+    SetMDCY(Pycomp(3222),1,0);
+    SetMDCY(Pycomp(3312),1,0);
+    SetMDCY(Pycomp(3322),1,0);
+    SetMDCY(Pycomp(3334),1,0);
+    // 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//
 //
@@ -81,13 +127,64 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
 //  select charm production
     switch (process) 
     {
+    case kPyOldUEQ2ordered:  //Old underlying events with Q2 ordered QCD processes
+//        Multiple interactions on.
+       SetMSTP(81,1);
+// Double Gaussian matter distribution.
+       SetMSTP(82,4);
+       SetPARP(83,0.5);
+       SetPARP(84,0.4);
+//  pT0.
+       SetPARP(82,2.0);
+//  Reference energy for pT0 and energy rescaling pace.
+       SetPARP(89,1800);
+       SetPARP(90,0.25);
+//  String drawing almost completely minimizes string length.
+       SetPARP(85,0.9);
+       SetPARP(86,0.95);
+// ISR and FSR activity.
+       SetPARP(67,4);
+       SetPARP(71,4);
+// Lambda_FSR scale.
+       SetPARJ(81,0.29);
+       break;
+    case kPyOldUEQ2ordered2:   
+// Old underlying events with Q2 ordered QCD processes
+// Multiple interactions on.
+       SetMSTP(81,1);
+// Double Gaussian matter distribution.
+       SetMSTP(82,4);
+       SetPARP(83,0.5);
+       SetPARP(84,0.4);
+// pT0.
+       SetPARP(82,2.0);
+// Reference energy for pT0 and energy rescaling pace.
+       SetPARP(89,1800);
+       SetPARP(90,0.16);  // here is the difference with  kPyOldUEQ2ordered
+// String drawing almost completely minimizes string length.
+       SetPARP(85,0.9);
+       SetPARP(86,0.95);
+// ISR and FSR activity.
+       SetPARP(67,4);
+       SetPARP(71,4);
+// Lambda_FSR scale.
+       SetPARJ(81,0.29);       
+       break;
+    case kPyOldPopcorn:  
+// Old production mechanism: Old Popcorn
+       SetMSEL(1);
+       SetMSTJ(12,3); 
+// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
+       SetMSTP(88,2); 
+// (D=1)see can be used to form  baryons (BARYON JUNCTION)
+       SetMSTJ(1,1);  
+       AtlasTuning();
+       break;
     case kPyCharm:
        SetMSEL(4);
-//
 //  heavy quark masses
 
        SetPMAS(4,1,1.2);
-       SetMSTU(16,2);
 //
 //    primordial pT
        SetMSTP(91,1);
@@ -98,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);
@@ -145,21 +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,7);             // 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
@@ -169,34 +288,44 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        SetMSEL(0);
        SetMSUB(95,1);             // low pt production
 
-//
-// ATLAS Tuning
-//
-       
-       SetMSTP(51,7);             // 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:
 //
 //  QCD Jets
 //
        SetMSEL(1);
-       break;
+ // Pythia Tune A (CDF)
+ //
+       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
+       SetPARP(85,0.90) ;         // Regulates gluon prod. mechanism
+       SetPARP(86,0.95);          // Regulates gluon prod. mechanism
+       SetPARP(89,1800.);         // [GeV]   Ref. energy
+       SetPARP(90,0.25);          // 2*epsilon (exponent in power law)
+       break;
     case kPyDirectGamma:
        SetMSEL(10);
        break;
     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.
@@ -204,37 +333,18 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
-
-      // All QCD processes
-      SetMSEL(1);
-
-      // No multiple interactions
-      SetMSTP(81,0);
-      SetPARP(81,0.0);
-      SetPARP(82,0.0);
-
-      // Initial/final parton shower on (Pythia default)
-      SetMSTP(61,1);
-      SetMSTP(71,1);
-
-      // 2nd order alpha_s
-      SetMSTP(2,2);
-
-      // QCD scales
-      SetMSTP(32,2);
-      SetPARP(34,1.0);
-
+       ConfigHeavyFlavor();
       // Intrinsic <kT>
       SetMSTP(91,1);
       SetPARP(91,1.304);
       SetPARP(93,6.52);
-
       // Set c-quark mass
       SetPMAS(4,1,1.2);
-
       break;
     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.
@@ -242,37 +352,19 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
-
-      // All QCD processes
-      SetMSEL(1);
-
-      // No multiple interactions
-      SetMSTP(81,0);
-      SetPARP(81,0.0);
-      SetPARP(82,0.0);
-
-      // Initial/final parton shower on (Pythia default)
-      SetMSTP(61,1);
-      SetMSTP(71,1);
-
-      // 2nd order alpha_s
-      SetMSTP(2,2);
-
-      // QCD scales
-      SetMSTP(32,2);
-      SetPARP(34,1.0);
-
+       ConfigHeavyFlavor();
       // Intrinsic <kT>
-      SetMSTP(91,1);
-      SetPARP(91,1.16);
-      SetPARP(93,5.8);
-
+       SetMSTP(91,1);
+       SetPARP(91,1.16);
+       SetPARP(93,5.8);
+       
       // Set c-quark mass
-      SetPMAS(4,1,1.2);
-
+       SetPMAS(4,1,1.2);
       break;
     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.
@@ -280,35 +372,42 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
-
-      // All QCD processes
-      SetMSEL(1);
-
-      // No multiple interactions
-      SetMSTP(81,0);
-      SetPARP(81,0.0);
-      SetPARP(82,0.0);
-
-      // Initial/final parton shower on (Pythia default)
-      SetMSTP(61,1);
-      SetMSTP(71,1);
-
-      // 2nd order alpha_s
-      SetMSTP(2,2);
-
-      // QCD scales
-      SetMSTP(32,2);
-      SetPARP(34,1.0);
-
+       ConfigHeavyFlavor();
       // Intrinsic <kT^2>
-      SetMSTP(91,1);
-      SetPARP(91,1.);
-      SetPARP(93,5.);
-
+       SetMSTP(91,1);
+       SetPARP(91,1.);
+       SetPARP(93,5.);
+       
       // Set c-quark mass
-      SetPMAS(4,1,1.2);
-
+       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
@@ -317,36 +416,16 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
-
-      // All QCD processes
-      SetMSEL(1);
-
-      // No multiple interactions
-      SetMSTP(81,0);
-      SetPARP(81,0.0);
-      SetPARP(82,0.0);
-
-      // Initial/final parton shower on (Pythia default)
-      SetMSTP(61,1);
-      SetMSTP(71,1);
-
-      // 2nd order alpha_s
-      SetMSTP(2,2);
-
+       ConfigHeavyFlavor();
       // QCD scales
-      SetMSTP(32,2);
-      SetPARP(34,1.0);
-      SetPARP(67,1.0);
-      SetPARP(71,1.0);
-
+       SetPARP(67,1.0);
+       SetPARP(71,1.0);
       // Intrinsic <kT>
-      SetMSTP(91,1);
-      SetPARP(91,2.035);
-      SetPARP(93,10.17);
-
+       SetMSTP(91,1);
+       SetPARP(91,2.035);
+       SetPARP(93,10.17);
       // Set b-quark mass
-      SetPMAS(5,1,4.75);
-
+       SetPMAS(5,1,4.75);
       break;
     case kPyBeautypPbMNR:
       // Tuning of Pythia parameters aimed to get a resonable agreement
@@ -356,36 +435,16 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
-
-      // All QCD processes
-      SetMSEL(1);
-
-      // No multiple interactions
-      SetMSTP(81,0);
-      SetPARP(81,0.0);
-      SetPARP(82,0.0);
-
-      // Initial/final parton shower on (Pythia default)
-      SetMSTP(61,1);
-      SetMSTP(71,1);
-
-      // 2nd order alpha_s
-      SetMSTP(2,2);
-
+       ConfigHeavyFlavor();
       // QCD scales
-      SetMSTP(32,2);
-      SetPARP(34,1.0);
-      SetPARP(67,1.0);
-      SetPARP(71,1.0);
-
+       SetPARP(67,1.0);
+       SetPARP(71,1.0);
       // Intrinsic <kT>
-      SetMSTP(91,1);
-      SetPARP(91,1.60);
-      SetPARP(93,8.00);
-
+       SetMSTP(91,1);
+       SetPARP(91,1.60);
+       SetPARP(93,8.00);
       // Set b-quark mass
-      SetPMAS(5,1,4.75);
-
+       SetPMAS(5,1,4.75);
       break;
     case kPyBeautyppMNR:
       // Tuning of Pythia parameters aimed to get a resonable agreement
@@ -395,44 +454,105 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
-
-      // All QCD processes
-      SetMSEL(1);
-
-      // No multiple interactions
-      SetMSTP(81,0);
-      SetPARP(81,0.0);
-      SetPARP(82,0.0);
-
-      // Initial/final parton shower on (Pythia default)
-      SetMSTP(61,1);
-      SetMSTP(71,1);
-
-      // 2nd order alpha_s
-      SetMSTP(2,2);
-
+       ConfigHeavyFlavor();
       // QCD scales
-      SetMSTP(32,2);
-      SetPARP(34,1.0);
-      SetPARP(67,1.0);
-      SetPARP(71,1.0);
-
-      // Intrinsic <kT>
-      SetMSTP(91,1);
-      SetPARP(91,1.);
-      SetPARP(93,5.);
+       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);
+      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);
+        SetPMAS(5,1,4.75);
+
+        AtlasTuning();
+        break; 
+    case kPyW:
+
+      //Inclusive production of W+/-
+      SetMSEL(0);
+      //f fbar -> W+ 
+      SetMSUB(2,1);
+      //       //f fbar -> g W+
+      //       SetMSUB(16,1);
+      //       //f fbar -> gamma W+
+      //       SetMSUB(20,1);
+      //       //f g -> f W+  
+      //       SetMSUB(31,1);
+      //       //f gamma -> f W+
+      //       SetMSUB(36,1);
+      
+      // Initial/final parton shower on (Pythia default)
+      // With parton showers on we are generating "W inclusive process"
+      SetMSTP(61,1); //Initial QCD & QED showers on
+      SetMSTP(71,1); //Final QCD & QED showers on
+      
+      break;  
+
+    case kPyZ:
+
+      //Inclusive production of Z
+      SetMSEL(0);
+      //f fbar -> Z/gamma
+      SetMSUB(1,1);
+      
+      //       // f fbar -> g Z/gamma
+      //       SetMSUB(15,1);
+      //       // f fbar -> gamma Z/gamma
+      //       SetMSUB(19,1);
+      //       // f g -> f Z/gamma
+      //       SetMSUB(30,1);
+      //       // f gamma -> f Z/gamma
+      //       SetMSUB(35,1);
+      
+      //only Z included, not gamma
+      SetMSTP(43,2);
+      
+      // 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;  
 
-      break;
     }
 //
 //  Initialize PYTHIA
     SetMSTP(41,1);   // all resonance decays switched on
-
     Initialize("CMS","p","p",fEcms);
-
+    
 }
 
 Int_t AliPythia::CheckedLuComp(Int_t kf)
@@ -541,24 +661,23 @@ 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 qTransport, Float_t maxLength, 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); 
 
     fQuenchingWeights = new AliQuenchingWeights();
     fQuenchingWeights->InitMult();
-    fQuenchingWeights->SetQTransport(qTransport);
+    fQuenchingWeights->SetK(k);
     fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
-    fQuenchingWeights->SetLengthMax(Int_t(maxLength));
-    fQuenchingWeights->SampleEnergyLoss();
+    fNGmax = ngmax;
+    fZmax  = zmax;
     
 }
 
@@ -585,95 +704,129 @@ void  AliPythia::Quench()
     static Float_t eMean = 0.;
     static Int_t   icall = 0;
     
-    Double_t p0[2][5];
-    Double_t p1[2][5];
-    Double_t p2[2][5];
-    Int_t   klast[2] = {-1, -1};
+    Double_t p0[4][5];
+    Double_t p1[4][5];
+    Double_t p2[4][5];
+    Int_t   klast[4] = {-1, -1, -1, -1};
 
     Int_t numpart   = fPyjets->N;
-    Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0.;
-    Double_t pxq[2], pyq[2], pzq[2], eq[2], yq[2], mq[2], pq[2], phiq[2], thetaq[2], ptq[2];
-    Bool_t  quenched[2];
-    Double_t phi;
-    Double_t zInitial[2], wjtKick[2];
-    Int_t nGluon[2];
-    
+    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];
+    Int_t nGluon[4];
+    Int_t qPdg[4];
     Int_t   imo, kst, pdg;
+    
 //
-//  Primary partons
+//  Sore information about Primary partons
+//
+//  j =
+//  0, 1 partons from hard scattering
+//  2, 3 partons from initial state radiation
+// 
+    for (Int_t i = 2; i <= 7; i++) {
+       Int_t j = 0;
+       // Skip gluons that participate in hard scattering
+       if (i == 4 || i == 5) continue;
+       // Gluons from hard Scattering
+       if (i == 6 || i == 7) {
+           j = i - 6;
+           pxq[j]    = fPyjets->P[0][i];
+           pyq[j]    = fPyjets->P[1][i];
+           pzq[j]    = fPyjets->P[2][i];
+           eq[j]     = fPyjets->P[3][i];
+           mq[j]     = fPyjets->P[4][i];
+       } else {
+           // Gluons from initial state radiation
+           //
+           // Obtain 4-momentum vector from difference between original parton and parton after gluon 
+           // radiation. Energy is calculated independently because initial state radition does not 
+           // conserve strictly momentum and energy for each partonic system independently.
+           //
+           // Not very clean. Should be improved !
+           //
+           //
+           j = i;
+           pxq[j]    = fPyjets->P[0][i] - fPyjets->P[0][i+2];
+           pyq[j]    = fPyjets->P[1][i] - fPyjets->P[1][i+2];
+           pzq[j]    = fPyjets->P[2][i] - fPyjets->P[2][i+2];
+           mq[j]     = fPyjets->P[4][i];
+           eq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
+       }
+//
+//  Calculate some kinematic variables
 //
-
-    
-    
-    for (Int_t i = 6; i <= 7; i++) {
-       Int_t j = i - 6;
-       
-       pxq[j]    = fPyjets->P[0][i];
-       pyq[j]    = fPyjets->P[1][i];
-       pzq[j]    = fPyjets->P[2][i];
-       eq[j]     = fPyjets->P[3][i];
-       mq[j]     = fPyjets->P[4][i];
        yq[j]     = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
        pq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
        phiq[j]   = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
        ptq[j]    = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
        thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
-       phi       = phiq[j];
-       
-       // Quench only central jets
-       if (TMath::Abs(yq[j]) > 2.5) {
-           zInitial[j] = 0.;
+       qPdg[j]   =  fPyjets->K[1][i];
+    }
+  
+    Double_t int0[4];
+    Double_t int1[4];
+    
+    fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
+
+    for (Int_t j = 0; j < 4; j++) {
+       //
+       // Quench only central jets and with E > 10.
+       //
+
+
+       Int_t itype = (qPdg[j] == 21) ? 2 : 1;
+       Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
+
+       if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
+           fZQuench[j] = 0.;
        } else {
-           pdg =  fPyjets->K[1][i];
-           
-           // Get length in nucleus
-           Double_t l;
-           fGlauber->GetLengthsForPythia(1, &phi, &l, -1.);
-           //
-           // Energy loss for given length and parton typr 
-           Int_t itype = (pdg == 21) ? 2 : 1;
-       
-           Double_t eloss   = fQuenchingWeights->GetELossRandom(itype, l, eq[j]);
-           if (eq[j] > 80. && TMath::Abs(yq[j]) < 0.5) {
+           if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
                icall ++;
                eMean += eloss;
            }
-           
            //
            // Extra pt
-           
-           wjtKick[j] = TMath::Sqrt(l *  fQuenchingWeights->GetQTransport());
+           Double_t l =   fQuenchingWeights->CalcLk(int0[j], int1[j]);     
+           wjtKick[j] = TMath::Sqrt(l *  fQuenchingWeights->CalcQk(int0[j], int1[j]));
            //
            // Fractional energy loss
-           zInitial[j] = eloss / eq[j];
+           fZQuench[j] = eloss / eq[j];
            //
            // Avoid complete loss
            //
-           if (zInitial[j] == 1.) zInitial[j] = 0.95;
+           if (fZQuench[j] == 1.) fZQuench[j] = fZmax;
            //
            // Some debug printing
-           printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n", 
-                  j, itype, eq[j], phi, l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
+
+           
+//         printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n", 
+//                j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
            
-           zInitial[j] = 1.;
-           while (zInitial[j] >= 0.95)  zInitial[j] = gRandom->Exp(0.2);
+//         fZQuench[j] = 0.8;
+//         while (fZQuench[j] >= 0.95)  fZQuench[j] = gRandom->Exp(0.2);
        }
        
-       quenched[j] = (zInitial[j] > 0.01);
+       quenched[j] = (fZQuench[j] > 0.01);
     } // primary partons
-  
+    
+    
+
     Double_t pNew[1000][4];
     Int_t    kNew[1000];
     Int_t icount = 0;
+    Double_t zquench[4];
+    
 //
 //  System Loop    
-    for (Int_t isys = 0; isys < 2; isys++) {
+    for (Int_t isys = 0; isys < 4; isys++) {
 //      Skip to next system if not quenched.
        if (!quenched[isys]) continue;
        
-       nGluon[isys]   = 1 + Int_t(zInitial[isys] / (1. - zInitial[isys]));
-       if (nGluon[isys] > 6) nGluon[isys] = 6;
-       zInitial[isys] = 1. - TMath::Power(1. - zInitial[isys], 1./Double_t(nGluon[isys]));
+       nGluon[isys]   = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
+       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]));
 
 
@@ -686,9 +839,6 @@ void  AliPythia::Quench()
 // Loop on radiation events
 
        for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
-           Double_t zHeavy = zInitial[isys];
-//
-
            while (1) {
                icount = 0;
                for (Int_t k = 0; k < 4; k++)
@@ -709,8 +859,12 @@ void  AliPythia::Quench()
 //      Quarks and gluons only
                    if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
 //      Particles from hard scattering only
+                   
                    if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
-                   if (imo != (isys + 7) && (imo % 1000)  != (isys + 7)) continue;
+                   Int_t imom = imo % 1000;
+                   if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
+                   if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;               
+                   
                    
 //      Skip comment lines
                    if (kst != 1 && kst != 2) continue;
@@ -727,9 +881,8 @@ void  AliPythia::Quench()
                    theta = TMath::ATan2(pt, pz);
                
 //
-//      Save 4-momentum sum for balancing      
-                   Int_t index = imo - 7;
-                   if (index >=  1000) index = imo % 1000 - 7;
+//      Save 4-momentum sum for balancing
+                   Int_t index = isys;
                    
                    p0[index][0] += px;
                    p0[index][1] += py;
@@ -740,7 +893,8 @@ void  AliPythia::Quench()
                    
 //
 //      Fractional energy loss
-                   Double_t z = zInitial[index];
+                   Double_t z = zquench[index];
+                   
                    
 //      Don't fully quench radiated gluons
 //
@@ -748,12 +902,12 @@ void  AliPythia::Quench()
 //      This small factor makes sure that the gluons are not too close in phase space to avoid recombination
 //
 
-                       z = 0.05;
+                       z = 0.02;
                    }
+//                 printf("z: %d %f\n", imo, z);
+                   
 
 //
-
-                   if (m > 0.) z = zHeavy;
                    
                    //
                    //
@@ -854,7 +1008,7 @@ void  AliPythia::Quench()
                    p2[isys][4] = TMath::Sqrt(p2[isys][4]);
                    break;
                } else {
-                   printf("Warning negative mass squared in system %d %f ! \n", isys, zInitial[isys]);
+                   printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
                    printf("4-Momentum: %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[isys][0], p2[isys][1], p2[isys][2], p2[isys][3], p2[isys][4]);
                    if (p2[isys][4] < -0.01) {
                        printf("Negative mass squared !\n");
@@ -1051,3 +1205,99 @@ void  AliPythia::Quench()
     }
 //     this->Pylist(1);
 } // end 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
+    pyevnw();
+}
+
+void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
+{
+    // Return event specific quenching parameters
+    xp = fXJet;
+    yp = fYJet;
+    for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
+
+}
+
+void AliPythia::ConfigHeavyFlavor()
+{
+    //
+    // Default configuration for Heavy Flavor production
+    //
+    // All QCD processes
+    //
+    SetMSEL(1);
+    
+    // No multiple interactions
+    SetMSTP(81,0);
+    SetPARP(81, 0.);
+    SetPARP(82, 0.);    
+    // Initial/final parton shower on (Pythia default)
+    SetMSTP(61,1);
+    SetMSTP(71,1);
+    
+    // 2nd order alpha_s
+    SetMSTP(2,2);
+    
+    // QCD scales
+    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");
+}