updated parametrisation for pi0, eta, omega, phi for different energies, allow select...
authorfbock <friederike.bock@cern.ch>
Tue, 21 Oct 2014 15:31:54 +0000 (17:31 +0200)
committermorsch <andreas.morsch@cern.ch>
Wed, 22 Oct 2014 07:46:43 +0000 (09:46 +0200)
EVGEN/AliGenEMlib.cxx
EVGEN/AliGenEMlib.h

index 46719f5..05d3d2a 100644 (file)
 ClassImp(AliGenEMlib)
 
 //Initializers for static members
-Int_t AliGenEMlib::fgSelectedPtParam=AliGenEMlib::kPizero7TeVpp;
-Int_t AliGenEMlib::fgSelectedCentrality=AliGenEMlib::kpp;
+Int_t AliGenEMlib::fgSelectedCollisionsSystem=AliGenEMlib::kpp7TeV; 
+Int_t AliGenEMlib::fgSelectedPtParamPi0=AliGenEMlib::kPizeroParam; 
+Int_t AliGenEMlib::fgSelectedPtParamEta=AliGenEMlib::kEtaParamRatiopp; 
+Int_t AliGenEMlib::fgSelectedPtParamOmega=AliGenEMlib::kOmegaParampp; 
+Int_t AliGenEMlib::fgSelectedPtParamPhi=AliGenEMlib::kPhiParampp;
+Int_t AliGenEMlib::fgSelectedCentrality=AliGenEMlib::kpp; 
 Int_t AliGenEMlib::fgSelectedV2Systematic=AliGenEMlib::kNoV2Sys;
 
 Double_t AliGenEMlib::CrossOverLc(double a, double b, double x){
-  if(x<b-a/2) return 1.0;
-  else if(x>b+a/2) return 0.0;
-  else return cos(((x-b)/a+0.5)*TMath::Pi())/2+0.5;
+       if(x<b-a/2) return 1.0;
+       else if(x>b+a/2) return 0.0;
+       else return cos(((x-b)/a+0.5)*TMath::Pi())/2+0.5;
 }
 Double_t AliGenEMlib::CrossOverRc(double a, double b, double x){
-  return 1-CrossOverLc(a,b,x);
+       return 1-CrossOverLc(a,b,x);
 }
 
 const Double_t AliGenEMlib::fgkV2param[kCentralities][16] = {
@@ -95,8 +99,8 @@ const Double_t AliGenEMlib::fgkPtParam[kCentralities][10] = {
   ,{ 1.823377e+02, 8.144309e-01, 4.291562e-01, 1.022767e+01, 3.585469e+00, 5.275078e+00, 3.144351e+00, 5.259097e+00, 2.675708e+01, 5.892506e+00 } // 10-20
   ,{ 4.851407e+02, 9.341151e-01, 4.716673e-01, 1.058090e+01, 4.681218e+00, 7.261284e+00, 3.883227e+00, 6.638627e+00, 1.562806e+01, 5.772127e+00 } // 20-30
   ,{ 3.157060e+01, 6.849451e-01, 4.868669e-01, 8.394558e+00, 3.539142e+00, 5.495280e+00, 4.102638e+00, 3.722991e+00, 1.638622e+01, 5.935963e+00 } // 30-40
-  ,{ 1.857919e+01, 6.185989e-01, 5.878869e-01, 7.035064e+00, 2.892415e+00, 4.339383e+00, 3.549679e+00, 2.821061e+00, 1.529318e+01, 6.091388e+00 } // 40-50
-  ,{ 1.069397e+01, 5.816587e-01, 6.542961e-01, 6.472858e+00, 2.643870e+00, 3.929020e+00, 3.339224e+00, 2.410371e+00, 9.606748e+00, 6.116685e+00 } // 50-60
+  ,{ 1.069397e+01, 5.816587e-01, 6.542961e-01, 6.472858e+00, 2.643870e+00, 3.929020e+00, 3.339224e+00, 2.410371e+00, 9.606748e+00, 6.116685e+00 } // 40-50
+  ,{ 1.857919e+01, 6.185989e-01, 5.878869e-01, 7.035064e+00, 2.892415e+00, 4.339383e+00, 3.549679e+00, 2.821061e+00, 1.529318e+01, 6.091388e+00} // 50-60
   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-10 
   ,{ 1.271594e+02, 7.790165e-01, 5.793214e-01, 8.050008e+00, 3.211312e+00, 4.825258e+00, 3.840509e+00, 3.046231e+00, 2.172177e+01, 5.983496e+00 } // 20-40
   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-60
@@ -126,22 +130,100 @@ const Double_t AliGenEMlib::fgkThermPtParam[kCentralities][2] = {
   ,{ 0.0000000000, 0.0000000000 } // 40-80
 };
 
-// MASS   0=>PIZERO, 1=>ETA, 2=>RHO, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI, 6=>JPSI
-const Double_t AliGenEMlib::fgkHM[8] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946, 3.0969, 0.0};
-
-const Double_t AliGenEMlib::fgkMtFactor[2][8] = { 
-  // {1.0, 0.5, 1.0, 0.9, 0.4, 0.23, 0.054},  // factor for pp from arXiv:1110.3929
-  // {1.0, 0.55, 1.0, 0.9, 0.4, 0.25, 0.004}    // factor for PbPb from arXiv:1110.3929
-  //{1., 0.48, 1.0, 0.9, 0.25, 0.4}, (old values)
-  //{1., 0.48, 1.0, 0.9, 0.4, 0.25}, (nlo values)
-  //{1., 0.48, 1.0, 0.8, 0.4, 0.2, 0.06} (combination of nlo and LHC measurements)
-  //https://aliceinfo.cern.ch/Figure/node/2634
-  //https://aliceinfo.cern.ch/Figure/node/2788
-  //https://aliceinfo.cern.ch/Figure/node/4403
-  //https://aliceinfo.cern.ch/Notes/node/87
-  //best guess:
-  {1., 0.48, 1.0, 0.9, 0.4, 0.25, 0., 0.}, //pp
-  {1., 0.48, 1.0, 0.9, 0.4, 0.25, 0., 0.}  //PbPb
+const Double_t AliGenEMlib::fgkModTsallisParamPi0PbPb[kCentralities][7] = { 
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // pp
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 0-5
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 5-10
+       {2.09114,               2.7482,         6.911,          51.1383,        -10.6896,       1.30818,        -1.59137        }, // 10-20
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 20-30
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 30-40
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 40-50
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 50-60
+       {970.684,               0.46451,        5.52064,        21.7707,        -4.2495,        4.62292,        -3.47699        }, // 00-10
+       {3.22534,               2.83717,        7.50505,        38.8015,        -8.93242,       0.9842,         -0.821508       }, // 20-40
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 40-60
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 60-80
+       {44.3972,               1.0644,         5.92254,        40.9254,        -8.07759,       3.30333,        -2.25078        }, // 00-20
+       {38.187,                1.58985,        6.81705,        31.2526,        -8.35251,       3.39482,        -1.56172        }, // 00-40
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 20-80
+       {89.2412,               0.150509,       4.97424,        112.986,        643.257,        3.41969,        -2.46034        }, // 40-80 
+};
+
+const Double_t AliGenEMlib::fgkModTsallisParamPiChargedPbPb[kCentralities][7] = { 
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // pp
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 0-5
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 5-10
+       {117.693,               1.20567,        6.57362,        29.2275,        -7.71065,       4.50046,        -1.91093        }, // 10-20
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 20-30
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 30-40
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 40-50
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 50-60
+       {229.035,               1.11283,        6.53404,        28.9739,        -8.15381,       5.05703,        -2.32825        }, // 00-10
+       {45.5686,               1.39268,        6.72519,        29.4513,        -7.23335,       3.80987,        -1.18599        }, // 20-40
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 40-60
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 60-80
+       {171.812,               1.14849,        6.54742,        29.0473,        -7.96679,       4.82915,        -2.15967        }, // 00-20
+       {108.083,               1.21125,        6.59362,        28.9651,        -7.70072,       4.56583,        -1.89318        }, // 00-40
+       {0.,                    0.,                     0.,                     0.,                     0.,                     0.,                     0.                      }, // 20-80
+       {3.14057,               0.329224,       4.8235,         491.145,        186.041,        2.91138,        -2.20281        }, // 40-80 
+};
+
+const Double_t AliGenEMlib::fgkParamSetPi07TeV[kNPi0Param][7] = { 
+       {0.134977,              2.31335/(2*TMath::Pi()),        0.1433,                 7.003,                  0,                      0,                      0                       }, //kPizeroParam
+       {-0.0580977,    0.580804,                                       5.0964,                 -669.913,               -160.2,         4.45906,        -0.465125       }, //kPizeroParamlow
+       {-7.09454,              0.218554,                                       5.2278,                 77.9877,                -359.457,       4.67862,        -0.996938       }, //kPizeroParamhigh
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPichargedParam
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPichargedParamlow
+       {-0.000632204,  0.371249,                                       4.56778,                -111488,                -15573.4,       4.33064,        -1.5506         }, //kPichargedParamhigh
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPizeroParamAlter
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPizeroParamAlterlow
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPizeroParamAlterhigh
+};
+
+const Double_t AliGenEMlib::fgkParamSetPi02760GeV[kNPi0Param][7] = { 
+       {0.134977,              1.7/(2*TMath::Pi()),            0.135,                  7.1,                    0,                      0,                      0                       }, //kPizeroParam
+       {-1.4583,               0.188108,                                       5.34499,                546.328,                -2142.93,       4.55572,        -1.82455        }, //kPizeroParamlow
+       {-0.0648943,    0.231029,                                       4.39238,                -3705.4,                -35.9761,       4.87127,        -2.00983        }, //kPizeroParamhigh
+       {0.755554,              0.20772,                                        5.24167,                11.8279,                1492.53,        3.93863,        -1.72404        }, //kPichargedParam
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPichargedParamlow
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPichargedParamhigh
+       {0.0610774,             5.86289,                                        -7.83516,               1.29301,                2.62416,        0.,                     0.                      }, //kPizeroParamAlter
+       {0.065338,              5.86705,                                        -9.05494,               1.38435,                3.58848,        0.,                     0.                      }, //kPizeroParamAlterlow
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPizeroParamAlterhigh
+};
+
+const Double_t AliGenEMlib::fgkParamSetPi0900GeV[kNPi0Param][7] = { 
+       {0.134977,              1.5/(2*TMath::Pi()),            0.132,                  7.8,                    0,                      0,                      0                       }, //kPizeroParam
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPizeroParamlow
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPizeroParamhigh
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPichargedParam
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPichargedParamlow
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPichargedParamhigh
+       {0.0467285,             59.6495,                                        -212.865,               0.0584012,              2.83816,        0.,                     0.                      }, //kPizeroParamAlter
+       {0.0468173,             66.9511,                                        -136.287,               0.0298099,              1.17147,        0.,                     0.                      }, //kPizeroParamAlterlow
+       {0.,                    0.,                                                     0.,                             0.,                             0.,                     0.,                     0.                      }, //kPizeroParamAlterhigh
+};
+
+
+// MASS   0=>PIZERO, 1=>ETA, 2=>RHO0, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI, 6=>JPSI, 7=>SIGMA, 8=>K0s, 9=>DELTA++, 10=>DELTA+, 11=>DELTA-, 12=>DELTA0, 13=>Rho+, 14=>Rho-, 15=>K0*
+const Double_t AliGenEMlib::fgkHM[16] = {0.1349766, 0.547853, 0.77549, 0.78265, 0.95778, 1.019455, 3.096916, 1.192642, 0.497614, 1.2311, 1.2349, 1.2349, 1.23340, 0.77549, 0.77549, 0.896};
+
+const Double_t AliGenEMlib::fgkMtFactor[3][16] = { 
+       // {1.0, 0.5, 1.0, 0.9, 0.4, 0.23, 0.054},  // factor for pp from arXiv:1110.3929
+       // {1.0, 0.55, 1.0, 0.9, 0.4, 0.25, 0.004}    // factor for PbPb from arXiv:1110.3929
+       //{1., 0.48, 1.0, 0.9, 0.25, 0.4}, (old values)
+       //{1., 0.48, 1.0, 0.9, 0.4, 0.25}, (nlo values)
+       //{1., 0.48, 1.0, 0.8, 0.4, 0.2, 0.06} (combination of nlo and LHC measurements)
+       //https://aliceinfo.cern.ch/Figure/node/2634
+       //https://aliceinfo.cern.ch/Figure/node/2788
+       //https://aliceinfo.cern.ch/Figure/node/4403
+       //https://aliceinfo.cern.ch/Notes/node/87
+       /*best guess:
+               - pp values for eta/pi0 [arXiv:1205.5724], omega/pi0 [arXiv:1210.5749], phi/(pi+/-) [arXiv:1208.5717] from measured 7 Tev data
+       */
+       {1., 0.476, 1.0, 0.85, 0.4, 0.13, 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0}, //pp
+       {1., 0.476, 1.0, 0.85, 0.4, 0.25, 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0}, //pPb
+       {1., 0.476, 1.0, 0.85, 0.4, 0.25, 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0}  //PbPb
 };
 
 //==========================================================================
@@ -163,11 +245,11 @@ Double_t AliGenEMlib::PtModifiedHagedornThermal(Double_t pt,
                                                 Double_t cT,
                                                 Double_t T)
 {
-  // Modified Hagedorn Thermal fit to Picharged for PbPb:
-  Double_t invYield;
-  invYield = c/TMath::Power(p0+pt/p1,n) + cT*exp(-1.0*pt/T);
+       // Modified Hagedorn Thermal fit to Picharged for PbPb:
+       Double_t invYield;
+       invYield = c/TMath::Power(p0+pt/p1,n) + cT*exp(-1.0*pt/T);
 
-  return invYield*(2*TMath::Pi()*pt);
+       return invYield*(2*TMath::Pi()*pt);
 }
 
 
@@ -179,11 +261,11 @@ Double_t AliGenEMlib::PtModifiedHagedornExp(Double_t pt,
                                            Double_t p0,
                                            Double_t n)
 {
-  // Modified Hagedorn exponentiel fit to Pizero for PbPb:
-  Double_t invYield;
-  invYield = c*TMath::Power(exp(-1*(p1*pt-p2*pt*pt))+pt/p0,-n);
+       // Modified Hagedorn exponentiel fit to Pizero for PbPb:
+       Double_t invYield;
+       invYield = c*TMath::Power(exp(-1*(p1*pt-p2*pt*pt))+pt/p0,-n);
 
-  return invYield*(2*TMath::Pi()*pt);
+       return invYield*(2*TMath::Pi()*pt);
 }
 
 
@@ -196,11 +278,11 @@ Double_t AliGenEMlib::PtModifiedHagedornExp2(Double_t pt,
                                              Double_t d,
                                              Double_t n)
 {
-  // Modified Hagedorn exponential fit to charged pions for pPb:
-  Double_t invYield;
-  invYield = c*TMath::Power(exp(-a*pt-b*pt*pt)+pt/p0+TMath::Power(pt/p1,d),-n);
+       // Modified Hagedorn exponential fit to charged pions for pPb:
+       Double_t invYield;
+       invYield = c*TMath::Power(exp(-a*pt-b*pt*pt)+pt/p0+TMath::Power(pt/p1,d),-n);
 
-  return invYield*(2*TMath::Pi()*pt);
+       return invYield*(2*TMath::Pi()*pt);
 }
 
 Double_t AliGenEMlib::PtTsallis(Double_t pt,
@@ -209,45 +291,115 @@ Double_t AliGenEMlib::PtTsallis(Double_t pt,
                                 Double_t T,
                                 Double_t n)
 {
-  // Tsallis fit to Pizero for pp:
-  Double_t mt;
-  Double_t invYield;
-  mt = sqrt(m*m + pt*pt);
-  invYield = c*((n-1.)*(n-2.))/(n*T*(n*T+m*(n-2.)))*pow(1.+(mt-m)/(n*T),-n);
+       // Tsallis fit to Pizero for pp:
+       Double_t mt;
+       Double_t invYield;
+       
+       mt = sqrt(m*m + pt*pt);
+       invYield = c*((n-1.)*(n-2.))/(n*T*(n*T+m*(n-2.)))*pow(1.+(mt-m)/(n*T),-n);
+
+       return invYield*(2*TMath::Pi()*pt);
+}
+
 
-  return invYield*(2*TMath::Pi()*pt);
+Double_t AliGenEMlib::PtXQCD(  Double_t pt,
+                                Double_t a,
+                                Double_t b,
+                                Double_t c,
+                                Double_t d,
+                                                           Double_t e,
+                                                               Double_t f)
+{
+       // QCD inspired function by Martin Wilde
+       // DISCLAIMER: Please be careful outside of the measured pT range
+       Double_t invYield = 0;
+       if(pt>0.05){
+               invYield = a*pow(pt,-1*(b+c/(pow(pt,d)+e*pow(pt,f))));
+       } else invYield = 100;
+       return invYield*(2*TMath::Pi()*pt);
 }
 
+Double_t AliGenEMlib::PtQCD(   Double_t pt,
+                                Double_t a,
+                                Double_t b,
+                                Double_t c,
+                                Double_t d,
+                                                           Double_t e)
+{
+       // QCD inspired function by Martin Wilde
+       // DISCLAIMER: Please be careful outside of the measured pT range
+       Double_t invYield = 0;
+       if(pt>0.05){
+               invYield = a*pow(pt,-1*(b+c/(pow(pt,d)+e)));
+       } else invYield = 100;
+       return invYield*(2*TMath::Pi()*pt);
+}
+
+Double_t AliGenEMlib::PtModTsallis(    Double_t pt,
+                                Double_t a,
+                                Double_t b,
+                                Double_t c,
+                                Double_t d,
+                                                           Double_t e,
+                                                               Double_t f,
+                                                               Double_t g,
+                                                               Double_t mass)
+{
+
+       Double_t invYield = 0;
+       Double_t mt = sqrt(mass*mass + pt*pt);
+       Double_t pt2 = pt*pt;
+       if(pt>0.05){
+               invYield = a*TMath::Power((1.+(mt-mass)/(b)),-c)*(d+e*pt+pt2)/(f+g*pt+pt2);
+       } else invYield = 100;
+       return invYield*(2*TMath::Pi()*pt);
+}
+
+Double_t AliGenEMlib::PtParticleRatiopp(Double_t pt,
+                                                                               Double_t m1,
+                                                                               Double_t m2,
+                                                                               Double_t c1,
+                                                                               Double_t c2,
+                                                                               Double_t T1,
+                                                                               Double_t T2,
+                                                                               Double_t n)
+{
+       
+       Double_t ratio = 0;
+       if (PtTsallis (pt, m2, c2, T2, n)>0) ratio = PtTsallis (pt, m1, c1, T1, n)/ PtTsallis (pt, m2, c2, T2, n);
+       return ratio;
+}                                                                              
+
+
 // Exponential
 Double_t AliGenEMlib::PtExponential(const Double_t *px, const Double_t *c){
-  const double &pt=px[0];
-  Double_t invYield = c[0]*exp(-pt*c[1]);
-  
-  return invYield*(2*TMath::Pi()*pt);
+       const double &pt=px[0];
+       Double_t invYield = c[0]*exp(-pt*c[1]);
+       
+       return invYield*(2*TMath::Pi()*pt);
 }
 
 // Hagedorn with additional Powerlaw
 Double_t AliGenEMlib::PtModifiedHagedornPowerlaw(const Double_t *px, const Double_t *c){
-  const double &pt=px[0];
-  Double_t invYield = c[0]*pow(c[1]+pt*c[2],-c[3])*CrossOverLc(c[5],c[4],pt)+CrossOverRc(c[7],c[6],pt)*c[8]*pow(pt+0.001,-c[9]); //pt+0.001: prevent powerlaw from exploding for pt->0
-  
-  return invYield*(2*TMath::Pi()*pt+0.001); //pt+0.001: be sure to be > 0
+       const double &pt=px[0];
+       Double_t invYield = c[0]*pow(c[1]+pt*c[2],-c[3])*CrossOverLc(c[5],c[4],pt)+CrossOverRc(c[7],c[6],pt)*c[8]*pow(pt+0.001,-c[9]); //pt+0.001: prevent powerlaw from exploding for pt->0
+       
+       return invYield*(2*TMath::Pi()*pt+0.001); //pt+0.001: be sure to be > 0
 }
 
 // double powerlaw for J/Psi yield
 Double_t AliGenEMlib::PtDoublePowerlaw(const Double_t *px, const Double_t *c){
-  const double &pt=px[0];
-  Double_t yield = c[0]*pt*pow(1+pow(pt*c[1],2),-c[2]);
-  
-  return yield;
+       const double &pt=px[0];
+       Double_t yield = c[0]*pt*pow(1+pow(pt*c[1],2),-c[2]);
+       
+       return yield;
 }
 
 // integral over krollwada with S=1^2*(1-mee^2/mh^2)^3 from mee=0 up to mee=mh
 // approximation is perfect for mh>20MeV
 Double_t AliGenEMlib::IntegratedKrollWada(const Double_t *mh, const Double_t *){
-  if(*mh<0.002941) return 0;
-  return 2*log(*mh/0.000511/exp(1.75))/411.11/TMath::Pi();
+       if(*mh<0.002941) return 0;
+       return 2*log(*mh/0.000511/exp(1.75))/411.11/TMath::Pi();
 }
 
 //--------------------------------------------------------------------------
@@ -257,48 +409,48 @@ Double_t AliGenEMlib::IntegratedKrollWada(const Double_t *mh, const Double_t *){
 //--------------------------------------------------------------------------
 Double_t AliGenEMlib::PtPromptRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  const static Double_t promptGammaPtParam[10] = { 8.715017e-02, 4.439243e-01, 1.011650e+00, 5.193789e+00, 2.194442e+01, 1.062124e+01, 2.469876e+01, 6.052479e-02, 5.611410e-02, 5.169743e+00 };
+       const static Double_t promptGammaPtParam[10] = { 8.715017e-02, 4.439243e-01, 1.011650e+00, 5.193789e+00, 2.194442e+01, 1.062124e+01, 2.469876e+01, 6.052479e-02, 5.611410e-02, 5.169743e+00 };
  
-  return PtModifiedHagedornPowerlaw(px,promptGammaPtParam)*GetTAA(fgSelectedCentrality);
+       return PtModifiedHagedornPowerlaw(px,promptGammaPtParam)*GetTAA(fgSelectedCentrality);
 }
 
 Double_t AliGenEMlib::PtThermalRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return PtExponential(px,fgkThermPtParam[fgSelectedCentrality]);
+       return PtExponential(px,fgkThermPtParam[fgSelectedCentrality]);
 }
 
 Double_t AliGenEMlib::PtDirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return PtPromptRealGamma(px,px)+PtThermalRealGamma(px,px);
+       return PtPromptRealGamma(px,px)+PtThermalRealGamma(px,px);
 }
 
 Int_t AliGenEMlib::IpDirectRealGamma(TRandom *)
 {
-  return 22;
+       return 22;
 }
 
 Double_t AliGenEMlib::YDirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return YFlat(*px);
+       return YFlat(*px);
 }
 
 Double_t AliGenEMlib::V2DirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  const static Double_t v2Param[3][16] = {
-    { 1.004245e-01, 1.057645e+00, 0.000000e+00, 2.836492e+00, 2.819767e+00, -6.231529e-02, 1.173054e+00, 2.836492e+00, 1.881590e-01, 1.183293e-02, 1.252249e+00, 0, 1, 4.876263e-03, 1.518526e+00, 4.5 } // 00-20, based on: https://aliceinfo.cern.ch/Notes/node/249
-    ,{ 1.619000e-01, 1.868201e+00, 6.983303e-15, 2.242170e+00, 4.484339e+00, -1.695734e-02, 2.301359e+00, 2.871469e+00, 1.619000e-01, 2.264320e-02, 1.028641e+00, 0, 1, 8.172203e-03, 1.271637e+00, 4.5 } // 20-40
-    ,{ 1.335000e-01, 1.076916e+00, 1.462605e-08, 2.785732e+00, 5.571464e+00, -2.356156e-02, 2.745437e+00, 2.785732e+00, 1.335000e-01, 1.571589e-02, 1.001131e+00, 0, 1, 5.179715e-03, 1.329344e+00, 4.5 } // 00-40
-  };
-  switch(fgSelectedCentrality){
-  case k0020: return V2Param(px,v2Param[0]); break;
-  case k2040: return V2Param(px,v2Param[1]); break;
-  case k0040: return V2Param(px,v2Param[2]); break;
-    // case k0010: return 0.43*V2Param(px,v2Param[1]); break;  //V2Pizero(0010)/V2Pizero(2040)=0.43 +-0.025
-    // case k1020: return 0.75*V2Param(px,v2Param[1]); break;  //V2Pizero(1020)/V2Pizero(2040)=0.75 +-0.04
-  case k0010: return 0.53*V2Param(px,v2Param[2]); break;  //V2Pizero(0010)/V2Pizero(0040)=0.53 +-0.03
-  case k1020: return 0.91*V2Param(px,v2Param[2]); break;  //V2Pizero(1020)/V2Pizero(0040)=0.91 +-0.04
-  }
-  return 0;
+       const static Double_t v2Param[3][16] = {
+               { 1.004245e-01, 1.057645e+00, 0.000000e+00, 2.836492e+00, 2.819767e+00, -6.231529e-02, 1.173054e+00, 2.836492e+00, 1.881590e-01, 1.183293e-02, 1.252249e+00, 0, 1, 4.876263e-03, 1.518526e+00, 4.5 } // 00-20, based on: https://aliceinfo.cern.ch/Notes/node/249
+               ,{ 1.619000e-01, 1.868201e+00, 6.983303e-15, 2.242170e+00, 4.484339e+00, -1.695734e-02, 2.301359e+00, 2.871469e+00, 1.619000e-01, 2.264320e-02, 1.028641e+00, 0, 1, 8.172203e-03, 1.271637e+00, 4.5 } // 20-40
+               ,{ 1.335000e-01, 1.076916e+00, 1.462605e-08, 2.785732e+00, 5.571464e+00, -2.356156e-02, 2.745437e+00, 2.785732e+00, 1.335000e-01, 1.571589e-02, 1.001131e+00, 0, 1, 5.179715e-03, 1.329344e+00, 4.5 } // 00-40
+       };
+       switch(fgSelectedCentrality){
+               case k0020: return V2Param(px,v2Param[0]); break;
+               case k2040: return V2Param(px,v2Param[1]); break;
+               case k0040: return V2Param(px,v2Param[2]); break;
+               // case k0010: return 0.43*V2Param(px,v2Param[1]); break;  //V2Pizero(0010)/V2Pizero(2040)=0.43 +-0.025
+               // case k1020: return 0.75*V2Param(px,v2Param[1]); break;  //V2Pizero(1020)/V2Pizero(2040)=0.75 +-0.04
+               case k0010: return 0.53*V2Param(px,v2Param[2]); break;  //V2Pizero(0010)/V2Pizero(0040)=0.53 +-0.03
+               case k1020: return 0.91*V2Param(px,v2Param[2]); break;  //V2Pizero(1020)/V2Pizero(0040)=0.91 +-0.04
+       }
+       return 0;
 }
 
 
@@ -314,27 +466,27 @@ Double_t AliGenEMlib::PtPromptVirtGamma( const Double_t *px, const Double_t */*d
 
 Double_t AliGenEMlib::PtThermalVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return IntegratedKrollWada(px,px)*PtThermalRealGamma(px,px);
+       return IntegratedKrollWada(px,px)*PtThermalRealGamma(px,px);
 }
 
 Double_t AliGenEMlib::PtDirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return IntegratedKrollWada(px,px)*PtDirectRealGamma(px,px);
+       return IntegratedKrollWada(px,px)*PtDirectRealGamma(px,px);
 }
 
 Int_t AliGenEMlib::IpDirectVirtGamma(TRandom *)
 {
-  return 220000;
+       return 220000;
 }
 
 Double_t AliGenEMlib::YDirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return YFlat(*px);
+       return YFlat(*px);
 }
 
 Double_t AliGenEMlib::V2DirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return V2DirectRealGamma(px,px);
+       return V2DirectRealGamma(px,px);
 }
 
 //--------------------------------------------------------------------------
@@ -345,278 +497,361 @@ Double_t AliGenEMlib::V2DirectVirtGamma( const Double_t *px, const Double_t */*d
 Int_t AliGenEMlib::IpPizero(TRandom *)
 {
   // Return pizero pdg code
-  return 111;     
+       return 111;     
 }
 
 Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
 {
-  // double pigammacorr=1; //misuse pion for direct gammas, tuned for 0040, iteration 0
-  // pigammacorr*=2.258900e-01*log(*px+0.001)+1.591291e+00;  //iteration 1
-  // pigammacorr*=6.601943e-03*log(*px+0.001)+9.593698e-01;  //iteration 2
-  // pigammacorr*=4.019933e-03*log(*px+0.001)+9.843412e-01;  //iteration 3
-  // pigammacorr*=-4.543991e-03*log(*px+0.001)+1.010886e+00; //iteration 4
-  // return pigammacorr*PtPromptRealGamma(px,px); //now the gammas from the pi->gg decay have the pt spectrum of prompt real gammas
-  
-  // fit functions and corresponding parameter of Pizero pT for pp @ 2.76 TeV and @ 7 TeV and for PbPb @ 2.76 TeV 
-
-  Double_t km=0.;
-  Double_t kc=0.;
-  Double_t kn=0.;
-  Double_t kcT=0.;
-  Double_t kT=0.;
-  Double_t kp0=0.;
-  Double_t kp1=0.;
-  Double_t kp2=0.;
-  Double_t ka=0.;
-  Double_t kb=0.;
-  Double_t kd=0.;
-
-  double n1,n2,n3,n4;
-  int oldCent;
-
-  switch(fgSelectedPtParam|fgSelectedCentrality) {
-    // fit to pi charged, same data like in kPiOldChargedPbPb,
-    // but tested and compared against newest (2014) neutral pi measurement
-  case kPichargedPbPb|k0005:
-  case kPichargedPbPb|k0510:
-  case kPichargedPbPb|k1020:
-  case kPichargedPbPb|k2030:
-  case kPichargedPbPb|k3040:
-  case kPichargedPbPb|k4050:
-  case kPichargedPbPb|k5060:
-  case kPichargedPbPb|k2040:
-    return PtModifiedHagedornPowerlaw(px,fgkPtParam[fgSelectedCentrality]);
-    break;
-  case kPichargedPbPb|k0010:
-    n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0005]);
-    n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0510]);
-    return (n1+n2)/2;
-    break;
-  case kPichargedPbPb|k0020:
-    n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0005]);
-    n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0510]);
-    n3=PtModifiedHagedornPowerlaw(px,fgkPtParam[k1020]);
-    return (n1+n2+2*n3)/4;
-    break;
-  case kPichargedPbPb|k0040:
-    n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0005]);
-    n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0510]);
-    n3=PtModifiedHagedornPowerlaw(px,fgkPtParam[k1020]);
-    n4=PtModifiedHagedornPowerlaw(px,fgkPtParam[k2040]);
-    return (n1+n2+2*n3+4*n4)/8;
-    break;
-  case kPichargedPbPb|k4060:
-    n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k4050]);
-    n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k5060]);
-    return (n1+n2)/2;
-    break;
-
-
-    // fit to pi charged v1
-    // charged pion from ToF, unidentified hadrons scaled with pion from TPC
-    // for Pb-Pb @ 2.76 TeV
-  case kPiOldChargedPbPb|k0005:
-    kc=1347.5; kp0=0.9393; kp1=2.254; kn=11.294; kcT=0.002537; kT=2.414;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k0510:
-    kc=1256.1; kp0=0.9545; kp1=2.248; kn=11.291; kcT=0.002662; kT=2.326;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k2030:
-    kc=7421.6; kp0=1.2059; kp1=1.520; kn=10.220; kcT=0.002150; kT=2.196;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k3040:
-    kc=1183.2; kp0=1.0478; kp1=1.623; kn=9.8073; kcT=0.00198333; kT=2.073;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-    // the following is what went into the Pb-Pb preliminary approval (0-10%)
-  case kPiOldChargedPbPb|k0010:
-    kc=1296.0; kp0=0.968; kp1=2.567; kn=12.27; kcT=0.004219; kT=2.207;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k1020:
-    kc=986.0; kp0=0.9752; kp1=2.376; kn=11.62; kcT=0.003116; kT=2.213;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k2040:
-    kc=17337.0; kp0=1.337; kp1=1.507; kn=10.629; kcT=0.00184; kT=2.234;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k4050:
-    kc=6220.0; kp0=1.322; kp1=1.224; kn=9.378; kcT=0.000595; kT=2.383;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k5060:
-    kc=2319.0; kp0=1.267; kp1=1.188; kn=9.044; kcT=0.000437; kT=2.276;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k4060:
-    kc=4724.0; kp0=1.319; kp1=1.195; kn=9.255; kcT=0.000511; kT=2.344;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k6080:
-    kc=2842.0; kp0=1.465; kp1=0.8324; kn=8.167; kcT=0.0001049; kT=2.29;
-    return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
-    break;
-  case kPiOldChargedPbPb|k0020:
-    oldCent=fgSelectedCentrality;
-    fgSelectedCentrality=k0010;
-    n1=PtPizero(px,px);
-    fgSelectedCentrality=k1020;
-    n2=PtPizero(px,px);
-    fgSelectedCentrality=oldCent;
-    return (n1+n2)/2;
-    break;
-  case kPiOldChargedPbPb|k0040:
-    oldCent=fgSelectedCentrality;
-    fgSelectedCentrality=k0010;
-    n1=PtPizero(px,px);
-    fgSelectedCentrality=k1020;
-    n2=PtPizero(px,px);
-    fgSelectedCentrality=k2040;
-    n3=PtPizero(px,px);
-    fgSelectedCentrality=oldCent;
-    return (n1+n2+2*n3)/4;
-    break;
-
-    // fit to pizero from conversion analysis
-    // for PbPb @ 2.76 TeV
-    // Pi0 spectra --> not final results 
-  case kPizeroPbPb|k0005:
-       kc=1952.832; kp1=0.264; kp2=0.069; kp0=1.206; kn=9.732;
-       return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
-       break;
-  case kPizeroPbPb|k0010:
-       kc=1810.029; kp1=0.291; kp2=0.059; kp0=1.170; kn=9.447;
-       return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
-       break;
-  case kPizeroPbPb|k0020:
-       kc=856.241; kp1=-0.409; kp2=-0.127; kp0=1.219; kn=9.030;
-       return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
-       break;     
-  case kPizeroPbPb|k1020:
-       kc=509.781; kp1=-0.784; kp2=-0.120; kp0=0.931; kn=7.299;
-       return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
-       break;
-  case kPizeroPbPb|k2040:
-       kc=541.049; kp1=0.542; kp2=-0.069; kp0=0.972; kn=7.866;
-       return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
-       break;
-  case kPizeroPbPb|k2080:
-       kc=222.577; kp1=0.634; kp2=0.009; kp0=0.915; kn=7.431;
-       return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
-       break;
-  case kPizeroPbPb|k4080:
-       kc=120.153; kp1=0.7; kp2=-0.14; kp0=0.835; kn=6.980;
-       return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
-       break;
-  case kPizeroPbPb|k0040:
-       kc=560.532; kp1=0.548; kp2=-0.048; kp0=1.055; kn=8.132;
-       return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
-       break;  
-  
-  
-    // fit to charged pions for p-Pb @ 5.02TeV     
-  case kPichargedPPb:
-       kc=235.5; ka=0.6903; kb=0.06864; kp0=2.289; kp1=0.5872; kd=0.6474; kn=7.842; 
-       return PtModifiedHagedornExp2(*px,kc,ka,kb,kp0,kp1,kd,kn);
-       break;
-
-
-    // Tsallis fit to final pizero (PHOS+PCM) -> used for publication
-    // for pp @ 7 TeV
-  case kPizero7TeVpp:
-  case kPizeroEta7TeVpp:
-    km=0.13498; kc=28.01; kT=0.139; kn=6.875;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-  case kPizero7TeVpplow:
-  case kPizeroEta7TeVpplow: 
-    km=0.13498; kc=23.84; kT=0.147; kn=7.025;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-  case kPizero7TeVpphigh:
-  case kPizeroEta7TeVpphigh:
-    km=0.13498; kc=32.47; kT=0.132; kn=6.749;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-    // Tsallis fit to pizero: preliminary result from PCM and PHOS (QM'11)
-    // for pp @ 2.76 TeV
-  case kPizero2760GeVpp:
-  case kPizeroEta2760GeVpp:     
-    km = 0.13498; kc = 19.75; kT = 0.130; kn = 7.051;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-  case kPizero2760GeVpplow:
-  case kPizeroEta2760GeVpplow:
-    km = 0.13498; kc = 16.12; kT = 0.142; kn = 7.327;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-  case kPizero2760GeVpphigh:
-  case kPizeroEta2760GeVpphigh:
-    km = 0.13498; kc = 25.18; kT = 0.118; kn = 6.782;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-
-  default:
-    return NULL;
-  }
+       // double pigammacorr=1; //misuse pion for direct gammas, tuned for 0040, iteration 0
+       // pigammacorr*=2.258900e-01*log(*px+0.001)+1.591291e+00;  //iteration 1
+       // pigammacorr*=6.601943e-03*log(*px+0.001)+9.593698e-01;  //iteration 2
+       // pigammacorr*=4.019933e-03*log(*px+0.001)+9.843412e-01;  //iteration 3
+       // pigammacorr*=-4.543991e-03*log(*px+0.001)+1.010886e+00; //iteration 4
+       // return pigammacorr*PtPromptRealGamma(px,px); //now the gammas from the pi->gg decay have the pt spectrum of prompt real gammas
+       
+       // fit functions and corresponding parameter of Pizero pT for pp @ 2.76 TeV and @ 7 TeV and for PbPb @ 2.76 TeV 
+
+       Double_t km=0.;
+       Double_t kc=0.;
+       Double_t kn=0.;
+       Double_t kcT=0.;
+       Double_t kT=0.;
+       Double_t kp0=0.;
+       Double_t kp1=0.;
+       Double_t kp2=0.;
+       Double_t ka=0.;
+       Double_t kb=0.;
+       Double_t kd=0.;
+
+       switch(fgSelectedCollisionsSystem) {
+               case kPbPb:
+                       switch (fgSelectedPtParamPi0){
+                               case kPichargedParam:
+                                       // fit to pi charged v1
+                                       // charged pion from ToF, unidentified hadrons scaled with pion from TPC
+                                       // for Pb-Pb @ 2.76 TeV
+                                       switch (fgSelectedCentrality){
+                                               case k0005:
+                                                       kc=1347.5; kp0=0.9393; kp1=2.254; kn=11.294; kcT=0.002537; kT=2.414;
+                                                       return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);     
+                                                       break;
+                                               case k0510:
+                                                       break;
+                                                       kc=1256.1; kp0=0.9545; kp1=2.248; kn=11.291; kcT=0.002662; kT=2.326;
+                                                       return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
+                                                       break;
+                                               case k2030:
+                                                       kc=7421.6; kp0=1.2059; kp1=1.520; kn=10.220; kcT=0.002150; kT=2.196;
+                                                       return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
+                                                       break;
+                                               case k3040:
+                                                       kc=1183.2; kp0=1.0478; kp1=1.623; kn=9.8073; kcT=0.00198333; kT=2.073;
+                                                       return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
+                                                       break;
+                                                       // the following is what went into the Pb-Pb preliminary approval (0-10%)
+                                               case k0010:
+                                                       return PtModTsallis(    *px,
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][0],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][1],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][2],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][3],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][4],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][5],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][6],
+                                                                       0.135);
+                                                       break;
+                                               case k1020:
+                                                       return PtModTsallis(    *px,
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][0],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][1],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][2],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][3],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][4],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][5],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][6],
+                                                                       0.135);
+                                                       break;
+                                               case k2040:
+                                                       return PtModTsallis(    *px,
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][0],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][1],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][2],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][3],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][4],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][5],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][6],
+                                                                       0.135);
+                                                       break;
+                                               case k4050:
+                                                       kc=6220.0; kp0=1.322; kp1=1.224; kn=9.378; kcT=0.000595; kT=2.383;
+                                                       return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
+                                                       break;
+                                               case k5060:
+                                                       kc=2319.0; kp0=1.267; kp1=1.188; kn=9.044; kcT=0.000437; kT=2.276;
+                                                       return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
+                                                       break;
+                                               case k4060:
+                                                       kc=4724.0; kp0=1.319; kp1=1.195; kn=9.255; kcT=0.000511; kT=2.344;
+                                                       return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
+                                                       break;
+                                               case k6080:
+                                                       kc=2842.0; kp0=1.465; kp1=0.8324; kn=8.167; kcT=0.0001049; kT=2.29;
+                                                       return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
+                                                       break;
+                                               case k0020:
+                                                       return PtModTsallis(    *px,
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][0],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][1],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][2],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][3],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][4],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][5],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][6],
+                                                                       0.135);
+                                                       break;
+                                               case k0040:
+                                                       return PtModTsallis(    *px,
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][0],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][1],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][2],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][3],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][4],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][5],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][6],
+                                                                       0.135);
+                                                       break;
+                                               case k4080:
+                                                       return PtModTsallis(    *px,
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][0],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][1],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][2],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][3],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][4],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][5],
+                                                                       fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][6],
+                                                                       0.135);
+                                                       break;
+                                               default:
+                                                       return 0;                                                       
+                                       }
+                               case kPizeroParam:      
+                                       return PtModTsallis(    *px,
+                                                                       fgkModTsallisParamPi0PbPb[fgSelectedCentrality][0],
+                                                                       fgkModTsallisParamPi0PbPb[fgSelectedCentrality][1],
+                                                                       fgkModTsallisParamPi0PbPb[fgSelectedCentrality][2],
+                                                                       fgkModTsallisParamPi0PbPb[fgSelectedCentrality][3],
+                                                                       fgkModTsallisParamPi0PbPb[fgSelectedCentrality][4],
+                                                                       fgkModTsallisParamPi0PbPb[fgSelectedCentrality][5],
+                                                                       fgkModTsallisParamPi0PbPb[fgSelectedCentrality][6],
+                                                                       0.135);
+                               default:
+                                       return 0;
+                               
+                       }
+               case kpPb:
+                       // fit to charged pions for p-Pb @ 5.02TeV     
+                       switch (fgSelectedPtParamPi0){
+                               case kPichargedParam:
+                                       kc=235.5; ka=0.6903; kb=0.06864; kp0=2.289; kp1=0.5872; kd=0.6474; kn=7.842; 
+                                       return PtModifiedHagedornExp2(*px,kc,ka,kb,kp0,kp1,kd,kn);
+                                       break;
+                               default:
+                                       return 0;
+       
+                       }
+               case kpp7TeV:
+                       switch (fgSelectedPtParamPi0){
+                               // Tsallis fit to final pizero (PHOS+PCM) -> used for publication
+                               // for pp @ 7 TeV                               
+                               case kPizeroParam: // fit to combined spectrum with stat errors only
+                                       return PtTsallis(*px,fgkParamSetPi07TeV[kPizeroParam][0],fgkParamSetPi07TeV[kPizeroParam][1],fgkParamSetPi07TeV[kPizeroParam][2],fgkParamSetPi07TeV[kPizeroParam][3]);
+                                       break;    
+                               case kPizeroParamlow:
+                                       return PtModTsallis(      *px, 
+                                                                                         fgkParamSetPi07TeV[kPizeroParamlow][0],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamlow][1],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamlow][2],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamlow][3],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamlow][4],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamlow][5],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamlow][6],
+                                                                                         0.135);
+                                       break;
+                               case kPizeroParamhigh:
+                                       return PtModTsallis(      *px, 
+                                                                                         fgkParamSetPi07TeV[kPizeroParamhigh][0],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamhigh][1],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamhigh][2],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamhigh][3],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamhigh][4],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamhigh][5],
+                                                                                         fgkParamSetPi07TeV[kPizeroParamhigh][6],
+                                                                                         0.135);
+                                       break;
+                               case kPichargedParamhigh:       
+                                       return PtModTsallis(      *px, 
+                                                                                         fgkParamSetPi07TeV[kPichargedParamhigh][0],
+                                                                                         fgkParamSetPi07TeV[kPichargedParamhigh][1],
+                                                                                         fgkParamSetPi07TeV[kPichargedParamhigh][2],
+                                                                                         fgkParamSetPi07TeV[kPichargedParamhigh][3],
+                                                                                         fgkParamSetPi07TeV[kPichargedParamhigh][4],
+                                                                                         fgkParamSetPi07TeV[kPichargedParamhigh][5],
+                                                                                         fgkParamSetPi07TeV[kPichargedParamhigh][6],
+                                                                                         0.135);
+                                       break;
+                                       
+                               default:
+                                       return 0;
+                       }
+
+                               
+               case kpp2760GeV:
+                       switch (fgSelectedPtParamPi0){
+                               // Tsallis fit to pizero: published pi0
+                               // for pp @ 2.76 TeV
+                               case kPizeroParam: //published fit parameters
+                                       return PtTsallis(*px,fgkParamSetPi02760GeV[kPizeroParam][0],fgkParamSetPi02760GeV[kPizeroParam][1],fgkParamSetPi02760GeV[kPizeroParam][2],fgkParamSetPi02760GeV[kPizeroParam][3]);
+                                       break;
+                               case kPizeroParamlow:
+                                       return PtModTsallis(    *px, 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamlow][0], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamlow][1], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamlow][2], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamlow][3], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamlow][4], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamlow][5], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamlow][6], 
+                                                                                       0.135);
+                                       break;
+                               case kPizeroParamhigh:
+                                       return PtModTsallis(    *px, 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamhigh][0], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamhigh][1], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamhigh][2], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamhigh][3], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamhigh][4], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamhigh][5], 
+                                                                                       fgkParamSetPi02760GeV[kPizeroParamhigh][6], 
+                                                                                       0.135);
+                                       break;
+                               case kPichargedParam:   
+                                       return PtModTsallis(      *px, 
+                                                                                         fgkParamSetPi02760GeV[kPichargedParam][0],
+                                                                                         fgkParamSetPi02760GeV[kPichargedParam][1],
+                                                                                         fgkParamSetPi02760GeV[kPichargedParam][2],
+                                                                                         fgkParamSetPi02760GeV[kPichargedParam][3],
+                                                                                         fgkParamSetPi02760GeV[kPichargedParam][4],
+                                                                                         fgkParamSetPi02760GeV[kPichargedParam][5],
+                                                                                         fgkParamSetPi02760GeV[kPichargedParam][6],
+                                                                                         0.135);
+                                       break;
+                               case kPizeroParamAlter: 
+                                       return PtQCD(   *px, 
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][0], 
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][1], 
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][2], 
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][3], 
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][4]);
+                                       break;
+                               case kPizeroParamAlterlow:
+                                       return PtQCD(   *px, 
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][0], 
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][1], 
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][2], 
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][3],
+                                                                       fgkParamSetPi02760GeV[kPizeroParamAlter][4]);
+                                       break;
+                               default:
+                                       return 0;                       
+                       }               
+               case kpp900GeV:
+                       switch (fgSelectedPtParamPi0){
+                               // Tsallis fit to pizero: published pi0
+                               // for pp @ 0.9 TeV
+                               case kPizeroParam: //published fit parameters
+                                       return PtTsallis(       *px,
+                                                                               fgkParamSetPi0900GeV[kPizeroParam][0],
+                                                                               fgkParamSetPi0900GeV[kPizeroParam][1],
+                                                                               fgkParamSetPi0900GeV[kPizeroParam][2],
+                                                                               fgkParamSetPi0900GeV[kPizeroParam][3]);
+                                       break;
+                               case kPizeroParamAlter:
+                                       return PtQCD(   *px, 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlter][0], 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlter][1], 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlter][2], 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlter][3], 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlter][4]);
+                                       break;
+                               case kPizeroParamhigh:
+                                       return PtQCD(   *px, 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlterlow][0], 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlterlow][1], 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlterlow][2], 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlterlow][3], 
+                                                                       fgkParamSetPi0900GeV[kPizeroParamAlterlow][4]);
+                                       break;
+                               default:
+                                       return 0;                       
+                       }               
+                       
+               default:
+                       return 0;
+       }
 
 }
 
 Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
 {
-  return YFlat(*py);
+       return YFlat(*py);
 
 }
 
 Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 {
-  double n1,n2,n3,n4,n5;
-  double v1,v2,v3,v4,v5;
-  switch(fgSelectedCentrality) {
-  case k0010:
-    n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
-    v1=V2Param(px,fgkV2param[k0005]);
-    n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
-    v2=V2Param(px,fgkV2param[k0510]);
-    return (n1*v1+n2*v2)/(n1+n2);
-    break;
-  case k0020:
-    n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
-    v1=V2Param(px,fgkV2param[k0005]);
-    n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
-    v2=V2Param(px,fgkV2param[k0510]);
-    n3=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k1020]);
-    v3=V2Param(px,fgkV2param[k1020]);
-    return (n1*v1+n2*v2+2*n3*v3)/(n1+n2+2*n3);
-    break;
-  case k2040:
-    n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k2030]);
-    v1=V2Param(px,fgkV2param[k2030]);
-    n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k3040]);
-    v2=V2Param(px,fgkV2param[k3040]);
-    return (n1*v1+n2*v2)/(n1+n2);
-    break;
-  case k0040:
-    n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
-    v1=V2Param(px,fgkV2param[k0005]);
-    n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
-    v2=V2Param(px,fgkV2param[k0510]);
-    n3=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k1020]);
-    v3=V2Param(px,fgkV2param[k1020]);
-    n4=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k2030]);
-    v4=V2Param(px,fgkV2param[k2030]);
-    n5=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k3040]);
-    v5=V2Param(px,fgkV2param[k3040]);
-    return (n1*v1+n2*v2+2*n3*v3+2*n4*v4+2*n5*v5)/(n1+n2+2*n3+2*n4+2*n5);
-    break;
-
-  default:
-    return V2Param(px,fgkV2param[fgSelectedCentrality]);
-  }
+       double n1,n2,n3,n4,n5;
+       double v1,v2,v3,v4,v5;
+       switch(fgSelectedCentrality) {
+               case k0010:
+                       n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
+                       v1=V2Param(px,fgkV2param[k0005]);
+                       n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
+                       v2=V2Param(px,fgkV2param[k0510]);
+                       return (n1*v1+n2*v2)/(n1+n2);
+                       break;
+               case k0020:
+                       n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
+                       v1=V2Param(px,fgkV2param[k0005]);
+                       n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
+                       v2=V2Param(px,fgkV2param[k0510]);
+                       n3=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k1020]);
+                       v3=V2Param(px,fgkV2param[k1020]);
+                       return (n1*v1+n2*v2+2*n3*v3)/(n1+n2+2*n3);
+                       break;
+               case k2040:
+                       n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k2030]);
+                       v1=V2Param(px,fgkV2param[k2030]);
+                       n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k3040]);
+                       v2=V2Param(px,fgkV2param[k3040]);
+                       return (n1*v1+n2*v2)/(n1+n2);
+                       break;
+               case k0040:
+                       n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
+                       v1=V2Param(px,fgkV2param[k0005]);
+                       n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
+                       v2=V2Param(px,fgkV2param[k0510]);
+                       n3=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k1020]);
+                       v3=V2Param(px,fgkV2param[k1020]);
+                       n4=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k2030]);
+                       v4=V2Param(px,fgkV2param[k2030]);
+                       n5=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k3040]);
+                       v5=V2Param(px,fgkV2param[k3040]);
+                       return (n1*v1+n2*v2+2*n3*v3+2*n4*v4+2*n5*v5)/(n1+n2+2*n3+2*n4+2*n5);
+                       break;
+
+               default:
+                       return V2Param(px,fgkV2param[fgSelectedCentrality]);
+       }
 }
 
 //--------------------------------------------------------------------------
@@ -627,7 +862,7 @@ Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 Int_t AliGenEMlib::IpEta(TRandom *)
 {
   // Return eta pdg code
-  return 221;     
+       return 221;     
 }
 
 Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
@@ -636,56 +871,87 @@ Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
   // fit functions and corresponding parameter of Eta pT for pp @ 2.76 TeV and @ 7 TeV
   // and mtscaled pT 
 
-  Double_t km = 0.;
-  Double_t kc = 0.;
-  Double_t kT = 0.;
-  Double_t kn = 0.;
-
-  switch(fgSelectedPtParam){ 
-    // Tsallis fit to final eta (PHOS+PCM) -> used for final publication
-    // for pp @ 7 TeV
-  case kPizeroEta7TeVpp:
-    km = 0.547853; kc = 2.496; kT = 0.229; kn = 6.985;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-  case kPizeroEta7TeVpplow:
-    km = 0.547853; kc = 1.970; kT = 0.253; kn = 7.591;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-  case kPizeroEta7TeVpphigh:
-    km = 0.547853; kc = 3.060; kT = 0.212; kn = 6.578;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-    // Tsallis fit to preliminary eta (QM'11)
-    // for pp @ 2.76 TeV
-  case kPizeroEta2760GeVpp:
-    km = 0.547853; kc = 1.971; kT = 0.188; kn = 6.308;
-    return PtTsallis(*px,km,kc,kT,kn);
-  case kPizeroEta2760GeVpplow:
-    km = 0.547853; kc = 1.228; kT = 0.220; kn = 7.030;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-  case kPizeroEta2760GeVpphigh:
-    km = 0.547853; kc = 2.802; kT = 0.164; kn = 5.815;
-    return PtTsallis(*px,km,kc,kT,kn);
-    break;
-
-  default:
-  return MtScal(*px,1);
-    break;
-
-  }
-
+       // parameters for Tsallis fit to eta
+       Double_t km = 0.;
+       Double_t kc = 0.;
+       Double_t kT = 0.;
+       Double_t kn = 0.;
+
+       // parameters for Tsallis fit to pi0
+       Double_t kmPi0 = 0.;
+       Double_t kcPi0 = 0.;
+       Double_t kTPi0 = 0.;
+       Double_t knPi0 = 0.;
+
+       // parameters for fit to eta/pi0
+       Double_t krm1 = 0.;
+       Double_t krm2 = 0.;
+       Double_t krc1 = 0.;
+       Double_t krc2 = 0.;
+       Double_t krT1 = 0.;
+       Double_t krT2 = 0.;
+       Double_t krn = 0.;
+       
+       switch(fgSelectedCollisionsSystem){
+               case kpp7TeV:
+                       switch(fgSelectedPtParamEta){ 
+                               // Tsallis fit to final eta (PHOS+PCM) -> used stat errors only for final publication
+                               // for pp @ 7 TeV
+                               case kEtaParamRatiopp:
+                                       krm1 = 0.547853; krm2 = 0.134977; krc1 = 1.44198e+11; krc2 = 2.06751e+12 ; krT1 = 0.154567 ; krT2 = 0.139634 ; krn=32.0715; 
+                                       kmPi0=0.134977; kcPi0=2.31335/(2*TMath::Pi()); kTPi0=0.1433; knPi0=7.003;
+                                       return PtParticleRatiopp(*px, krm1, krm2, krc1, krc2, krT1, krT2, krn) * PtTsallis(*px,kmPi0,kcPi0,kTPi0,knPi0);
+                                       break;
+                               case kEtaParampp:
+                                       km = 0.547853; kc = 0.290164/(2*TMath::Pi()); kT = 0.212; kn = 7.352;
+                                       return PtTsallis(*px,km,kc,kT,kn);
+                                       break;
+                               // NOTE: None of these parametrisations look right - no idea where they come from       
+                               case kEtaParampplow:
+                                       km = 0.547853; kc = 1.970; kT = 0.253; kn = 7.591;
+                                       return PtTsallis(*px,km,kc,kT,kn);
+                                       break;
+                               case kEtaParampphigh:
+                                       km = 0.547853; kc = 3.060; kT = 0.212; kn = 6.578;
+                                       return PtTsallis(*px,km,kc,kT,kn);
+                                       break;
+                               default:
+                               return MtScal(*px,kEta);
+                       }
+               case kpp2760GeV:        
+                       switch(fgSelectedPtParamEta){ 
+                               // Tsallis fit to preliminary eta (QM'11)
+                               // for pp @ 2.76 TeV
+                               // NOTE: None of these parametrisations look right - no idea where they come from
+                               case kEtaParampp:
+                                       km = 0.547853; kc = 1.971; kT = 0.188; kn = 6.308;
+                                       return PtTsallis(*px,km,kc,kT,kn);
+                               case kEtaParampplow:
+                                       km = 0.547853; kc = 1.228; kT = 0.220; kn = 7.030;
+                                       return PtTsallis(*px,km,kc,kT,kn);
+                                       break;
+                               case kEtaParampphigh:
+                                       km = 0.547853; kc = 2.802; kT = 0.164; kn = 5.815;
+                                       return PtTsallis(*px,km,kc,kT,kn);
+                                       break;
+                               default:
+                                       return MtScal(*px,kEta);
+                                       break;
+                       }
+               default:
+                       return MtScal(*px,kEta);
+                       break;  
+       }
 }
 
 Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
 {
-  return YFlat(*py);
+       return YFlat(*py);
 }
 
 Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return KEtScal(*px,1); //V2Param(px,fgkV2param[1][fgSelectedV2Param]);
+       return KEtScal(*px,kEta); //V2Param(px,fgkV2param[1][fgSelectedV2Param]);
 }
 
 //--------------------------------------------------------------------------
@@ -693,28 +959,29 @@ Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
 //                              Rho
 //
 //--------------------------------------------------------------------------
-Int_t AliGenEMlib::IpRho(TRandom *)
+Int_t AliGenEMlib::IpRho0(TRandom *)
 {
-  // Return rho pdg code
-  return 113;     
+       // Return rho pdg code
+       return 113;     
 }
 
-Double_t AliGenEMlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::PtRho0( const Double_t *px, const Double_t */*dummy*/ )
 {
-  // Rho pT
-  return MtScal(*px,2);
+       // Rho pT
+       return MtScal(*px,kRho0);
 }
 
-Double_t AliGenEMlib::YRho( const Double_t *py, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::YRho0( const Double_t *py, const Double_t */*dummy*/ )
 {
-  return YFlat(*py);
+       return YFlat(*py);
 }
 
-Double_t AliGenEMlib::V2Rho( const Double_t *px, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::V2Rho0( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return KEtScal(*px,2);
+       return KEtScal(*px,kRho0);
 }
 
+
 //--------------------------------------------------------------------------
 //
 //                              Omega
@@ -722,24 +989,72 @@ Double_t AliGenEMlib::V2Rho( const Double_t *px, const Double_t */*dummy*/ )
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpOmega(TRandom *)
 {
-  // Return omega pdg code
-  return 223;     
+       // Return omega pdg code
+       return 223;     
 }
 
 Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
 {
-  // Omega pT
-  return MtScal(*px,3);
+       // Omega pT
+       // fit functions and corresponding parameter of Omega pT for preliminary pp @ 7 TeV data
+       // and mtscaled pT 
+
+       // parameters for Tsallis fit to omega
+       Double_t km = 0.;
+       Double_t kc = 0.;
+       Double_t kT = 0.;
+       Double_t kn = 0.;
+
+       // parameters for Tsallis fit to pi0
+       Double_t kmPi0 = 0.;
+       Double_t kcPi0 = 0.;
+       Double_t kTPi0 = 0.;
+       Double_t knPi0 = 0.;
+
+       // parameters for fit to omega/pi0
+       Double_t krm1 = 0.;
+       Double_t krm2 = 0.;
+       Double_t krc1 = 0.;
+       Double_t krc2 = 0.;
+       Double_t krT1 = 0.;
+       Double_t krT2 = 0.;
+       Double_t krn = 0.;
+       
+       switch(fgSelectedCollisionsSystem){
+               case kpp7TeV:
+                       switch(fgSelectedPtParamOmega){ 
+                               // Tsallis fit to final omega (PHOS) -> stat errors only, preliminary QM12
+                               // for pp @ 7 TeV
+                               case kOmegaParamRatiopp:
+                                       krm1 = 0.78265; krm2 = 0.134977; krc1 = 21240028553.4600143433; krc2 = 168266377865.0805969238 ; krT1 = 0.21175 ; krT2 = 0.14328 ; krn=12.8831831756; 
+                                       kmPi0=0.134977; kcPi0=2.31335/(2*TMath::Pi()); kTPi0=0.1433; knPi0=7.003;
+                                       return PtParticleRatiopp(*px, krm1, krm2, krc1, krc2, krT1, krT2, krn) * PtTsallis(*px,kmPi0,kcPi0,kTPi0,knPi0);
+                                       break;
+                               case kOmegaParampp:
+                                       km = 0.78265; kc = 0.340051/(2*TMath::Pi()); kT = 0.206; kn = 6.31422;
+                                       return PtTsallis(*px,km,kc,kT,kn);
+                                       break;
+                               default:
+                               return MtScal(*px,kOmega);
+                       }
+               default:
+                       return MtScal(*px,kOmega);
+                       break;  
+       }
+       
+       return MtScal(*px,kOmega);
+
 }
 
 Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
 {
-  return YFlat(*py);
+       return YFlat(*py);
 }
 
 Double_t AliGenEMlib::V2Omega( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return KEtScal(*px,3);
+       return KEtScal(*px,kOmega);
+
 }
 
 
@@ -750,25 +1065,25 @@ Double_t AliGenEMlib::V2Omega( const Double_t *px, const Double_t */*dummy*/ )
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpEtaprime(TRandom *)
 {
-  // Return etaprime pdg code
-  return 331;     
+       // Return etaprime pdg code
+       return 331;     
 }
 
 Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
 {
-  // Eta pT
-  return MtScal(*px,4);
+       // Eta pT
+       return MtScal(*px,kEtaprime);
 }
 
 Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
 {
-  return YFlat(*py);
+       return YFlat(*py);
 
 }
 
 Double_t AliGenEMlib::V2Etaprime( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return KEtScal(*px,4);
+       return KEtScal(*px,kEtaprime);
 }
 
 //--------------------------------------------------------------------------
@@ -778,24 +1093,67 @@ Double_t AliGenEMlib::V2Etaprime( const Double_t *px, const Double_t */*dummy*/
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpPhi(TRandom *)
 {
-  // Return phi pdg code
-  return 333;     
+       // Return phi pdg code
+       return 333;     
 }
 
 Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
 {
   // Phi pT
-  return MtScal(*px,5);
+  // fit functions and corresponding parameter of Phi pT for preliminary pp @ 7 TeV data
+       // and PbPb collisions
+       // and mtscaled pT 
+
+       // parameters for Tsallis fit to phi
+       Double_t km = 0.;
+       Double_t kc = 0.;
+       Double_t kT = 0.;
+       Double_t kn = 0.;
+
+       
+       switch(fgSelectedCollisionsSystem){
+//             case kPbPb:
+//                     switch(fgSelectedCentrality){
+//                             // Tsallis fit to final phi->K+K- (TPC, ITS) -> stat+syst
+//                             case k0010:
+//                                     switch(fgSelectedPtParamPhi){ 
+//                                             case kPhiParamPbPb:
+//                                                     km = 0.78265; kc = 0.340051/(2*TMath::Pi()); kT = 0.206; kn = 6.31422;
+//                                                     return PtTsallis(*px,km,kc,kT,kn);
+//                                                     break;
+//                                             default:
+//                                                     return MtScal(*px,kPhi);
+//                                     }       
+//                             default:
+//                                     return MtScal(*px,kPhi);
+//                     }
+               case kpp7TeV:
+                       switch(fgSelectedPtParamPhi){ 
+                               // Tsallis fit to final phi->K+K- (TPC, ITS) -> stat+syst
+                               // for pp @ 7 TeV
+                               case kPhiParampp:
+                                       km = 1.01946; kc = 0.0269578/(2*TMath::Pi()); kT = 0.2718119311; kn = 6.6755739295;
+                                       return PtTsallis(*px,km,kc,kT,kn);
+                                       break;
+                               default:
+                               return MtScal(*px,kPhi);
+                       }
+               default:
+                       return MtScal(*px,kPhi);
+                       break;  
+       }
+       return MtScal(*px,kPhi);
+
 }
 
 Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
 {
-  return YFlat(*py);
+       return YFlat(*py);
 }
 
 Double_t AliGenEMlib::V2Phi( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return KEtScal(*px,5);
+       return KEtScal(*px,kPhi);
 }
 
 //--------------------------------------------------------------------------
@@ -805,59 +1163,363 @@ Double_t AliGenEMlib::V2Phi( const Double_t *px, const Double_t */*dummy*/ )
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpJpsi(TRandom *)
 {
-  // Return phi pdg code
-  return 443;
+       // Return phi pdg code
+       return 443;
 }
 
 Double_t AliGenEMlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/ )
 {
-  // Jpsi pT
-  // based on: //https://aliceinfo.cern.ch/Notes/node/242, https://aliceinfo.cern.ch/Figure/node/3457, www.sciencedirect.com/science/article/pii/S0370269312011446
-  const static Double_t jpsiPtParam[2][3] = {
-    {  9.686337e-03, 2.629441e-01, 4.552044e+00 }
-    ,{ 3.403549e-03, 2.897061e-01, 3.644278e+00 }
-  };
-  const double pt=px[0]*2.28/2.613;
-  switch(fgSelectedCentrality) {
-  case k0020: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[0]); break;
-  case k2040: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[1]); break;
-  case k0040: return 0.5*2.405*(PtDoublePowerlaw(&pt,jpsiPtParam[0])+PtDoublePowerlaw(&pt,jpsiPtParam[1])); break;
-  }
-  return 0;
+       // Jpsi pT
+       // based on: //https://aliceinfo.cern.ch/Notes/node/242, https://aliceinfo.cern.ch/Figure/node/3457, www.sciencedirect.com/science/article/pii/S0370269312011446
+       const static Double_t jpsiPtParam[2][3] = {
+               {  9.686337e-03, 2.629441e-01, 4.552044e+00 }
+               ,{ 3.403549e-03, 2.897061e-01, 3.644278e+00 }
+       };
+       const double pt=px[0]*2.28/2.613;
+       switch(fgSelectedCentrality) {
+               case k0020: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[0]); break;
+               case k2040: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[1]); break;
+               case k0040: return 0.5*2.405*(PtDoublePowerlaw(&pt,jpsiPtParam[0])+PtDoublePowerlaw(&pt,jpsiPtParam[1])); break;
+               default:
+                       return MtScal(*px,kJpsi);
+
+       }
+       return 0;
 }
 
 Double_t AliGenEMlib::YJpsi( const Double_t *py, const Double_t */*dummy*/ )
 {
-  return YFlat(*py);
+       return YFlat(*py);
 }
 
 Double_t AliGenEMlib::V2Jpsi( const Double_t *px, const Double_t */*dummy*/ )
 {
-  const static Double_t v2Param[16] = { 1.156000e-01, 8.936854e-01, 0.000000e+00, 4.000000e+00, 6.222375e+00, -1.600314e-01, 8.766676e-01, 7.824143e+00, 1.156000e-01, 3.484503e-02, 4.413685e-01, 0, 1, 3.484503e-02, 4.413685e-01, 7.2 };
-  switch(fgSelectedCentrality){
-  case k2040: return V2Param(px,v2Param); break;
-  case k0010: return 0.43*V2Param(px,v2Param); break;  //V2Pizero(0010)/V2Pizero(2040)=0.43 +-0.025
-  case k1020: return 0.75*V2Param(px,v2Param); break;  //V2Pizero(1020)/V2Pizero(2040)=0.75 +-0.04
-  case k0020: return 0.66*V2Param(px,v2Param); break;  //V2Pizero(0020)/V2Pizero(2040)=0.66 +-0.035
-  case k0040: return 0.82*V2Param(px,v2Param); break;  //V2Pizero(0040)/V2Pizero(2040)=0.82 +-0.05
-  }
-  return 0;
+       const static Double_t v2Param[16] = { 1.156000e-01, 8.936854e-01, 0.000000e+00, 4.000000e+00, 6.222375e+00, -1.600314e-01, 8.766676e-01, 7.824143e+00, 1.156000e-01, 3.484503e-02, 4.413685e-01, 0, 1, 3.484503e-02, 4.413685e-01, 7.2 };
+       switch(fgSelectedCentrality){
+               case k2040: return V2Param(px,v2Param); break;
+               case k0010: return 0.43*V2Param(px,v2Param); break;  //V2Pizero(0010)/V2Pizero(2040)=0.43 +-0.025
+               case k1020: return 0.75*V2Param(px,v2Param); break;  //V2Pizero(1020)/V2Pizero(2040)=0.75 +-0.04
+               case k0020: return 0.66*V2Param(px,v2Param); break;  //V2Pizero(0020)/V2Pizero(2040)=0.66 +-0.035
+               case k0040: return 0.82*V2Param(px,v2Param); break;  //V2Pizero(0040)/V2Pizero(2040)=0.82 +-0.05
+       }
+       return 0;
+}
+
+//--------------------------------------------------------------------------
+//
+//                              Sigma
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpSigma(TRandom *)
+{
+       // Return Sigma pdg code
+       return 3212;     
+}
+
+Double_t AliGenEMlib::PtSigma( const Double_t *px, const Double_t */*dummy*/ )
+{
+       // Sigma pT
+       return MtScal(*px,kSigma0);
+}
+
+Double_t AliGenEMlib::YSigma( const Double_t *py, const Double_t */*dummy*/ )
+{
+       return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::V2Sigma0( const Double_t *px, const Double_t */*dummy*/ )
+{
+       return KEtScal(*px,kSigma0);
+}
+
+
+//--------------------------------------------------------------------------
+//
+//                              K0short
+//
+//--------------------------------------------------------------------------
+
+Int_t AliGenEMlib::IpK0short(TRandom *)
+{
+       // Return kzeroshort pdg code
+       return 310;
+}
+
+Double_t AliGenEMlib::PtK0short( const Double_t *px, const Double_t */*dummy*/ )
+{
+       // K0short pT
+
+       Double_t ka = 0; 
+       Double_t kb = 0; 
+       Double_t kc = 0; 
+       Double_t kd = 0; 
+       Double_t ke = 0; 
+       Double_t kf = 0;
+               
+       switch (fgSelectedCentrality){
+               case k0010:
+                       ka =9.21859; kb=5.71299; kc=-3.34251; kd=0.48796; ke=0.0192272; kf=3.82224;
+                       return PtXQCD(  *px, ka, kb, kc, kd, ke, kf);
+                       break;
+               case k1020:
+                       ka=6.2377; kb=5.6133; kc=-117.295; kd=3.51154; ke=36.3047; kf=0.456243;
+                       return PtXQCD(  *px, ka, kb, kc, kd, ke, kf);
+                       break;
+               case k0020:
+                       ka=7.7278; kb=5.6686; kc=-3.29259; kd=0.475403; ke=0.0223951; kf=3.69326;
+                       return PtXQCD(  *px, ka, kb, kc, kd, ke, kf);
+                       break;
+               case k2040:
+                       ka=3.38301; kb= 5.5323; kc=-96.078; kd=3.30782; ke=31.085; kf=0.466908;
+                       return PtXQCD(  *px, ka, kb, kc, kd, ke, kf);
+                       break;
+               case k0040:
+                       ka=5.55478; kb=5.61919; kc=-125.635; kd=3.5637; ke=38.9668; kf=0.47068;
+                       return PtXQCD(  *px, ka, kb, kc, kd, ke, kf);
+                       break;
+               case k4080:
+                       ka=0.731606; kb=5.49931; kc=-25.3106; kd=2.2439; ke=8.25063; kf= 0.289288;
+                       return PtXQCD(  *px, ka, kb, kc, kd, ke, kf);
+                       break;
+               default:
+                       return MtScal(*px,kK0s);
+                       break;
+                                       
+       }
+}
+Double_t AliGenEMlib::YK0short( const Double_t *py, const Double_t */*dummy*/ )
+{
+       return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::V2K0sshort( const Double_t *px, const Double_t */*dummy*/ )
+{
+       return KEtScal(*px,kK0s);
+}
+
+
+//--------------------------------------------------------------------------
+//
+//                              Delta ++
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpDeltaPlPl(TRandom *)
+{
+       // Return Delta++ pdg code
+       return 2224;     
+}
+
+Double_t AliGenEMlib::PtDeltaPlPl( const Double_t *px, const Double_t */*dummy*/ )
+{
+       // Delta++ pT
+       return MtScal(*px,kDeltaPlPl);
+}
+
+Double_t AliGenEMlib::YDeltaPlPl( const Double_t *py, const Double_t */*dummy*/ )
+{
+       return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::V2DeltaPlPl( const Double_t *px, const Double_t */*dummy*/ )
+{
+       return KEtScal(*px,kDeltaPlPl);
+}
+
+
+//--------------------------------------------------------------------------
+//
+//                              Delta +
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpDeltaPl(TRandom *)
+{
+       // Return Delta+ pdg code
+       return 2214;     
+}
+
+Double_t AliGenEMlib::PtDeltaPl( const Double_t *px, const Double_t */*dummy*/ )
+{
+       // Delta+ pT
+       return MtScal(*px,kDeltaPl);
+}
+
+Double_t AliGenEMlib::YDeltaPl( const Double_t *py, const Double_t */*dummy*/ )
+{
+       return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::V2DeltaPl( const Double_t *px, const Double_t */*dummy*/ )
+{
+       return KEtScal(*px,kDeltaPl);
+}
+
+
+//--------------------------------------------------------------------------
+//
+//                              Delta -
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpDeltaMi(TRandom *)
+{
+       // Return Delta- pdg code
+       return 1114;     
+}
+
+Double_t AliGenEMlib::PtDeltaMi( const Double_t *px, const Double_t */*dummy*/ )
+{
+       // Delta- pT
+       return MtScal(*px,kDeltaMi);
+}
+
+Double_t AliGenEMlib::YDeltaMi( const Double_t *py, const Double_t */*dummy*/ )
+{
+       return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::V2DeltaMi( const Double_t *px, const Double_t */*dummy*/ )
+{
+       return KEtScal(*px,kDeltaMi);
+}
+
+
+
+//--------------------------------------------------------------------------
+//
+//                              Delta 0
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpDeltaZero(TRandom *)
+{
+       // Return Delta0 pdg code
+       return 2114;     
+}
+
+Double_t AliGenEMlib::PtDeltaZero( const Double_t *px, const Double_t */*dummy*/ )
+{
+       // Delta0 pT
+       return MtScal(*px,kDeltaZero);
+}
+
+Double_t AliGenEMlib::YDeltaZero( const Double_t *py, const Double_t */*dummy*/ )
+{
+       return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::V2DeltaZero( const Double_t *px, const Double_t */*dummy*/ )
+{
+       return KEtScal(*px,kDeltaZero);
+}
+
+
+//--------------------------------------------------------------------------
+//
+//                              rho +
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpRhoPl(TRandom *)
+{
+       // Return rho+ pdg code
+       return 213;     
+}
+
+Double_t AliGenEMlib::PtRhoPl( const Double_t *px, const Double_t */*dummy*/ )
+{
+       // rho + pT
+       return MtScal(*px,kRhoPl);
+}
+
+Double_t AliGenEMlib::YRhoPl( const Double_t *py, const Double_t */*dummy*/ )
+{
+       return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::V2RhoPl( const Double_t *px, const Double_t */*dummy*/ )
+{
+       return KEtScal(*px,kRhoPl);
+}
+
+
+//--------------------------------------------------------------------------
+//
+//                              rho -
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpRhoMi(TRandom *)
+{
+       // Return rho- pdg code
+       return -213;     
 }
 
+Double_t AliGenEMlib::PtRhoMi( const Double_t *px, const Double_t */*dummy*/ )
+{
+       // rho- pT
+       return MtScal(*px,kRhoMi);
+}
+
+Double_t AliGenEMlib::YRhoMi( const Double_t *py, const Double_t */*dummy*/ )
+{
+       return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::V2RhoMi( const Double_t *px, const Double_t */*dummy*/ )
+{
+       return KEtScal(*px,kRhoMi);
+}
+
+
+//--------------------------------------------------------------------------
+//
+//                             K0 *
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpK0star(TRandom *)
+{
+       // Return K0 * pdg code
+       return 313;     
+}
+
+Double_t AliGenEMlib::PtK0star( const Double_t *px, const Double_t */*dummy*/ )
+{
+       // K0 * pT
+       return MtScal(*px,kK0star);
+}
+
+Double_t AliGenEMlib::YK0star( const Double_t *py, const Double_t */*dummy*/ )
+{
+       return YFlat(*py);
+
+}
+
+Double_t AliGenEMlib::V2K0star( const Double_t *px, const Double_t */*dummy*/ )
+{
+       return KEtScal(*px,kK0star);
+}
+
+
+
 Double_t AliGenEMlib::YFlat(Double_t /*y*/)
 {
-  //--------------------------------------------------------------------------
-  //
-  //                    flat rapidity distribution 
-  //
-  //--------------------------------------------------------------------------
+       //--------------------------------------------------------------------------
+       //
+       //                    flat rapidity distribution 
+       //
+       //--------------------------------------------------------------------------
 
-  Double_t dNdy = 1.;   
+       Double_t dNdy = 1.;   
 
-  return dNdy;
+       return dNdy;
 
 }
 
+
 //============================================================= 
 //
 //                    Mt-scaling  
@@ -866,68 +1528,90 @@ Double_t AliGenEMlib::YFlat(Double_t /*y*/)
 //
 Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
 {
-  // Function for the calculation of the Pt distribution for a 
-  // given particle np, from the pizero Pt distribution using  
-  // mt scaling. 
-
-  Double_t scaledPt = sqrt(pt*pt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
-  Double_t scaledYield = PtPizero(&scaledPt, (Double_t*) 0);
-
-  //     VALUE MESON/PI AT 5 GeV/c
-  Double_t NormPt = 5.;
-  Double_t scaledNormPt = sqrt(NormPt*NormPt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
-
-  Double_t norm = fgkMtFactor[int(bool(fgSelectedCentrality))][np] * (PtPizero(&NormPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
-
-  return norm*(pt/scaledPt)*scaledYield;
+       // Function for the calculation of the Pt distribution for a 
+       // given particle np, from the pizero Pt distribution using  
+       // mt scaling. 
+
+       Double_t scaledPt = sqrt(pt*pt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
+       Double_t scaledYield = PtPizero(&scaledPt, (Double_t*) 0);
+
+       //     VALUE MESON/PI AT 5 GeV/c
+       Double_t NormPt = 5.;
+       Double_t scaledNormPt = sqrt(NormPt*NormPt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
+
+       Int_t selectedCol;
+       switch (fgSelectedCollisionsSystem){
+               case kpp900GeV:
+                       selectedCol=0;
+                       break;
+               case kpp2760GeV:
+                       selectedCol=0;
+                       break;
+               case kpp7TeV:
+                       selectedCol=0;
+                       break;
+               case kpPb:
+                       selectedCol=1;
+                       break;
+               case kPbPb:
+                       selectedCol=2;
+                       break;
+               default:
+                       selectedCol=0;
+                       printf("<AliGenEMlib::MtScal> no collision system has been given\n");
+       }
+       
+       Double_t norm = fgkMtFactor[selectedCol][np] * (PtPizero(&NormPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
+
+       return norm*(pt/scaledPt)*scaledYield;
 }
 
 Double_t AliGenEMlib::KEtScal(Double_t pt, Int_t np)
 {
-  const int nq=2; //number of quarks for particle np, here always 2
-  Double_t scaledPt = sqrt(pow(2.0/nq*(sqrt(pt*pt+fgkHM[np]*fgkHM[np])-fgkHM[np])+fgkHM[0],2)-fgkHM[0]*fgkHM[0]);
-  return V2Pizero(&scaledPt, (Double_t*) 0);
+       const int nq=2; //number of quarks for particle np, here always 2
+       Double_t scaledPt = sqrt(pow(2.0/nq*(sqrt(pt*pt+fgkHM[np]*fgkHM[np])-fgkHM[np])+fgkHM[0],2)-fgkHM[0]*fgkHM[0]);
+       return V2Pizero(&scaledPt, (Double_t*) 0);
 }
 
 Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
 {
-  // Very general parametrization of the v2
-
-  const double &pt=px[0];
-  double val=CrossOverLc(par[4],par[3],pt)*(2*par[0]/(1+TMath::Exp(par[1]*(par[2]-pt)))-par[0])+CrossOverRc(par[4],par[3],pt)*((par[8]-par[5])/(1+TMath::Exp(par[6]*(pt-par[7])))+par[5]);
-  double sys=0;
-  if(fgSelectedV2Systematic){
-    double syspt=pt>par[15]?par[15]:pt;
-    sys=fgSelectedV2Systematic*par[11+fgSelectedV2Systematic*2]*pow(syspt,par[12+fgSelectedV2Systematic*2]);
-  }
-  return std::max(val+sys,0.0);
+       // Very general parametrization of the v2
+
+       const double &pt=px[0];
+       double val=CrossOverLc(par[4],par[3],pt)*(2*par[0]/(1+TMath::Exp(par[1]*(par[2]-pt)))-par[0])+CrossOverRc(par[4],par[3],pt)*((par[8]-par[5])/(1+TMath::Exp(par[6]*(pt-par[7])))+par[5]);
+       double sys=0;
+       if(fgSelectedV2Systematic){
+               double syspt=pt>par[15]?par[15]:pt;
+               sys=fgSelectedV2Systematic*par[11+fgSelectedV2Systematic*2]*pow(syspt,par[12+fgSelectedV2Systematic*2]);
+       }
+       return std::max(val+sys,0.0);
 }
 
 Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
 {
-  // Flat v2
+       // Flat v2
 
-  return 0.0;
+       return 0.0;
 }
 
 Double_t AliGenEMlib::GetTAA(Int_t cent){
-  const static Double_t taa[16] = { 1.0,    // pp
-                                   26.32,  // 0-5
-                                   20.56,  // 5-10
-                                   14.39,  // 10-20
-                                   8.70,   // 20-30
-                                   5.001,  // 30-40
-                                   2.675,  // 40-50
-                                   1.317,  // 50-60
-                                   23.44,  // 0-10
-                                   6.85,   // 20-40
-                                   1.996,  // 40-60
-                                   0.4174, // 60-80
-                                   18.91,  // 0-20
-                                   12.88,  // 0-40
-                                   3.088,  // 20-80
-                                   1.207}; // 40-80
-  return taa[cent];  
+       const static Double_t taa[16] = { 1.0,    // pp
+                                               26.32,  // 0-5
+                                               20.56,  // 5-10
+                                               14.39,  // 10-20
+                                               8.70,   // 20-30
+                                               5.001,  // 30-40
+                                               2.675,  // 40-50
+                                               1.317,  // 50-60
+                                               23.44,  // 0-10
+                                               6.85,   // 20-40
+                                               1.996,  // 40-60
+                                               0.4174, // 60-80
+                                               18.91,  // 0-20
+                                               12.88,  // 0-40
+                                               3.088,  // 20-80
+                                               1.207}; // 40-80
+       return taa[cent];  
 }
 
 //==========================================================================
@@ -942,173 +1626,274 @@ typedef Int_t (*GenFuncIp) (TRandom *);
 
 GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
 {
-  // Return pointer to pT parameterisation
-   GenFunc func=0;
-   TString sname(tname);
-
-   switch (param) 
-    {
-     case kDirectRealGamma:
-       func=PtDirectRealGamma;
-      break;
-     case kDirectVirtGamma:
-       func=PtDirectVirtGamma;
-      break;
-    case kPizero:
-      func=PtPizero;
-      break;
-    case kEta:
-      func=PtEta;
-      break;
-    case kRho:
-      func=PtRho;
-      break;
-    case kOmega:
-      func=PtOmega;
-      break;
-    case kEtaprime:
-      func=PtEtaprime;
-      break;
-    case kPhi:
-      func=PtPhi;
-      break;
-    case kJpsi:
-      func=PtJpsi;
-      break;
-
-    default:
-      func=0;
-      printf("<AliGenEMlib::GetPt> unknown parametrisation\n");
+       // Return pointer to pT parameterisation
+       GenFunc func=0;
+       TString sname(tname);
+
+       switch (param) {
+               case kDirectRealGamma:
+                       func=PtDirectRealGamma;
+                       break;
+               case kDirectVirtGamma:
+                       func=PtDirectVirtGamma;
+                       break;
+               case kPizero:
+                       func=PtPizero;
+                       break;
+               case kEta:
+                       func=PtEta;
+                       break;
+               case kRho0:
+                       func=PtRho0;
+                       break;
+               case kOmega:
+                       func=PtOmega;
+                       break;
+               case kEtaprime:
+                       func=PtEtaprime;
+                       break;
+               case kPhi:
+                       func=PtPhi;
+                       break;
+               case kJpsi:
+                       func=PtJpsi;
+                       break;
+               case kSigma0:
+                       func= PtSigma;
+                       break;
+               case kK0s:
+                       func= PtK0short;
+                       break;
+               case kDeltaPlPl:
+                       func= PtDeltaPlPl;
+                       break;
+               case kDeltaPl:
+                       func= PtDeltaPlPl;
+                       break;
+               case kDeltaMi:
+                       func= PtDeltaMi;
+                       break;
+               case kDeltaZero:
+                       func= PtDeltaZero;
+                       break;
+               case kRhoPl:
+                       func= PtRhoPl;
+                       break;
+               case kRhoMi:
+                       func= PtRhoMi;
+                       break;
+               case kK0star:
+                       func= PtK0star;
+                       break;
+               default:
+                       func=0;
+                       printf("<AliGenEMlib::GetPt> unknown parametrisation\n");
     }
-   return func;
+       return func;
 }
 
 GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
 {
-  // Return pointer to y- parameterisation
-   GenFunc func=0;
-   TString sname(tname);
-
-   switch (param) 
-    {
-    case kDirectRealGamma:
-      func=YDirectRealGamma;
-      break;
-    case kDirectVirtGamma:
-      func=YDirectVirtGamma;
-      break;
-    case kPizero:
-         func=YPizero;
-         break;
-    case kEta:
-         func=YEta;
-         break;
-    case kRho:
-         func=YRho;
-         break;
-    case kOmega:
-         func=YOmega;
-         break;
-    case kEtaprime:
-         func=YEtaprime;
-         break;
-    case kPhi:
-         func=YPhi;
-         break;
-    case kJpsi:
-      func=YJpsi;
-      break;
-
-    default:
-        func=0;
-        printf("<AliGenEMlib::GetY> unknown parametrisation\n");
-    }
-    return func;
+       // Return pointer to y- parameterisation
+       GenFunc func=0;
+       TString sname(tname);
+
+       switch (param) {
+               case kDirectRealGamma:
+                       func=YDirectRealGamma;
+                       break;
+               case kDirectVirtGamma:
+                       func=YDirectVirtGamma;
+                       break;
+               case kPizero:
+                       func=YPizero;
+                       break;
+               case kEta:
+                       func=YEta;
+                       break;
+               case kRho0:
+                       func=YRho0;
+                       break;
+               case kOmega:
+                       func=YOmega;
+                       break;
+               case kEtaprime:
+                       func=YEtaprime;
+                       break;
+               case kPhi:
+                       func=YPhi;
+                       break;
+               case kJpsi:
+                       func=YJpsi;
+                       break;
+               case kSigma0:
+                       func=YSigma;
+                       break;           
+               case kK0s:
+                       func=YK0short;
+                       break;
+               case kDeltaPlPl:
+                       func=YDeltaPlPl;
+                       break;
+               case kDeltaPl:
+                       func=YDeltaPl;
+                       break;
+               case kDeltaMi:
+                       func=YDeltaMi;
+                       break;
+               case kDeltaZero:
+                       func=YDeltaZero;
+                       break;
+               case kRhoPl:
+                       func=YRhoPl;
+                       break;
+               case kRhoMi:
+                       func=YRhoMi;
+                       break;
+               case kK0star:
+                       func=YK0star;
+                       break;
+               default:
+                       func=0;
+                       printf("<AliGenEMlib::GetY> unknown parametrisation\n");
+       }
+       return func;
 }
 
 GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
 {
-  // Return pointer to particle type parameterisation
-   GenFuncIp func=0;
-   TString sname(tname);
-
-   switch (param) 
-    {
-    case kDirectRealGamma:
-      func=IpDirectRealGamma;
-      break;
-    case kDirectVirtGamma:
-      func=IpDirectVirtGamma;
-      break;
-    case kPizero:
-         func=IpPizero;
-         break;
-    case kEta:
-         func=IpEta;
-         break;
-    case kRho:
-         func=IpRho;
-         break;
-    case kOmega:
-         func=IpOmega;
-         break;
-    case kEtaprime:
-         func=IpEtaprime;
-         break;
-    case kPhi:
-         func=IpPhi;
-         break;
-    case kJpsi:
-      func=IpJpsi;
-      break;
-
-    default:
-        func=0;
-        printf("<AliGenEMlib::GetIp> unknown parametrisation\n");
-    }
-    return func;
+       // Return pointer to particle type parameterisation
+       GenFuncIp func=0;
+       TString sname(tname);
+
+       switch (param) {
+               case kDirectRealGamma:
+                       func=IpDirectRealGamma;
+                       break;
+               case kDirectVirtGamma:
+                       func=IpDirectVirtGamma;
+                       break;
+               case kPizero:
+                       func=IpPizero;
+                       break;
+               case kEta:
+                       func=IpEta;
+                       break;
+               case kRho0:
+                       func=IpRho0;
+                       break;
+               case kOmega:
+                       func=IpOmega;
+                       break;
+               case kEtaprime:
+                       func=IpEtaprime;
+                       break;
+               case kPhi:
+                       func=IpPhi;
+                       break;
+               case kJpsi:
+                       func=IpJpsi;
+                       break;
+               case kSigma0:
+                       func=IpSigma;
+                       break; 
+               case kK0s:
+                       func=IpK0short;
+                       break;
+               case kDeltaPlPl:
+                       func=IpDeltaPlPl;
+                       break;
+               case kDeltaPl:
+                       func=IpDeltaPl;
+                       break;
+               case kDeltaMi:
+                       func=IpDeltaMi;
+                       break;
+               case kDeltaZero:
+                       func=IpDeltaZero;
+                       break;
+               case kRhoPl:
+                       func=IpRhoPl;
+                       break;
+               case kRhoMi:
+                       func=IpRhoMi;
+                       break;
+               case kK0star:
+                       func=IpK0star;
+                       break;
+               default:
+                       func=0;
+                       printf("<AliGenEMlib::GetIp> unknown parametrisation\n");
+       }
+       return func;
 }
 
 GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
 {
-  // Return pointer to v2-parameterisation
-  GenFunc func=0;
-  TString sname(tname);
-
-  switch (param) 
-    {
-    case kDirectRealGamma:
-      func=V2DirectRealGamma;
-      break;
-    case kDirectVirtGamma:
-      func=V2DirectVirtGamma;
-      break;
-    case kPizero:
-      func=V2Pizero;
-      break;
-    case kEta:
-      func=V2Eta;
-      break;
-    case kRho:
-      func=V2Pizero;
-      break;
-    case kOmega:
-      func=V2Pizero;
-      break;
-    case kEtaprime:
-      func=V2Pizero;
-      break;
-    case kPhi:
-      func=V2Pizero;
-      break;
-    case kJpsi:
-      func=V2Jpsi;
-      break;
-
-    default:
-      func=0;
-      printf("<AliGenEMlib::GetV2> unknown parametrisation\n");
-    }
-  return func;
+       // Return pointer to v2-parameterisation
+       GenFunc func=0;
+       TString sname(tname);
+
+       switch (param) {
+               case kDirectRealGamma:
+                       func=V2DirectRealGamma;
+                       break;
+               case kDirectVirtGamma:
+                       func=V2DirectVirtGamma;
+                       break;
+               case kPizero:
+                       func=V2Pizero;
+                       break;
+               case kEta:
+                       func=V2Eta;
+                       break;
+               case kRho0:
+                       func=V2Rho0;
+                       break;
+               case kOmega:
+                       func=V2Omega;
+                       break;
+               case kEtaprime:
+                       func=V2Etaprime;
+                       break;
+               case kPhi:
+                       func=V2Phi;
+                       break;
+               case kJpsi:
+                       func=V2Jpsi;
+                       break;
+               case kSigma0:
+                       func=V2Sigma0;
+                       break;    
+               case kK0s:
+                       func=V2K0sshort;
+                       break;
+               case kDeltaPlPl:
+                       func=V2DeltaPlPl;
+                       break;
+               case kDeltaPl:
+                       func=V2DeltaPl;
+                       break;
+               case kDeltaMi:
+                       func=V2DeltaMi;
+                       break;
+               case kDeltaZero:
+                       func=V2DeltaZero;
+                       break;
+               case kRhoPl:
+                       func=V2RhoPl;
+                       break;
+               case kRhoMi:
+                       func=V2RhoMi;
+                       break;
+               case kK0star:
+                       func=V2K0star;
+                       break;
+
+               default:
+                       func=0;
+                       printf("<AliGenEMlib::GetV2> unknown parametrisation\n");
+       }
+       return func;
 }
 
index 2e4742e..a02e7e9 100644 (file)
 class TRandom;
 
 class AliGenEMlib :public AliGenLib {
-public:
-    
-  enum Particle_t{kPizero=0, kEta, kRho, kOmega, kEtaprime, kPhi, kJpsi, kDirectRealGamma, kDirectVirtGamma };
-  enum Centrality_t{kpp=0x0, k0005=0x1, k0510=0x2, k1020=0x3, k2030=0x4, k3040=0x5, k4050=0x6, k5060=0x7, k0010=0x8, k2040=0x9, k4060=0xA, k6080=0xB, k0020=0xC, k0040=0xD, k2080=0xE, k4080=0xF, kCentralities=0x10};
-  enum PtParamSet_t{kPizero7TeVpp=0x00, kPizeroEta7TeVpp=0x10, kPizero7TeVpplow=0x20, kPizeroEta7TeVpplow=0x30, kPizero7TeVpphigh=0x40, kPizeroEta7TeVpphigh=0x50, kPizero2760GeVpp=0x60, kPizeroEta2760GeVpp=0x70, kPizero2760GeVpplow=0x80, kPizeroEta2760GeVpplow=0x90, kPizero2760GeVpphigh=0xA0, kPizeroEta2760GeVpphigh=0xB0, kPiOldChargedPbPb=0xC0, kPichargedPbPb=0xD0, kPizeroPbPb=0xE0, kPichargedPPb=0xF0 };
-  enum v2Sys_t{kLoV2Sys=-1, kNoV2Sys=0, kUpV2Sys=+1};
-  
-  AliGenEMlib() { } ;
-
-  static void SelectParams(PtParamSet_t ptSelect, Centrality_t centSelect=kpp, Int_t v2sys=kNoV2Sys) 
-  { fgSelectedPtParam=ptSelect; fgSelectedCentrality=centSelect; fgSelectedV2Systematic=v2sys; }
-
-    GenFunc   GetPt(Int_t param, const char * tname=0) const;
-    GenFunc   GetY(Int_t param, const char * tname=0) const;
-    GenFuncIp GetIp(Int_t param, const char * tname=0) const;    
-  GenFunc   GetV2(Int_t param, const char * tname=0) const;
-
-  //private:
-
-  // General functions
-
-  static Int_t fgSelectedPtParam; // selected pT parameter
-  static Int_t fgSelectedCentrality; // selected Centrality
-  static Int_t fgSelectedV2Systematic; // selected v2 systematics, usefully values: -1,0,1
-
-
-  static Double_t PtModifiedHagedornThermal(Double_t pt, 
-                                           Double_t c, 
-                                           Double_t p0, 
-                                           Double_t p1, 
-                                           Double_t n,
-                                           Double_t cT,
-                                           Double_t T);
-
-
-  static Double_t PtModifiedHagedornExp(Double_t pt,
-                                       Double_t c,
-                                       Double_t p0,
-                                       Double_t p1,
-                                       Double_t p2,
-                                       Double_t n); 
-
-
-  static Double_t PtModifiedHagedornExp2(Double_t pt,
-                                           Double_t c,
-                                           Double_t a,
-                                           Double_t b,
-                                           Double_t p0,
-                                           Double_t p1,
-                                           Double_t d,
-                                           Double_t n);
-
-
-  static Double_t PtTsallis(Double_t pt,
-                           Double_t m,
-                           Double_t c,
-                           Double_t T,
-                           Double_t n);
-
-  static Double_t PtExponential(const Double_t *pt, const Double_t *param);
-  static Double_t PtModifiedHagedornPowerlaw(const Double_t *pt, const Double_t *param);
-  static Double_t PtDoublePowerlaw(const Double_t *pt, const Double_t *param);
-  static Double_t IntegratedKrollWada(const Double_t *mh, const Double_t *);
-
-  // direct gamma
-  static Double_t PtPromptRealGamma(const Double_t *px, const Double_t *dummy);
-  static Double_t PtPromptVirtGamma(const Double_t *px, const Double_t *dummy);
-  static Double_t PtThermalRealGamma(const Double_t *px, const Double_t *dummy);
-  static Double_t PtThermalVirtGamma(const Double_t *px, const Double_t *dummy);
-
-  static Int_t    IpDirectRealGamma(TRandom *ran);
-  static Double_t PtDirectRealGamma(const Double_t *px, const Double_t *dummy);
-  static Double_t YDirectRealGamma(const Double_t *py, const Double_t *dummy);
-  static Double_t V2DirectRealGamma(const Double_t *px, const Double_t *dummy);
-
-  static Int_t    IpDirectVirtGamma(TRandom *ran);
-  static Double_t PtDirectVirtGamma(const Double_t *px, const Double_t *dummy);
-  static Double_t YDirectVirtGamma(const Double_t *py, const Double_t *dummy);
-  static Double_t V2DirectVirtGamma(const Double_t *px, const Double_t *dummy);
-
-  // Pizero
-    static Int_t    IpPizero(TRandom *ran);
-  static Double_t PtPizero(const Double_t *px, const Double_t *dummy);
-    static Double_t YPizero(const Double_t *py, const Double_t *dummy);
-  static Double_t V2Pizero(const Double_t *px, const Double_t *dummy);
-
-  // Eta
-    static Int_t    IpEta(TRandom *ran);
-  static Double_t PtEta(const Double_t *px, const Double_t *dummy);
-    static Double_t YEta(const Double_t *py, const Double_t *dummy);
-  static Double_t V2Eta(const Double_t *px, const Double_t *dummy);
-
-  // Rho
-    static Int_t    IpRho(TRandom *ran);
-  static Double_t PtRho(const Double_t *px, const Double_t *dummy);
-    static Double_t YRho(const Double_t *py, const Double_t *dummy);
-  static Double_t V2Rho(const Double_t *py, const Double_t *dummy);
-
-  // Omega
-    static Int_t    IpOmega(TRandom *ran);
-  static Double_t PtOmega(const Double_t *px, const Double_t *dummy);
-    static Double_t YOmega(const Double_t *py, const Double_t *dummy);
-  static Double_t V2Omega(const Double_t *py, const Double_t *dummy);
-
-  // Etaprime
-    static Int_t    IpEtaprime(TRandom *ran);
-  static Double_t PtEtaprime(const Double_t *px, const Double_t *dummy);
-    static Double_t YEtaprime(const Double_t *py, const Double_t *dummy);
-  static Double_t V2Etaprime(const Double_t *py, const Double_t *dummy);
-
-  // Phi
-    static Int_t    IpPhi(TRandom *ran);
-  static Double_t PtPhi(const Double_t *px, const Double_t *dummy);
-    static Double_t YPhi(const Double_t *py, const Double_t *dummy);
-  static Double_t V2Phi(const Double_t *py, const Double_t *dummy);
-
-  // Jpsi
-  static Int_t    IpJpsi(TRandom *ran);
-  static Double_t PtJpsi(const Double_t *px, const Double_t *dummy);
-  static Double_t YJpsi(const Double_t *py, const Double_t *dummy);
-  static Double_t V2Jpsi(const Double_t *py, const Double_t *dummy);
-
-  // General
-  //static Double_t PtFlat(const Double_t *px, const Double_t *dummy);
-    static Double_t YFlat(Double_t y);
-  static Double_t MtScal(Double_t pt, Int_t np);
-  static Double_t V2Param(const Double_t *px, const Double_t *param);
-  static Double_t V2Flat(const Double_t *px, const Double_t *param);
-  static Double_t KEtScal(Double_t pt, Int_t np);
-  static Double_t GetTAA(Int_t cent);
-
-  static Double_t CrossOverLc(double a, double b, double x);
-  static Double_t CrossOverRc(double a, double b, double x);
-
-  static const Double_t fgkPtParam[kCentralities][10];          // parameters of pi pt spectrum
-  static const Double_t fgkV2param[kCentralities][16];          // parameters of pi v2 spectrum
-  static const Double_t fgkRawPtOfV2Param[kCentralities][10];   // parameters of the raw pt spectrum of v2 analysis
-  static const Double_t fgkThermPtParam[kCentralities][2];      // parameters of thermal gamma pt
-  static const Double_t fgkHM[8];                    // particle masses
-  static const Double_t fgkMtFactor[2][8];           // mt scaling factor
-
-  ClassDef(AliGenEMlib,0)
+       public:
+
+               enum Particle_t{kPizero=0x0, kEta=0x1, kRho0=0x2, kOmega=0x3, kEtaprime=0x4, kPhi=0x5, kJpsi=0x6, kSigma0=0x7, kK0s=0x8, kDeltaPlPl=0x9, kDeltaPl=0xA, kDeltaMi=0xB, kDeltaZero=0xC, kRhoPl = 0xD, kRhoMi = 0xE, kK0star = 0xF,  kDirectRealGamma=0x10, kDirectVirtGamma=0x11  };
+               enum Centrality_t{ kpp = 0x0, 
+                                                  k0005=0x1, k0510=0x2, k1020=0x3, k2030=0x4, k3040=0x5, k4050=0x6, k5060=0x7, k0010=0x8, k2040=0x9, k4060=0xA, k6080=0xB, k0020=0xC, k0040=0xD, k2080=0xE, k4080=0xF, kCentralities=0x10};
+               enum CollisionSystem_t {kpp900GeV=0x0, kpp2760GeV=0x1, kpp7TeV=0x2, kpPb=0x3, kPbPb=0x4};
+               enum PtParamSetPi0_t{kPizeroParam=0x0, kPizeroParamlow=0x1, kPizeroParamhigh=0x2, 
+                                                       kPichargedParam=0x3, kPichargedParamlow=0x4, kPichargedParamhigh=0x5, 
+                                                       kPizeroParamAlter=0x6, kPizeroParamAlterlow=0x7, kPizeroParamAlterhigh=0x8,
+                                                       kNPi0Param=0x9
+               };
+               enum PtParamSetEta_t{kEtaParampp=0x0, kEtaParampplow=0x1, kEtaParampphigh=0x2,
+                                                        kEtaParamRatiopp=0x3, kEtaParamRatiopplow=0x4, kEtaParamRatiopphigh=0x5,
+                                                        kEtaParamPbPb=0x6,
+                                                        kEtaParamPPb=0x7
+               };
+               enum PtParamSetOmega_t{kOmegaParampp=0x0, kOmegaParampplow=0x1, kOmegaParampphigh=0x2,
+                                                          kOmegaParamRatiopp=0x3, kOmegaParamRatiopplow=0x4, kOmegaParamRatiopphigh=0x5,
+                                                        kOmegaParamPbPb=0x6,
+                                                        kOmegaParamPPb=0x7
+               };
+               enum PtParamSetPhi_t{kPhiParampp=0x0, kPhiParampplow=0x1, kPhiParampphigh=0x2,
+                                                        kPhiParamPbPb=0x3,
+                                                        kPhiParamPPb=0x4
+               };
+               
+               
+               enum v2Sys_t{kLoV2Sys=-1, kNoV2Sys=0, kUpV2Sys=+1};
+       
+               AliGenEMlib() { } ;
+
+               static void SelectParams( Int_t collisionSystem,
+                                                                 Int_t ptSelectPi0, 
+                                                                 Int_t ptSelectEta, 
+                                                                 Int_t ptSelectOmega,
+                                                                 Int_t ptSelectPhi, 
+                                                             Int_t centSelect=kpp, 
+                                                             Int_t v2sys=kNoV2Sys) {   
+                       fgSelectedCollisionsSystem = collisionSystem;
+                       fgSelectedPtParamPi0=ptSelectPi0; 
+                       fgSelectedPtParamEta=ptSelectEta; 
+                       fgSelectedPtParamOmega=ptSelectOmega; 
+                       fgSelectedPtParamPhi=ptSelectPhi; 
+                       fgSelectedCentrality=centSelect; 
+                       fgSelectedV2Systematic=v2sys; 
+                                                                         
+               }
+
+               GenFunc   GetPt(Int_t param, const char * tname=0) const;
+               GenFunc   GetY(Int_t param, const char * tname=0) const;
+               GenFuncIp GetIp(Int_t param, const char * tname=0) const;    
+               GenFunc   GetV2(Int_t param, const char * tname=0) const;
+
+       //private:
+
+       // General functions
+
+               // General functions
+               static Int_t fgSelectedCollisionsSystem; // selected pT parameter
+               static Int_t fgSelectedPtParamPi0; // selected pT parameter
+               static Int_t fgSelectedPtParamEta; // selected pT parameter
+               static Int_t fgSelectedPtParamOmega; // selected pT parameter
+               static Int_t fgSelectedPtParamPhi; // selected pT parameter
+               static Int_t fgSelectedCentrality; // selected Centrality
+               static Int_t fgSelectedV2Systematic; // selected v2 systematics, usefully values: -1,0,1
+
+
+               static Double_t PtModifiedHagedornThermal(Double_t pt, 
+                                                               Double_t c, 
+                                                               Double_t p0, 
+                                                               Double_t p1, 
+                                                               Double_t n,
+                                                               Double_t cT,
+                                                               Double_t T);
+
+
+               
+               static Double_t PtModifiedHagedornExp(Double_t pt,
+                                                       Double_t c,
+                                                       Double_t p0,
+                                                       Double_t p1,
+                                                       Double_t p2,
+                                                       Double_t n); 
+
+
+               static Double_t PtModifiedHagedornExp2(Double_t pt,
+                                                                                               Double_t c,
+                                                                                               Double_t a,
+                                                                                               Double_t b,
+                                                                                               Double_t p0,
+                                                                                               Double_t p1,
+                                                                                               Double_t d,
+                                                                                               Double_t n);
+
+
+               static Double_t PtTsallis(Double_t pt,
+                                               Double_t m,
+                                               Double_t c,
+                                               Double_t T,
+                                               Double_t n);
+
+               static Double_t PtParticleRatiopp(Double_t pt,
+                                                                               Double_t m1,
+                                                                               Double_t m2,
+                                                                               Double_t c1,
+                                                                               Double_t c2,
+                                                                               Double_t T1,
+                                                                               Double_t T2,
+                                                                               Double_t n);
+               
+               static Double_t PtXQCD( Double_t pt,
+                                Double_t a,
+                                Double_t b,
+                                Double_t c,
+                                Double_t d,
+                                                               Double_t e,
+                                                               Double_t f);
+               
+               static Double_t PtModTsallis(   Double_t pt,
+                                Double_t a,
+                                Double_t b,
+                                Double_t c,
+                                Double_t d,
+                                                           Double_t e,
+                                                               Double_t f,
+                                                               Double_t g,
+                                                               Double_t mass);
+               
+               static Double_t PtQCD(  Double_t pt,
+                                Double_t a,
+                                Double_t b,
+                                Double_t c,
+                                Double_t d,
+                                                           Double_t e);
+
+               static Double_t PtExponential(const Double_t *pt, const Double_t *param);
+               static Double_t PtModifiedHagedornPowerlaw(const Double_t *pt, const Double_t *param);
+               static Double_t PtDoublePowerlaw(const Double_t *pt, const Double_t *param);
+               static Double_t IntegratedKrollWada(const Double_t *mh, const Double_t *);
+
+               // direct gamma
+               static Double_t PtPromptRealGamma(const Double_t *px, const Double_t *dummy);
+               static Double_t PtPromptVirtGamma(const Double_t *px, const Double_t *dummy);
+               static Double_t PtThermalRealGamma(const Double_t *px, const Double_t *dummy);
+               static Double_t PtThermalVirtGamma(const Double_t *px, const Double_t *dummy);
+
+               static Int_t    IpDirectRealGamma(TRandom *ran);
+               static Double_t PtDirectRealGamma(const Double_t *px, const Double_t *dummy);
+               static Double_t YDirectRealGamma(const Double_t *py, const Double_t *dummy);
+               static Double_t V2DirectRealGamma(const Double_t *px, const Double_t *dummy);
+
+               static Int_t    IpDirectVirtGamma(TRandom *ran);
+               static Double_t PtDirectVirtGamma(const Double_t *px, const Double_t *dummy);
+               static Double_t YDirectVirtGamma(const Double_t *py, const Double_t *dummy);
+               static Double_t V2DirectVirtGamma(const Double_t *px, const Double_t *dummy);
+
+               // Pizero
+               static Int_t    IpPizero(TRandom *ran);
+               static Double_t PtPizero(const Double_t *px, const Double_t *dummy);
+               static Double_t YPizero(const Double_t *py, const Double_t *dummy);
+               static Double_t V2Pizero(const Double_t *px, const Double_t *dummy);
+
+               // Eta
+               static Int_t    IpEta(TRandom *ran);
+               static Double_t PtEta(const Double_t *px, const Double_t *dummy);
+               static Double_t YEta(const Double_t *py, const Double_t *dummy);
+               static Double_t V2Eta(const Double_t *px, const Double_t *dummy);
+
+               // Rho
+               static Int_t    IpRho0(TRandom *ran);
+               static Double_t PtRho0(const Double_t *px, const Double_t *dummy);
+               static Double_t YRho0(const Double_t *py, const Double_t *dummy);
+               static Double_t V2Rho0(const Double_t *py, const Double_t *dummy);
+
+
+               // Omega
+               static Int_t    IpOmega(TRandom *ran);
+               static Double_t PtOmega(const Double_t *px, const Double_t *dummy);
+               static Double_t YOmega(const Double_t *py, const Double_t *dummy);
+               static Double_t V2Omega(const Double_t *py, const Double_t *dummy);
+
+               // Etaprime
+               static Int_t    IpEtaprime(TRandom *ran);
+               static Double_t PtEtaprime(const Double_t *px, const Double_t *dummy);
+               static Double_t YEtaprime(const Double_t *py, const Double_t *dummy);
+               static Double_t V2Etaprime(const Double_t *py, const Double_t *dummy);
+
+               // Phi
+               static Int_t    IpPhi(TRandom *ran);
+               static Double_t PtPhi(const Double_t *px, const Double_t *dummy);
+               static Double_t YPhi(const Double_t *py, const Double_t *dummy);
+               static Double_t V2Phi(const Double_t *py, const Double_t *dummy);
+
+               // Jpsi
+               static Int_t    IpJpsi(TRandom *ran);
+               static Double_t PtJpsi(const Double_t *px, const Double_t *dummy);
+               static Double_t YJpsi(const Double_t *py, const Double_t *dummy);
+               static Double_t V2Jpsi(const Double_t *py, const Double_t *dummy);
+
+               // Sigma
+               static Int_t    IpSigma(TRandom *ran);
+               static Double_t PtSigma( const Double_t *px, const Double_t *dummy );
+               static Double_t YSigma(const Double_t *py, const Double_t *dummy);
+               static Double_t V2Sigma0( const Double_t *px, const Double_t *dummy );
+               
+               // K0short
+               static Int_t    IpK0short(TRandom *ran);
+               static Double_t PtK0short( const Double_t *px, const Double_t *dummy );
+               static Double_t YK0short(const Double_t *py, const Double_t *dummy);
+               static Double_t V2K0sshort( const Double_t *px, const Double_t *dummy );
+
+               // Delta++
+               static Int_t    IpDeltaPlPl(TRandom *ran);
+               static Double_t PtDeltaPlPl( const Double_t *px, const Double_t *dummy );
+               static Double_t YDeltaPlPl(const Double_t *py, const Double_t *dummy);
+               static Double_t V2DeltaPlPl( const Double_t *px, const Double_t *dummy );
+               
+               // Delta+
+               static Int_t    IpDeltaPl(TRandom *ran);
+               static Double_t PtDeltaPl( const Double_t *px, const Double_t *dummy );
+               static Double_t YDeltaPl(const Double_t *py, const Double_t *dummy);
+               static Double_t V2DeltaPl( const Double_t *px, const Double_t *dummy );
+               
+               // Delta-
+               static Int_t    IpDeltaMi(TRandom *ran);
+               static Double_t PtDeltaMi( const Double_t *px, const Double_t *dummy );
+               static Double_t YDeltaMi(const Double_t *py, const Double_t *dummy);
+               static Double_t V2DeltaMi( const Double_t *px, const Double_t *dummy );
+               
+               // Delta0
+               static Int_t    IpDeltaZero(TRandom *ran);
+               static Double_t PtDeltaZero( const Double_t *px, const Double_t *dummy );
+               static Double_t YDeltaZero(const Double_t *py, const Double_t *dummy);
+               static Double_t V2DeltaZero( const Double_t *px, const Double_t *dummy );
+
+               // Rho+
+               static Int_t    IpRhoPl(TRandom *ran);
+               static Double_t PtRhoPl( const Double_t *px, const Double_t *dummy );
+               static Double_t YRhoPl(const Double_t *py, const Double_t *dummy);
+               static Double_t V2RhoPl( const Double_t *px, const Double_t *dummy );
+
+               // Rho-
+               static Int_t    IpRhoMi(TRandom *ran);
+               static Double_t PtRhoMi( const Double_t *px, const Double_t *dummy );
+               static Double_t YRhoMi(const Double_t *py, const Double_t *dummy);
+               static Double_t V2RhoMi( const Double_t *px, const Double_t *dummy );
+
+               // K0*
+               static Int_t    IpK0star(TRandom *ran);
+               static Double_t PtK0star( const Double_t *px, const Double_t *dummy );
+               static Double_t YK0star(const Double_t *py, const Double_t *dummy);
+               static Double_t V2K0star( const Double_t *px, const Double_t *dummy );
+               
+
+               // General
+               //static Double_t PtFlat(const Double_t *px, const Double_t *dummy);
+                       static Double_t YFlat(Double_t y);
+               static Double_t MtScal(Double_t pt, Int_t np);
+               static Double_t V2Param(const Double_t *px, const Double_t *param);
+               static Double_t V2Flat(const Double_t *px, const Double_t *param);
+               static Double_t KEtScal(Double_t pt, Int_t np);
+               static Double_t GetTAA(Int_t cent);
+
+               static Double_t CrossOverLc(double a, double b, double x);
+               static Double_t CrossOverRc(double a, double b, double x);
+
+               static const Double_t fgkPtParam[kCentralities][10];                                    // parameters of pi pt spectrum
+               static const Double_t fgkModTsallisParamPi0PbPb[kCentralities][7];                      // parameters for ModTsallis function for pi0 in PbPb 
+               static const Double_t fgkModTsallisParamPiChargedPbPb[kCentralities][7];        // parameters for ModTsallis function for pi+- in PbPb 
+               static const Double_t fgkV2param[kCentralities][16];                                    // parameters of pi v2
+               static const Double_t fgkRawPtOfV2Param[kCentralities][10];                             // parameters of the raw pt spectrum of v2 analysys
+               static const Double_t fgkThermPtParam[kCentralities][2];                                        // parameters of thermal gamma pt
+               static const Double_t fgkHM[16];                                                                                        // particle masses
+               static const Double_t fgkMtFactor[3][16];                                                                       // mt scaling factor
+               static const Double_t fgkParamSetPi07TeV[kNPi0Param][7];                                        // parameters for pi0 in 7 TeV
+               static const Double_t fgkParamSetPi02760GeV[kNPi0Param][7];                                     // parameters for pi0 in 2.76 TeV
+               static const Double_t fgkParamSetPi0900GeV[kNPi0Param][7];                                      // parameters for pi0 in 0.9 TeV
+
+               ClassDef(AliGenEMlib,1)
+
 };