]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliPythia.cxx
making FillQA consistent with AOD
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
index 6520a1350466c364ae90561aac3f688493b89978..360cd8e59948cae1653b481bd63b2ccdbb5f80f7 100644 (file)
@@ -1,4 +1,3 @@
-
 /**************************************************************************
  * 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 "AliOmegaDalitz.h"
+#include "AliLog.h"
 #include "TVector3.h"
+#include "TLorentzVector.h"
+#include "PyquenCommon.h"
 
 ClassImp(AliPythia)
 
@@ -31,6 +34,10 @@ ClassImp(AliPythia)
 # define pyrobo pyrobo_
 # define pyquen pyquen_
 # define pyevnw pyevnw_
+# define pyshowq pyshowq_
+# define qpygin0 qpygin0_
+# define pytune  pytune_
+# define py2ent  py2ent_
 # define type_of_call
 #else
 # define pyclus PYCLUS
@@ -38,6 +45,10 @@ ClassImp(AliPythia)
 # define pyrobo PYROBO
 # define pyquen PYQUEN
 # define pyevnw PYEVNW
+# define pyshowq PYSHOWQ
+# define qpygin0 QPYGIN0
+# define pytune  PYTUNE
+# define py2ent  PY2ENT
 # define type_of_call _stdcall
 #endif
 
@@ -46,8 +57,11 @@ 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(){;}
-
+extern "C" void type_of_call pyevnw();
+extern "C" void type_of_call pyshowq(Int_t &, Int_t &, Double_t &);
+extern "C" void type_of_call pytune(Int_t &);
+extern "C" void type_of_call py2ent(Int_t &, Int_t&, Int_t&, Double_t&);
+extern "C" void type_of_call qpygin0(); 
 //_____________________________________________________________________________
 
 AliPythia* AliPythia::fgAliPythia=NULL;
@@ -56,10 +70,16 @@ AliPythia::AliPythia():
     fProcess(kPyMb),
     fEcms(0.),
     fStrucFunc(kCTEQ5L),
+    fProjectile("p"),
+    fTarget("p"),
     fXJet(0.),
     fYJet(0.),
+    fNGmax(30),
+    fZmax(0.97),
     fGlauber(0),
-    fQuenchingWeights(0)
+    fQuenchingWeights(0),
+    fItune(-1),
+    fOmegaDalitz()
 {
 // Default Constructor
 //
@@ -68,6 +88,10 @@ AliPythia::AliPythia():
       AliPythiaRndm::SetPythiaRandom(GetRandom());
     fGlauber          = 0;
     fQuenchingWeights = 0;
+    Int_t i = 0;
+    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;
 }
 
 AliPythia::AliPythia(const AliPythia& pythia):
@@ -76,37 +100,48 @@ AliPythia::AliPythia(const AliPythia& pythia):
     fProcess(kPyMb),
     fEcms(0.),
     fStrucFunc(kCTEQ5L),
+    fProjectile("p"),
+    fTarget("p"),
     fXJet(0.),
     fYJet(0.),
+    fNGmax(30),
+    fZmax(0.97),
     fGlauber(0),
-    fQuenchingWeights(0)
+    fQuenchingWeights(0),
+    fItune(-1),
+    fOmegaDalitz()
 {
     // 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 AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
+void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc, Int_t itune)
 {
 // Initialise the process to generate 
     if (!AliPythiaRndm::GetPythiaRandom()) 
       AliPythiaRndm::SetPythiaRandom(GetRandom());
     
+    fItune = itune;
+    
     fProcess = process;
     fEcms = energy;
     fStrucFunc = strucfunc;
 //...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);
+    SetMDCY(Pycomp(111) ,1,0); // pi0
+    SetMDCY(Pycomp(310) ,1,0); // K0S
+    SetMDCY(Pycomp(3122),1,0); // kLambda
+    SetMDCY(Pycomp(3112),1,0); // sigma -
+    SetMDCY(Pycomp(3222),1,0); // sigma +
+    SetMDCY(Pycomp(3312),1,0); // xi - 
+    SetMDCY(Pycomp(3322),1,0); // xi 0
+    SetMDCY(Pycomp(3334),1,0); // omega-
     // Select structure function 
     SetMSTP(52,2);
-    SetMSTP(51,strucfunc);
+    SetMSTP(51, AliStructFuncType::PDFsetIndex(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);
@@ -238,7 +273,8 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
 
        AtlasTuning();
        break;
-    case kPyMbDefault:
+       
+    case kPyMbAtlasTuneMC09:
 // Minimum Bias pp-Collisions
 //
 //   
@@ -249,6 +285,40 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        SetMSUB(94,1);             // double diffraction
        SetMSUB(95,1);             // low pt production
 
+       AtlasTuningMC09();
+       break;
+
+    case kPyMbWithDirectPhoton:
+// Minimum Bias pp-Collisions with direct photon processes added 
+//
+//   
+//      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
+
+       SetMSUB(14,1);             //
+       SetMSUB(18,1);             //
+       SetMSUB(29,1);             //
+       SetMSUB(114,1);            //
+       SetMSUB(115,1);            //
+
+
+       AtlasTuning();
+       break;
+
+    case kPyMbDefault:
+// 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
        break;
     case kPyLhwgMb:
 // Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
@@ -260,7 +330,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        SetMSUB(94,1);             // double diffraction
        SetMSUB(95,1);             // low pt production
 
-        SetMSTP(51,kCTEQ6ll);      // CTEQ6ll pdf
+        SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ6ll));      // CTEQ6ll pdf
        SetMSTP(52,2);
        SetMSTP(68,1);
        SetMSTP(70,2);
@@ -305,15 +375,17 @@ 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
-       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;
+       if (fItune < 0) {
+         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;
@@ -360,6 +432,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     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.
@@ -462,6 +535,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        // 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
@@ -541,14 +615,65 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       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);
+            
+      break;
     }
 //
 //  Initialize PYTHIA
+//
+//  Select the tune
+    if (itune > -1) {
+      Pytune(itune);
+      if (GetMSTP(192) > 1 || GetMSTP(193) > 1) {
+       AliWarning(Form("Structure function for tune %5d set to %5s\n", 
+                       itune,  AliStructFuncType::PDFsetName(strucfunc).Data()));
+       SetMSTP(52,2);
+       SetMSTP(51, AliStructFuncType::PDFsetIndex(strucfunc));
+      }
+    }
+//  
     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);
+    }
+    fOmegaDalitz.Init();
 }
 
 Int_t AliPythia::CheckedLuComp(Int_t kf)
@@ -559,7 +684,7 @@ Int_t AliPythia::CheckedLuComp(Int_t kf)
     return kc;
 }
 
-void AliPythia::SetNuclei(Int_t a1, Int_t a2)
+void AliPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdf)
 {
 // Treat protons as inside nuclei with mass numbers a1 and a2  
 //    The MSTP array in the PYPARS common block is used to enable and 
@@ -570,9 +695,11 @@ void AliPythia::SetNuclei(Int_t a1, Int_t a2)
 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
 //    MSTP(192) : Mass number of nucleus side 1
 //    MSTP(193) : Mass number of nucleus side 2
+//    MSTP(194) : Nuclear structure function: 0: EKS98 8:EPS08 9:EPS09LO 19:EPS09NLO
     SetMSTP(52,2);
     SetMSTP(192, a1);
-    SetMSTP(193, a2);  
+    SetMSTP(193, a2); 
+    SetMSTP(194, pdf);
 }
        
 
@@ -655,16 +782,58 @@ void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_
     pyrobo(imi, ima, the, phi, bex, bey, bez);
 }
 
+void AliPythia::Pytune(Int_t itune)
+{
+/*
+C
+C ITUNE    NAME (detailed descriptions below)
+C     0 Default : No settings changed => linked Pythia version's defaults.
+C ====== Old UE, Q2-ordered showers ==========================================
+C   100       A : Rick Field's CDF Tune A 
+C   101      AW : Rick Field's CDF Tune AW
+C   102      BW : Rick Field's CDF Tune BW
+C   103      DW : Rick Field's CDF Tune DW
+C   104     DWT : Rick Field's CDF Tune DW with slower UE energy scaling
+C   105      QW : Rick Field's CDF Tune QW (NB: needs CTEQ6.1M pdfs externally)
+C   106 ATLAS-DC2: Arthur Moraes' (old) ATLAS tune (ATLAS DC2 / Rome)
+C   107     ACR : Tune A modified with annealing CR
+C   108      D6 : Rick Field's CDF Tune D6 (NB: needs CTEQ6L pdfs externally)
+C   109     D6T : Rick Field's CDF Tune D6T (NB: needs CTEQ6L pdfs externally)
+C ====== Intermediate Models =================================================
+C   200    IM 1 : Intermediate model: new UE, Q2-ordered showers, annealing CR
+C   201     APT : Tune A modified to use pT-ordered final-state showers
+C ====== New UE, interleaved pT-ordered showers, annealing CR ================
+C   300      S0 : Sandhoff-Skands Tune 0 
+C   301      S1 : Sandhoff-Skands Tune 1
+C   302      S2 : Sandhoff-Skands Tune 2
+C   303     S0A : S0 with "Tune A" UE energy scaling
+C   304    NOCR : New UE "best try" without colour reconnections
+C   305     Old : New UE, original (primitive) colour reconnections
+C   306 ATLAS-CSC: Arthur Moraes' (new) ATLAS tune (needs CTEQ6L externally)
+C ======= The Uppsala models =================================================
+C   ( NB! must be run with special modified Pythia 6.215 version )
+C   ( available from http://www.isv.uu.se/thep/MC/scigal/        )
+C   400   GAL 0 : Generalized area-law model. Old parameters
+C   401   SCI 0 : Soft-Colour-Interaction model. Old parameters
+C   402   GAL 1 : Generalized area-law model. Tevatron MB retuned (Skands)
+*/
+    pytune(itune);
+}
 
+void AliPythia::Py2ent(Int_t idx, Int_t pdg1, Int_t pdg2, Double_t p){
+  // Inset 2-parton system at line idx
+  py2ent(idx, pdg1, pdg2, p);
+}
 
-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 = AliFastGlauber::Instance();
     fGlauber->Init(2);
     fGlauber->SetCentralityClass(cMin, cMax); 
 
@@ -672,6 +841,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;
+    
 }
 
 
@@ -706,7 +878,7 @@ void  AliPythia::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;
@@ -789,7 +961,7 @@ void  AliPythia::Quench()
            //
            // Avoid complete loss
            //
-           if (fZQuench[j] == 1.) fZQuench[j] = 0.97;
+           if (fZQuench[j] > fZmax) fZQuench[j] = fZmax;
            //
            // Some debug printing
 
@@ -818,7 +990,7 @@ void  AliPythia::Quench()
        if (!quenched[isys]) continue;
        
        nGluon[isys]   = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
-       if (nGluon[isys] > 30) nGluon[isys] = 30;
+       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]));
 
@@ -954,9 +1126,13 @@ void  AliPythia::Quench()
                    //
                    //     Calculate new px, py
                    //
-                   Double_t pxNew   = jtNew / jt * pxs;
-                   Double_t pyNew   = jtNew / jt * pys;        
+                   Double_t pxNew   = 0;
+                   Double_t pyNew   = 0;
                    
+                   if (jt>0) {
+                     pxNew = jtNew / jt * pxs;
+                     pyNew = jtNew / jt * pys;
+                   }   
 //                 Double_t dpx = pxs - pxNew;
 //                 Double_t dpy = pys - pyNew;
 //                 Double_t dpz = pl  - plNew;
@@ -1094,17 +1270,17 @@ void  AliPythia::Quench()
                //
                // Isotropic decay ????
                Double_t cost = 2. * gRandom->Rndm() - 1.;
-               Double_t sint = TMath::Sqrt(1. - cost * cost);
-               Double_t phi =  2. * TMath::Pi() * gRandom->Rndm();
+               Double_t sint = TMath::Sqrt((1.-cost)*(1.+cost));
+               Double_t phis =  2. * TMath::Pi() * gRandom->Rndm();
                
                Double_t pz1 =   pst * cost;
                Double_t pz2 =  -pst * cost;
                Double_t pt1 =   pst * sint;
                Double_t pt2 =  -pst * sint;
-               Double_t px1 =   pt1 * TMath::Cos(phi);
-               Double_t py1 =   pt1 * TMath::Sin(phi);     
-               Double_t px2 =   pt2 * TMath::Cos(phi);
-               Double_t py2 =   pt2 * TMath::Sin(phi);     
+               Double_t px1 =   pt1 * TMath::Cos(phis);
+               Double_t py1 =   pt1 * TMath::Sin(phis);            
+               Double_t px2 =   pt2 * TMath::Cos(phis);
+               Double_t py2 =   pt2 * TMath::Sin(phis);            
                
                fPyjets->P[0][iGlu] = px1;
                fPyjets->P[1][iGlu] = py1;
@@ -1203,15 +1379,43 @@ 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
     pyevnw();
 }
 
+void  AliPythia::Pyshowq(Int_t ip1, Int_t ip2, Double_t qmax)
+{
+    //  Call medium-modified Pythia jet reconstruction algorithm
+    //
+    pyshowq(ip1, ip2, qmax);
+}
+ void AliPythia::Qpygin0()                                                                      
+ {                                                                                              
+     // New multiple interaction scenario                                                       
+     qpygin0();                                                                                 
+ }
+
 void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
 {
     // Return event specific quenching parameters
@@ -1230,10 +1434,13 @@ void AliPythia::ConfigHeavyFlavor()
     //
     SetMSEL(1);
     
-    // No multiple interactions
-    SetMSTP(81,0);
-    SetPARP(81, 0.);
-    SetPARP(82, 0.);    
+
+    if (fItune < 0) {
+      // 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);
@@ -1250,18 +1457,53 @@ 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
+    if (fItune > -1) return;
+    printf("ATLAS TUNE \n");
+    
+    SetMSTP(51, AliStructFuncType::PDFsetIndex(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
+}
+
+void AliPythia::AtlasTuningMC09()
+{
+    //
+    // Configuration for the ATLAS tuning
+    if (fItune > -1) return;
+    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
 }
 
 AliPythia& AliPythia::operator=(const  AliPythia& rhs)
@@ -1278,3 +1520,31 @@ AliPythia& AliPythia::operator=(const  AliPythia& rhs)
     //
     Fatal("Copy","Not implemented!\n");
 }
+
+void AliPythia::DalitzDecays()
+{
+
+  //
+  // Replace all omega dalitz decays with the correct matrix element decays
+  //
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+    if (fPyjets->K[1][i] != 223)                       continue;       
+    Int_t fd = fPyjets->K[3][i] - 1;
+    Int_t ld = fPyjets->K[4][i] - 1;
+    if (fd < 0)                                        continue;
+    if ((ld - fd) != 2)                                continue;
+    if ((fPyjets->K[1][fd] != 111) ||
+       ((TMath::Abs(fPyjets->K[1][fd+1]) != 11) && (TMath::Abs(fPyjets->K[1][fd+1]) != 13)))
+       continue;
+    TLorentzVector omega(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+    Int_t pdg = TMath::Abs(fPyjets->K[1][fd+1]);
+    fOmegaDalitz.Decay(pdg, &omega);
+    for (Int_t j = 0; j < 3; j++) {
+      for (Int_t k = 0; k < 4; k++) {
+       TLorentzVector vec = (fOmegaDalitz.Products())[2-j];
+       fPyjets->P[k][fd+j] = vec[k];
+      }
+    }
+  }
+}