Changes concerning the cocktail simulations used for heavy flavour v2 analyses.
authormorsch <andreas.morsch@cern.ch>
Thu, 6 Feb 2014 15:36:37 +0000 (16:36 +0100)
committermorsch <andreas.morsch@cern.ch>
Thu, 6 Feb 2014 15:36:37 +0000 (16:36 +0100)
The update does the following things:
PYTHIA6/AliDecayerPythia.cxx:
- enable J/Psi decays into electrons and gammas
EVGEN/*:
- use well defined pT and v2 spectra parametrizations for the pions
- implement KET scaling
- implement "forced decay" of gammas

Theodor Rascanu

EVGEN/AliGenEMCocktail.cxx
EVGEN/AliGenEMCocktail.h
EVGEN/AliGenEMlib.cxx
EVGEN/AliGenEMlib.h
EVGEN/AliGenParam.cxx
EVGEN/AliGenParam.h
PYTHIA6/AliDecayerPythia.cxx

index 011c074..5bb28da 100644 (file)
@@ -15,7 +15,7 @@
 
 /* $Id: AliGenEMCocktail.cxx 40702 2010-04-26 13:09:52Z morsch $ */
 
-// Class to create cocktails for physics with electrons, di-electrons,
+// Class to create the cocktail for physics with electrons, di-electrons,
 // and photons from the decay of the following sources:
 // pizero, eta, rho, omega, etaprime, phi
 // Kinematic distributions of the sources are taken from AliGenEMlib.
@@ -25,6 +25,7 @@
 // or they are generated according to the pT distributions themselves
 // (weighting mode: kAnalog)  
  
 #include <TObjArray.h>
 #include <TParticle.h>
 #include <TF1.h>
@@ -49,33 +50,38 @@ ClassImp(AliGenEMCocktail)
   
 //________________________________________________________________________
 AliGenEMCocktail::AliGenEMCocktail()
-  :AliGenCocktail(),
+:AliGenCocktail(),
    fDecayer(0),
    fDecayMode(kAll),
    fWeightingMode(kNonAnalog),
    fNPart(1000),
-   fYieldArray()
+  fYieldArray(),
+  fPtSelect(0),
+  fCentrality(0),
+  fV2Systematic(0),
+  fForceConv(kFALSE),
+  fHeaviestParticle(kGENs)
 {
-// Constructor
+  // Constructor
 
 }
 
 //_________________________________________________________________________
 AliGenEMCocktail::~AliGenEMCocktail()
 {
-// Destructor
+  // Destructor
 
 }
 
 //_________________________________________________________________________
 void AliGenEMCocktail::CreateCocktail()
 {
-// create and add sources to the cocktail
+  // create and add sources to the cocktail
 
   fDecayer->SetForceDecay(fDecayMode);
   fDecayer->ForceDecay();
 
-// Set kinematic limits
+  // Set kinematic limits
   Double_t ptMin  = fPtMin;
   Double_t ptMax  = fPtMax;
   Double_t yMin   = fYMin;;
@@ -85,80 +91,108 @@ void AliGenEMCocktail::CreateCocktail()
   AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax));
   AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
 
-// Create and add electron sources to the generator
+  //Initialize user selection for Pt Parameterization and centrality:
+  AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
+
+  // Create and add electron sources to the generator
 
-// pizero
-  AliGenParam * genpizero=0;
+  // pizero
+  AliGenParam *genpizero=0;
   Char_t namePizero[10];    
-  snprintf(namePizero,10, "Pizero");    
+  snprintf(namePizero,10,"Pizero");    
   genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
   AddSource2Generator(namePizero,genpizero);
   TF1 *fPtPizero = genpizero->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
   fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
+  fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
-// eta  
-  AliGenParam * geneta=0;
+
+  // eta  
+  if(fHeaviestParticle<kGenEta)return;
+  AliGenParam *geneta=0;
   Char_t nameEta[10];    
-  snprintf(nameEta,10, "Eta");    
+  snprintf(nameEta,10,"Eta");    
   geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
   AddSource2Generator(nameEta,geneta);
   TF1 *fPtEta = geneta->GetPt();
-#if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
-  fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
-#else
+#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
+#else
+  fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
-// rho  
-  AliGenParam * genrho=0;
+
+  // rho  
+  if(fHeaviestParticle<kGenRho)return;
+  AliGenParam *genrho=0;
   Char_t nameRho[10];    
-  snprintf(nameRho,10, "Rho");    
+  snprintf(nameRho,10,"Rho");    
   genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
   AddSource2Generator(nameRho,genrho);
   TF1 *fPtRho = genrho->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
+  fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
-// omega
-  AliGenParam * genomega=0;
+  
+  // omega
+  if(fHeaviestParticle<kGenOmega)return;
+  AliGenParam *genomega=0;
   Char_t nameOmega[10];    
-  snprintf(nameOmega,10, "Omega");    
+  snprintf(nameOmega,10,"Omega");    
   genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
   AddSource2Generator(nameOmega,genomega);
   TF1 *fPtOmega = genomega->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(double*)0,1.e-6);
+  fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
-// etaprime
-  AliGenParam * genetaprime=0;
+
+  // etaprime
+  if(fHeaviestParticle<kGenEtaprime)return;
+  AliGenParam *genetaprime=0;
   Char_t nameEtaprime[10];    
-  snprintf(nameEtaprime,10, "Etaprime");    
+  snprintf(nameEtaprime,10,"Etaprime");    
   genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
   AddSource2Generator(nameEtaprime,genetaprime);
   TF1 *fPtEtaprime = genetaprime->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
+  fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
-// phi  
-  AliGenParam * genphi=0;
+
+  // phi  
+  if(fHeaviestParticle<kGenPhi)return;
+  AliGenParam *genphi=0;
   Char_t namePhi[10];    
-  snprintf(namePhi, 10, "Phi");    
+  snprintf(namePhi,10,"Phi");    
   genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
   AddSource2Generator(namePhi,genphi);
   TF1 *fPtPhi = genphi->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
+  fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+
+  // jpsi  
+  if(fHeaviestParticle<kGenJpsi)return;
+  AliGenParam *genjpsi=0;
+  Char_t nameJpsi[10];    
+  snprintf(nameJpsi,10,"Jpsi");    
+  genjpsi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
+  AddSource2Generator(nameJpsi,genjpsi);
+  TF1 *fPtJpsi = genjpsi->GetPt();
+#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
+  fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
+#else
+  fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+#endif
+
 }
 
 //-------------------------------------------------------------------
@@ -173,6 +207,7 @@ void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
   genSource->SetYRange(fYMin, fYMax);
   genSource->SetPhiRange(phiMin, phiMax);
   genSource->SetWeighting(fWeightingMode);
+  genSource->SetForceGammaConversion(fForceConv);
   if (!gMC) genSource->SetDecayer(fDecayer);  
   genSource->Init();
     
@@ -182,7 +217,7 @@ void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
 //-------------------------------------------------------------------
 void AliGenEMCocktail::Init()
 {
-// Initialisation
+  // Initialisation
   TIter next(fEntries);
   AliGenCocktailEntry *entry;
   if (fStack) {
@@ -195,7 +230,7 @@ void AliGenEMCocktail::Init()
 //_________________________________________________________________________
 void AliGenEMCocktail::Generate()
 {
-// Generate event 
+  // Generate event 
   TIter next(fEntries);
   AliGenCocktailEntry *entry = 0;
   AliGenCocktailEntry *preventry = 0;
@@ -206,16 +241,16 @@ void AliGenEMCocktail::Generate()
 
   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
     
-// Generate the vertex position used by all generators    
+  // Generate the vertex position used by all generators    
   if(fVertexSmear == kPerEvent) Vertex();
 
-//Reseting stack
+  //Reseting stack
   AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader)
     if (runloader->Stack())
       runloader->Stack()->Clean();
   
-// Loop over generators and generate events
+  // Loop over generators and generate events
   Int_t igen = 0;
   Float_t evPlane;
   Rndm(&evPlane,1);
@@ -223,7 +258,6 @@ void AliGenEMCocktail::Generate()
   while((entry = (AliGenCocktailEntry*)next())) {
     gen = entry->Generator();
     gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
-    gen->SetTime(fTime);
     
     if (fNPart > 0) {
       igen++;  
@@ -238,7 +272,7 @@ void AliGenEMCocktail::Generate()
   }  
   next.Reset();
 
-// Setting weights for proper absolute normalization
+  // Setting weights for proper absolute normalization
   Int_t iPart, iMother;
   Int_t pdgMother = 0;
   Double_t weight = 0.;
@@ -273,11 +307,13 @@ void AliGenEMCocktail::Generate()
     case 333:
       dNdy = fYieldArray[kGenPhi];
       break;
+    case 443:
+      dNdy = fYieldArray[kGenJpsi];
+      break;
       
     default:
       dNdy = 0.;
     }
-    
     weight = dNdy*part->GetWeight();
     part->SetWeight(weight);
   }    
@@ -290,9 +326,6 @@ void AliGenEMCocktail::Generate()
   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
   
   fHeader->SetPrimaryVertex(eventVertex);
-  fHeader->SetInteractionTime(fTime);
 
   gAlice->SetGenEventHeader(fHeader);
 }
-
-    
index db823b9..947506c 100644 (file)
 //
 
 #include "AliGenCocktail.h"
-#include "AliGenParam.h"
-#include "AliDecayer.h"
+#include "AliGenEMlib.h"
 #include "AliPythia.h"
+#include "AliDecayer.h"
+#include "AliGenParam.h"
+
 
 class AliGenCocktailEntry;
 
 class AliGenEMCocktail : public AliGenCocktail
 {
- public:
+public:
 
     AliGenEMCocktail();
-    enum GeneratorCode { kGenPizero, kGenEta, kGenRho, kGenOmega, kGenEtaprime, kGenPhi, kGENs };    
+  enum GeneratorCode { kGenPizero=0, kGenEta, kGenRho, kGenOmega, kGenEtaprime, kGenPhi, kGenJpsi, kGENs };    
 
     virtual ~AliGenEMCocktail();    
     virtual void Init();
@@ -29,28 +31,40 @@ class AliGenEMCocktail : public AliGenCocktail
     virtual void Generate();    
     Float_t GetDecayMode()         const {return fDecayMode;}
     Float_t GetWeightingMode()     const {return fWeightingMode;}
-    
     void    SetDecayer(AliDecayer* const decayer){fDecayer = decayer;}
     void    SetDecayMode(Decay_t decay){ fDecayMode = decay;}
     void    SetWeightingMode(Weighting_t weight){ fWeightingMode = weight;}
-    void    SetNPart(Int_t npart){ fNPart = npart;}
+  void    SetNPart(Int_t npart){ fNPart = npart; }
+  void    SetPtParam(Int_t PtSelect){ fPtSelect = PtSelect; }
+  void    SetCentrality(Int_t cent){ fCentrality = cent; }
+  void    SetV2Systematic(Int_t v2sys){ fV2Systematic = v2sys; }
+  void    SetForceGammaConversion(Bool_t force=kTRUE){ fForceConv=force; }
+  void    SetHeaviestParticle(Int_t part){ fHeaviestParticle=part; }
     
- protected:
+protected:
 
-    //
- private:
+  TVector3& OrthogonalVector(TVector3 &inVec);
+    
+private:
     AliGenEMCocktail(const AliGenEMCocktail &cocktail); 
     AliGenEMCocktail & operator=(const AliGenEMCocktail &cocktail); 
 
     void AddSource2Generator(Char_t *nameReso, AliGenParam* const genReso);
     
     AliDecayer* fDecayer;        // External decayer
-    Decay_t fDecayMode;          // decay mode in which resonances are forced to decay, default: kAll
-    Weighting_t fWeightingMode;  // weighting mode: kAnalog or kNonAnalog
+  Decay_t fDecayMode;   //decay mode in which resonances are forced to decay, default: kAll
+  Weighting_t fWeightingMode; //weighting mode: kAnalog or kNonAnalog
     
     Int_t    fNPart;             // multiplicity of each source per event
     Double_t fYieldArray[kGENs]; // array of dN/dy for each source
 
+  Int_t fPtSelect; // selected pT parameterization
+  Int_t fCentrality; // selected centrality
+  Int_t fV2Systematic; //selected systematic error for v2 parameters
+
+  Bool_t fForceConv; //select whether you want to force all gammas to convert imidediately
+  Int_t fHeaviestParticle; //select up to which particle to simulate
+
     ClassDef(AliGenEMCocktail,1)  //  cocktail for EM physics
 };
 
index 1851e87..5be02c9 100644 (file)
@@ -26,6 +26,7 @@
 /////////////////////////////////////////////////////////////////////////////
 
 
+#include <Riostream.h>
 #include "TMath.h"
 #include "TRandom.h"
 #include "TString.h"
 
 ClassImp(AliGenEMlib)
 
-AliGenEMlib::Centbins_t AliGenEMlib::fCentbin=k2030;
+//Initializers for static members
+Int_t AliGenEMlib::fgSelectedPtParam=AliGenEMlib::kPizero7TeVpp;
+Int_t AliGenEMlib::fgSelectedCentrality=AliGenEMlib::kpp;
+Int_t AliGenEMlib::fgSelectedV2Systematic=AliGenEMlib::kNoV2Sys;
 
 Double_t AliGenEMlib::CrossOverLc(const double a, const double b, const double x){
-if(x<b-a/2) return 1.0;
+  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;
 }
@@ -45,38 +49,43 @@ Double_t AliGenEMlib::CrossOverRc(const double a, const double b, const double x
   return 1-CrossOverLc(a,b,x);
 }
 
-const Double_t AliGenEMlib::fpTparam[8][14] = {
-  // charged pion
-  { 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 0.0000000000, 1.0000000000, 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 9999.0000000, 1.0000000000, 0.0000000000, 0.0000000000 },   //
-  { 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 0.0000000000, 1.0000000000, 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 9999.0000000, 1.0000000000, 0.0000000000, 0.0000000000 },   //
-  { 3.217746e+03, 1.632570e+00, 9.340162e+00, 1.0000, 2.650000e+00, 4.200635e+00, 9.039589e+08, 1.191548e-01, 6.907293e+00, 1.0000, 6.750000e+00, 4.099970e+00, 6.036772e+01, 5.928279e+00 },   //10-20
-  { 2.125480e+03, 1.711882e+00, 9.585665e+00, 1.0000, 3.100000e+00, 5.041511e+00, 2.431808e+08, 1.155071e-01, 6.574663e+00, 1.0000, 6.250000e+00, 1.842070e+00, 3.928902e+01, 5.860970e+00 },   //20-30
-  { 1.577897e+03, 1.411948e+00, 8.638815e+00, 1.0000, 2.550000e+00, 4.066432e+00, 3.774439e+08, 1.000330e-01, 6.535971e+00, 1.0000, 6.750000e+00, 2.482514e+00, 3.495494e+01, 5.954321e+00 },   //30-40
-  { 7.061859e+02, 1.223810e+00, 8.532113e+00, 1.0000, 1.350000e+00, 1.956213e+00, 1.318169e+04, 5.658401e-01, 7.157575e+00, 1.0000, 5.250000e+00, 3.900000e+00, 1.958841e+01, 6.114787e+00 },   //40-50
-  { 7.061859e+02, 1.223810e+00, 8.532113e+00, 1.0000, 1.350000e+00, 1.956213e+00, 1.318169e+04, 5.658401e-01, 7.157575e+00, 1.0000, 5.250000e+00, 3.900000e+00, 1.958841e+01, 6.114787e+00 },   //
-  { 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 0.0000000000, 1.0000000000, 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 9999.0000000, 1.0000000000, 0.0000000000, 0.0000000000 } }; //pp
-
-const Double_t AliGenEMlib::fv2param[2][8][9] = { { 
-    // charged pion
-    {  3.971545e-02, 2.111261e+00, 6.499885e-02, 3.659836e+00, 7.812234e+00, 1.479471e-02, 1.241708e+00, 2.817950e+00, 1.910663e-01 }   //0-5
-    ,{ 7.424082e-02, 2.417588e+00, 7.891084e-02, 2.232841e+00, 3.398147e+00, 3.410206e-02, 4.138276e-01,-1.152430e+00, 5.729093e-01 }   //5-10
-    ,{ 1.094608e-01, 2.357420e+00, 7.614904e-02, 2.294245e+00, 3.538364e+00, 4.932739e-02, 4.557926e-01, 9.218702e-03, 6.428540e-01 }   //10-20
-    ,{ 1.456377e-01, 2.322408e+00, 7.747166e-02, 2.148424e+00, 3.238113e+00, 5.414722e-02, 4.042938e-01, 1.903040e-01, 6.816790e-01 }   //20-30
-    ,{ 1.745154e-01, 2.234489e+00, 7.962225e-02, 2.007075e+00, 2.925536e+00, 5.499138e-02, 3.958957e-01, 1.780793e+00, 4.852930e-01 }   //30-40
-    ,{ 2.384302e-01, 1.578935e+00, 9.234643e-02, 4.328926e+00, 1.020669e+01, 1.001340e-01, 2.433905e+00, 2.966673e+00, 2.646807e-01 }   //40-50
-    ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }   //
-    ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }   //pp no v2
-  },{
-    // charged kaon
-    {  0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
-    ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
-    ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
-    ,{ 1.432314e-01, 1.607297e+00, 2.296868e-01, 1.821156e+00, 2.951652e+00,-2.201385e-01, 1.110419e-01, 8.228701e+00, 4.469488e-01 }  //20-30
-    ,{ 1.691586e-01, 1.617165e+00, 2.159350e-01, 1.649338e+00, 2.607635e+00,-9.536279e+00, 3.860086e-01, 1.802996e+01, 2.509780e-01 }  //30-40
-    ,{ 1.733831e-01, 1.712705e+00, 1.935862e-01, 1.745523e+00, 2.845436e+00,-5.772356e-01, 6.861616e-02, 5.292878e+00, 9.058815e-01 }  //40-50
-    ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
-    ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //pp no v2
-  } };
+const Double_t AliGenEMlib::fgkV2param[16][15] = {
+  // charged pion                                                                                                                        cent, based on
+  {  0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000 }   // pp no V2
+  ,{ 6.792545e-02, 1.324227e+00, 2.071866e-02, 3.385157e+00, 5.183718e+00, 2.783082e-02, 5.000000e+00, 4.078812e+00, 5.444857e-02, 1.555569e-03, 1.099253e+00, 0, 1, 8.823293e-04, 1.310387e+00 }   // 0-5,  Francesco
+  ,{ 1.219655e-01, 1.237698e+00, 2.595132e-02, 3.157085e+00, 3.697252e+00, 4.576262e-02, 1.103851e+00, 5.102592e+00, 1.072299e-01, 2.338602e-03, 8.475128e-01, 0, 1, 1.631554e-03, 9.432292e-01 }   // 5-10, Francesco
+  ,{ 1.758981e-01, 1.253238e+00, 2.823197e-02, 3.455476e+00, 4.283204e+00, 6.053191e-02, 6.737159e-01, 4.420205e+00, 1.733567e-01, 3.292614e-03, 6.389341e-01, 0, 1, 3.234254e-03, 5.304084e-01 }   // 10-20,Francesco
+  ,{ 2.231108e-01, 1.316781e+00, 3.537289e-02, 3.228243e+00, 3.697226e+00, 7.131631e-02, 5.902192e-01, 4.938186e+00, 2.123481e-01, 4.775233e-03, 5.707533e-01, 0, 1, 3.906720e-03, 5.261247e-01 }   // 20-30,Francesco
+  ,{ 2.045341e-01, 1.324054e+00, 1.125395e-01, 2.179936e+00, 5.614603e+00, 8.946014e-02, 5.601938e-01, 2.959228e+00, 3.669442e-01, 5.318697e-03, 6.032741e-01, 0, 1, 4.231316e-03, 5.563485e-01 }   // 30-40,Francesco
+  ,{ 2.074778e-01, 1.894322e+00, 6.019224e-02, 2.414500e+00, 3.822721e+00, 9.891310e-02, 4.976126e-01, 1.256500e+00, 5.210237e-01, 5.687290e-03, 6.785355e-01, 0, 1, 3.954051e-03, 6.913428e-01 }   // 40-50,Francesco
+  ,{ 1.961099e-01, 2.122719e+00, 7.633636e-02, 1.731213e+00, 2.427536e+00, 1.026013e-01, 7.932261e-01, 4.082226e+00, 2.495494e-01, 5.677386e-03, 8.405024e-01, 0, 1, 2.961753e-03, 1.008892e+00 }   // 50-60,Francesco
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000 }   // 0-10
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000 }   // 20-40
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000 }   // 40-60
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000 }   // 60-80
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000 }   // 0-20
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000 }   // 0-40
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000 }   // 20-80
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000 }   // 40-80
+};
+
+// MASS   0=>PIZERO, 1=>ETA, 2=>RHO, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI, 6=>JPSI
+const Double_t AliGenEMlib::fgkHM[7] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946, 3.0969};
+
+const Double_t AliGenEMlib::fgkMtFactor[2][7] = { 
+  // {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
+  //J/Psi PbPb from Comparison with Julian Books J/Psi -> e+e-, might be contradicting with https://aliceinfo.cern.ch/Figure/node/3457
+  //best guess:
+  {1., 0.4, 1.0, 0.8, 0.4, 0.2, 0.0054}, //pp
+  {1., 0.4, 1.0, 0.8, 0.4, 0.25, 0.01}  //PbPb
+};
 
 //==========================================================================
 //
@@ -86,36 +95,236 @@ const Double_t AliGenEMlib::fv2param[2][8][9] = { {
 
 //--------------------------------------------------------------------------
 //
+//                              General functions
+//
+//--------------------------------------------------------------------------
+Double_t AliGenEMlib::PtModifiedHagedornThermal(const Double_t pt,
+                                                const Double_t c,
+                                                const Double_t p0,
+                                                const Double_t p1,
+                                                const Double_t n,
+                                                const Double_t cT,
+                                                const 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);
+
+  return invYield*(2*TMath::Pi()*pt);
+}
+
+
+
+Double_t AliGenEMlib::PtModifiedHagedornExp(const Double_t pt,
+                                           const Double_t c,
+                                           const Double_t p1,
+                                           const Double_t p2,
+                                           const Double_t p0,
+                                           const 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);
+
+  return invYield*(2*TMath::Pi()*pt);
+}
+
+
+Double_t AliGenEMlib::PtModifiedHagedornExp2(const Double_t pt,
+                                             const Double_t c,
+                                             const Double_t a,
+                                             const Double_t b,
+                                             const Double_t p0,
+                                             const Double_t p1,
+                                             const Double_t d,
+                                             const 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);
+
+  return invYield*(2*TMath::Pi()*pt);
+}
+
+Double_t AliGenEMlib::PtTsallis(const Double_t pt,
+                                const Double_t m,
+                                const Double_t c,
+                                const Double_t T,
+                                const 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);
+
+  return invYield*(2*TMath::Pi()*pt);
+}
+
+
+
+//--------------------------------------------------------------------------
+//
 //                              Pizero
 //
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpPizero(TRandom *)
 {
-// Return pizero pdg code
+  // Return pizero pdg code
   return 111;     
 }
 
 Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
 {
-// Generate pizero pT distribution from modified Hagedorn parameterization
-// taken from fit to unidentified hadrons in pp at 7 TeV
-  // const Double_t kc=0.000565;
-  // const Double_t kp0=0.2472;
-  // const Double_t kp1=4.354;
-  // const Double_t kn=7.007;
-  Double_t invYield;
-  Double_t x=*px;
-
-  // invYield = kc/TMath::Power(kp0+x/kp1,kn);
-
-  if(!x)return 0;
-  const Double_t *c=fpTparam[fCentbin];
+  // 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(fgSelectedPtParam|fgSelectedCentrality) {
+    // fit to pi charged v1
+    // charged pion from ToF, unidentified hadrons scaled with pion from TPC
+    // for Pb-Pb @ 2.76 TeV
+  case kPichargedPbPb|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 kPichargedPbPb|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 kPichargedPbPb|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 kPichargedPbPb|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 kPichargedPbPb|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 kPichargedPbPb|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 kPichargedPbPb|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 kPichargedPbPb|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 kPichargedPbPb|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 kPichargedPbPb|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 kPichargedPbPb|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;
    
-  invYield = CrossOverLc(c[5],c[4],x)*c[0]/TMath::Power(c[3]+x/c[1],c[2]);
-  invYield +=CrossOverRc(c[5],c[4],x)*c[6]/TMath::Power(c[9]+x/c[7],c[8])*CrossOverLc(c[11],c[10],x);
-  invYield +=CrossOverRc(c[11],c[10],x)*c[12]/TMath::Power(x,c[13]);
 
-  return invYield*(2*TMath::Pi()*x);
+    // 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_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
@@ -126,9 +335,39 @@ Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
 
 Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return V2Param(px,fv2param[0][fCentbin]);
+  int selectedCentrality;
+  double n1,n2,v1,v2;
+  switch(fgSelectedCentrality) {
+  case k0010:
+    selectedCentrality=fgSelectedCentrality;
+    fgSelectedCentrality=k0005;
+    n1=PtPizero(px,(Double_t*) 0);
+    v1=V2Param(px,fgkV2param[fgSelectedCentrality]);
+    fgSelectedCentrality=k0510;
+    n2=PtPizero(px,(Double_t*) 0);
+    v2=V2Param(px,fgkV2param[fgSelectedCentrality]);
+    fgSelectedCentrality=selectedCentrality;
+    return (n1*v1+n2*v2)/(n1+n2);
+    break;
+  case k2040:
+    selectedCentrality=fgSelectedCentrality;
+    fgSelectedCentrality=k2030;
+    n1=PtPizero(px,(Double_t*) 0);
+    v1=V2Param(px,fgkV2param[fgSelectedCentrality]);
+    fgSelectedCentrality=k3040;
+    n2=PtPizero(px,(Double_t*) 0);
+    v2=V2Param(px,fgkV2param[fgSelectedCentrality]);
+    fgSelectedCentrality=selectedCentrality;
+    return (n1*v1+n2*v2)/(n1+n2);
+    break;
+
+  default:
+    return V2Param(px,fgkV2param[fgSelectedCentrality]);
+  }
 }
 
+
+
 //--------------------------------------------------------------------------
 //
 //                              Eta
@@ -136,14 +375,56 @@ Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpEta(TRandom *)
 {
-// Return eta pdg code
+  // Return eta pdg code
   return 221;     
 }
 
 Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
 {
-// Eta pT
+
+  // 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;
+
+  }
+
 }
 
 Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
@@ -153,7 +434,7 @@ Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
 
 Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return V2Param(px,fv2param[1][fCentbin]);
+  return KEtScal(*px,1); //V2Param(px,fgkV2param[1][fgSelectedV2Param]);
 }
 
 //--------------------------------------------------------------------------
@@ -163,20 +444,24 @@ Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpRho(TRandom *)
 {
-// Return rho pdg code
+  // Return rho pdg code
   return 113;     
 }
 
 Double_t AliGenEMlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
 {
-// Rho pT
+  // Rho pT
   return MtScal(*px,2);
 }
 
 Double_t AliGenEMlib::YRho( const Double_t *py, const Double_t */*dummy*/ )
 {
   return YFlat(*py);
+}
 
+Double_t AliGenEMlib::V2Rho( const Double_t *px, const Double_t */*dummy*/ )
+{
+  return KEtScal(*px,2);
 }
 
 //--------------------------------------------------------------------------
@@ -186,22 +471,27 @@ Double_t AliGenEMlib::YRho( const Double_t *py, const Double_t */*dummy*/ )
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpOmega(TRandom *)
 {
-// Return omega pdg code
+  // Return omega pdg code
   return 223;     
 }
 
 Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
 {
-// Omega pT
+  // Omega pT
   return MtScal(*px,3);
 }
 
 Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
 {
   return YFlat(*py);
+}
 
+Double_t AliGenEMlib::V2Omega( const Double_t *px, const Double_t */*dummy*/ )
+{
+  return KEtScal(*px,3);
 }
 
+
 //--------------------------------------------------------------------------
 //
 //                              Etaprime
@@ -209,13 +499,13 @@ Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpEtaprime(TRandom *)
 {
-// Return etaprime pdg code
+  // Return etaprime pdg code
   return 331;     
 }
 
 Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
 {
-// Eta pT
+  // Eta pT
   return MtScal(*px,4);
 }
 
@@ -225,6 +515,11 @@ Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
 
 }
 
+Double_t AliGenEMlib::V2Etaprime( const Double_t *px, const Double_t */*dummy*/ )
+{
+  return KEtScal(*px,4);
+}
+
 //--------------------------------------------------------------------------
 //
 //                              Phi
@@ -232,29 +527,71 @@ Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
 //--------------------------------------------------------------------------
 Int_t AliGenEMlib::IpPhi(TRandom *)
 {
-// Return phi pdg code
+  // Return phi pdg code
   return 333;     
 }
 
 Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
 {
-// Phi pT
+  // Phi pT
   return MtScal(*px,5);
 }
 
 Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
 {
   return YFlat(*py);
-
 }
 
-Double_t AliGenEMlib::YFlat(Double_t /*y*/)
+Double_t AliGenEMlib::V2Phi( const Double_t *px, const Double_t */*dummy*/ )
 {
+  return KEtScal(*px,5);
+}
+
 //--------------------------------------------------------------------------
 //
-//                    flat rapidity distribution 
+//                              Jpsi
 //
 //--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpJpsi(TRandom *)
+{
+  // Return phi pdg code
+  return 443;
+}
+
+Double_t AliGenEMlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/ )
+{
+  // Jpsi pT
+  return MtScal(*px,6);
+}
+
+Double_t AliGenEMlib::YJpsi( const Double_t *py, const Double_t */*dummy*/ )
+{
+  return YFlat(*py);
+}
+
+Double_t AliGenEMlib::V2Jpsi( const Double_t *px, const Double_t */*dummy*/ )
+{
+  const int oldSys=fgSelectedV2Systematic;
+  fgSelectedV2Systematic=kNoV2Sys;
+  double ret=0;
+
+  switch(oldSys){
+  case kLoV2Sys: ret=0; break;
+  case kNoV2Sys: ret=KEtScal(*px,6)/2; break;
+  case kUpV2Sys: ret=KEtScal(*px,6); break;
+  }
+
+  fgSelectedV2Systematic=oldSys;
+  return ret;
+}
+
+Double_t AliGenEMlib::YFlat(Double_t /*y*/)
+{
+  //--------------------------------------------------------------------------
+  //
+  //                    flat rapidity distribution 
+  //
+  //--------------------------------------------------------------------------
 
   Double_t dNdy = 1.;   
 
@@ -268,35 +605,42 @@ Double_t AliGenEMlib::YFlat(Double_t /*y*/)
 //
 //============================================================= 
 //
- Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
+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. 
-
-// MASS   0=>PIZERO, 1=>ETA, 2=>RHO, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI
+  // Function for the calculation of the Pt distribution for a 
+  // given particle np, from the pizero Pt distribution using  
+  // mt scaling. 
 
-  const Double_t khm[6] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946};
-
-  Double_t scaledPt = sqrt(pt*pt + khm[np]*khm[np] - khm[0]*khm[0]);
+  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
-
-  Double_t normPt = 5.;
-  Double_t scaledNormPt = sqrt(normPt*normPt + khm[np]*khm[np] - khm[0]*khm[0]);
-  const Double_t kfmax[6]={1., 0.48, 1.0, 0.9, 0.25, 0.4};
+  //     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 = kfmax[np] * (PtPizero(&normPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
+  Double_t norm = fgkMtFactor[int(bool(fgSelectedCentrality))][np] * (PtPizero(&NormPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
 
   return norm*(pt/scaledPt)*scaledYield;
 }
 
+Double_t AliGenEMlib::KEtScal(const 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/nq*(sqrt(pt*pt+fgkHM[np]*fgkHM[np])-fgkHM[np])+fgkHM[0],2)-fgkHM[0]*fgkHM[0]);
+  double val=V2Pizero(&scaledPt, (Double_t*) 0);
+  static const double syserr[12]={0., 0.09, 0.07, 0.06, 0.04, 0.04, 0.04, 0.05, 0., 0., 0., 0.}; //based on pi vs kaon
+  double sys=fgSelectedV2Systematic*min(fgkV2param[fgSelectedCentrality][0],fgkV2param[fgSelectedCentrality][8])*syserr[fgSelectedCentrality];
+  return TMath::Max(val+sys,0.0);
+}
+
 Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
 {
   // Very general parametrization of the v2
 
-  return TMath::Max(CrossOverLc(par[4],par[3],px[0])*(2*par[0]/(1+TMath::Exp(par[1]*(par[2]-px[0])))-par[0])+CrossOverRc(par[4],par[3],px[0])*((par[8]-par[5])/(1+TMath::Exp(par[6]*(px[0]-par[7])))+par[5]),0.0);
+  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=fgSelectedV2Systematic*par[11+fgSelectedV2Systematic*2]*pow(pt,par[12+fgSelectedV2Systematic*2]);
+  return TMath::Max(val+sys,0.0);
 }
 
 Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
@@ -306,6 +650,8 @@ Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
   return 1.0;
 }
 
+
+
 //==========================================================================
 //
 //                     Set Getters 
@@ -318,7 +664,7 @@ typedef Int_t (*GenFuncIp) (TRandom *);
 
 GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
 {
-// Return pointer to pT parameterisation
+  // Return pointer to pT parameterisation
    GenFunc func=0;
    TString sname(tname);
 
@@ -342,6 +688,9 @@ GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
     case kPhi:
       func=PtPhi;
       break;
+    case kJpsi:
+      func=PtJpsi;
+      break;
 
     default:
       func=0;
@@ -352,7 +701,7 @@ GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
 
 GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
 {
-// Return pointer to y- parameterisation
+  // Return pointer to y- parameterisation
    GenFunc func=0;
    TString sname(tname);
 
@@ -376,6 +725,9 @@ GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
     case kPhi:
          func=YPhi;
          break;
+    case kJpsi:
+      func=YJpsi;
+      break;
 
     default:
         func=0;
@@ -386,7 +738,7 @@ GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
 
 GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
 {
-// Return pointer to particle type parameterisation
+  // Return pointer to particle type parameterisation
    GenFuncIp func=0;
    TString sname(tname);
 
@@ -410,6 +762,9 @@ GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
     case kPhi:
          func=IpPhi;
          break;
+    case kJpsi:
+      func=IpJpsi;
+      break;
 
     default:
         func=0;
@@ -444,6 +799,9 @@ GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
     case kPhi:
       func=V2Pizero;
       break;
+    case kJpsi:
+      func=V2Jpsi;
+      break;
 
     default:
       func=0;
@@ -451,3 +809,4 @@ GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
     }
   return func;
 }
+
index de51364..81c0864 100644 (file)
@@ -20,63 +20,127 @@ class TRandom;
 
 class AliGenEMlib :public AliGenLib {
 public:
+    
+  enum particles{kPizero, kEta, kRho, kOmega, kEtaprime, kPhi, kJpsi};
+  enum centrality{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{kPizero7TeVpp=0x10, kPizeroEta7TeVpp=0x20, kPizero7TeVpplow=0x30, kPizeroEta7TeVpplow=0x40, kPizero7TeVpphigh=0x50, kPizeroEta7TeVpphigh=0x60, kPizero2760GeVpp=0x70, kPizeroEta2760GeVpp=0x80, kPizero2760GeVpplow=0x90, kPizeroEta2760GeVpplow=0xA0, kPizero2760GeVpphigh=0xB0, kPizeroEta2760GeVpphigh=0xC0, kPichargedPbPb=0xD0, kPizeroPbPb=0xE0, kPichargedPPb=0xF0 };
+  enum v2Sys{kLoV2Sys=-1, kNoV2Sys=0, kUpV2Sys=+1};
+  
+  AliGenEMlib() { } ;
+
+  static void SelectParams(Int_t ptSelect, Int_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;
 
-  typedef enum {k0005=0, k0510, k1020, k2030, k3040, k4050, k5060, kpp} Centbins_t;
-  enum particles{kPizero, kEta, kRho, kOmega, kEtaprime, kPhi};
+  //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(const Double_t pt, 
+                                           const Double_t c, 
+                                           const Double_t p0, 
+                                           const Double_t p1, 
+                                           const Double_t n,
+                                           const Double_t cT,
+                                           const Double_t T);
+
+
+  static Double_t PtModifiedHagedornExp(const Double_t pt,
+                                       const Double_t c,
+                                       const Double_t p0,
+                                       const Double_t p1,
+                                       const Double_t p2,
+                                       const Double_t n); 
+
+
+  static Double_t PtModifiedHagedornExp2(const Double_t pt,
+                                           const Double_t c,
+                                           const Double_t a,
+                                           const Double_t b,
+                                           const Double_t p0,
+                                           const Double_t p1,
+                                           const Double_t d,
+                                           const Double_t n);
+
+
+  static Double_t PtTsallis(const Double_t pt,
+                           const Double_t m,
+                           const Double_t c,
+                           const Double_t T,
+                           const Double_t n);
+
+
 
-private:
 
   // Pizero
     static Int_t    IpPizero(TRandom *ran);
-    static Double_t PtPizero( const Double_t *px, const Double_t *dummy );
+  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 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 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 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 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 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 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 CrossOverLc(const double a, const double b, const double x);
   static Double_t CrossOverRc(const double a, const double b, const double x);
 
-  static const Double_t fv2param[2][8][9];
-  static const Double_t fpTparam[8][14];
-  static Centbins_t fCentbin;
+  static const Double_t fgkV2param[16][15]; // array of v2 parameter
+  static const Double_t fgkHM[7];
+  static const Double_t fgkMtFactor[2][7];
 
   ClassDef(AliGenEMlib,0)
 };
 
+
 #endif
index 6589d1a..be599f9 100644 (file)
@@ -68,7 +68,8 @@ AliGenParam::AliGenParam()
        fTrials(0),
        fDeltaPt(0.01),
        fSelectAll(kFALSE),
-       fDecayer(0)
+  fDecayer(0),
+  fForceConv(kFALSE)
 {
   // Default constructor
 }
@@ -91,7 +92,8 @@ AliGenParam::AliGenParam(Int_t npart, const AliGenLib * Library,  Int_t param, c
      fTrials(0),
      fDeltaPt(0.01),
      fSelectAll(kFALSE),
-     fDecayer(0)
+   fDecayer(0),
+   fForceConv(kFALSE)
 {
   // Constructor using number of particles parameterisation id and library
     fName = "Param";
@@ -118,7 +120,8 @@ AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char
     fTrials(0),
     fDeltaPt(0.01),
     fSelectAll(kFALSE),
-    fDecayer(0)
+  fDecayer(0),
+  fForceConv(kFALSE)
 {
   // Constructor using parameterisation id and number of particles
   //
@@ -166,7 +169,8 @@ AliGenParam::AliGenParam(Int_t npart, Int_t param,
      fTrials(0),
      fDeltaPt(0.01),
      fSelectAll(kFALSE),
-     fDecayer(0)
+  fDecayer(0),
+  fForceConv(kFALSE)
 {
   // Constructor
   // Gines Martinez 1/10/99 
@@ -194,6 +198,120 @@ AliGenParam::~AliGenParam()
   delete  fdNdPhi;
 }
 
+//-------------------------------------------------------------------
+TVector3 AliGenParam::OrthogonalVector(TVector3 &inVec){
+  double abc[]={inVec.x(), inVec.y(), inVec.z()}; 
+  double xyz[]={1,1,1};
+  int solvDim=0;
+  double tmp=abc[0];
+  for(int i=0; i<3; i++)
+    if(abs(abc[i])>tmp){
+      solvDim=i;
+      tmp=abs(abc[i]);
+    }
+  xyz[solvDim]=(-abc[(1+solvDim)%3]-abc[(2+solvDim)%3])/abc[(0+solvDim)%3];
+  
+  TVector3 res(xyz[0],xyz[1],xyz[2]);
+  return res;
+}
+
+double AliGenParam::ScreenFunc1(double d){
+  if(d>1)
+    return 21.12-4.184*log(d+0.952);
+  else
+    return 20.867-3.242*d+0.652*d*d;
+}
+
+double AliGenParam::ScreenFunc2(double d){
+  if(d>1)
+    return 21.12-4.184*log(d+0.952);
+  else
+    return 20.209-1.93*d-0.086*d*d;
+}
+
+double AliGenParam::EnergyFraction(double Z, double E){
+  double e0=0.000511/E;
+  double aZ=Z/137.036;
+  double dmin=ScreenVar(Z,e0,0.5);
+  double fcZ=(aZ*aZ)*(1/(1+aZ*aZ)+0.20206-0.0368*aZ*aZ+0.0083*aZ*aZ*aZ);
+  double Fz=8*log(Z)/3;
+  if(E>0.05)
+    Fz+=8*fcZ;
+  double dmax=exp((42.24-Fz)/8.368)-0.952;
+  if(!(dmax>dmin))
+    return fRandom->Uniform(e0,0.5);
+
+  double e1=0.5-0.5*sqrt(1-dmin/dmax);
+  double emin=TMath::Max(e0,e1);
+  
+  double normval=1/(0.5*(ScreenFunc1(dmin)-0.5*Fz)+0.1666667*(ScreenFunc2(dmin)-0.5*Fz));
+  while(true){
+    double y=fRandom->Uniform(emin,1);
+    double eps=y*(0.38453+y*(0.10234+y*(0.026072+y*(0.014367-0.027313*y)))); //inverse to the enveloping cumulative probability distribution
+    double val=fRandom->Uniform(0,2.12*(eps*eps-eps)+1.53); //enveloping probability density
+    double d=ScreenVar(Z,e0,eps);
+    double bh=((eps*eps)+(1-eps)*(1-eps))*(ScreenFunc1(d)-0.5*Fz)+0.6666667*eps*(1-eps)*(ScreenFunc2(d)-0.5*Fz);
+    bh*=normval;
+    if(val<bh)
+      return eps;
+  }
+}
+
+double AliGenParam::PolarAngle(double E){
+  float rand[3];
+  AliRndm rndm;
+  rndm.Rndm(rand,3);
+  double u=-8*log(rand[1]*rand[2])/5;
+  if(!(rand[0]<9.0/36))
+    u/=3;
+  return u*0.000511/E;
+}
+
+Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
+{
+  //based on: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node27.html
+  //     and: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node58.html
+  Int_t nPartNew=nPart;
+  for(int iPart=0; iPart<nPart; iPart++){      
+    TParticle *gamma = (TParticle *) particles->At(iPart);
+    if(gamma->GetPdgCode()!=22) continue;
+
+    TVector3 gammaV3(gamma->Px(),gamma->Py(),gamma->Pz());
+    Float_t az=fRandom->Uniform(TMath::Pi()*2);
+    double frac=EnergyFraction(1,gamma->Energy());
+    double Ee1=frac*gamma->Energy();
+    double Ee2=(1-frac)*gamma->Energy();
+    double Pe1=Ee1;//sqrt(Ee1*Ee1-0.000511*0.000511);
+    double Pe2=Ee2;//sqrt(Ee2*Ee2-0.000511*0.000511);
+    
+    TVector3 rotAxis(OrthogonalVector(gammaV3));
+    rotAxis.Rotate(az,gammaV3);
+    TVector3 e1V3(gammaV3);
+    e1V3.Rotate(PolarAngle(Ee1),rotAxis);
+    e1V3=e1V3.Unit();
+    e1V3*=Pe1;
+    TVector3 e2V3(gammaV3);
+    e2V3.Rotate(-PolarAngle(Ee2),rotAxis);
+    e2V3=e2V3.Unit();
+    e2V3*=Pe2;
+    // gamma = new TParticle(*gamma);
+    // particles->RemoveAt(iPart);
+    gamma->SetFirstDaughter(nPartNew+1);
+    gamma->SetLastDaughter(nPartNew+2);
+    // new((*particles)[iPart]) TParticle(*gamma);
+    // delete gamma;
+
+    TLorentzVector vtx;
+    gamma->ProductionVertex(vtx);
+    new((*particles)[nPartNew]) TParticle(11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e1V3,Ee1), vtx);
+    nPartNew++;
+    new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e2V3,Ee2), vtx);
+    nPartNew++;
+  }
+  // particles->Compress();
+  return particles->GetEntriesFast();
+}
+
 //____________________________________________________________
 void AliGenParam::Init()
 {
@@ -370,14 +488,14 @@ void AliGenParam::GenerateN(Int_t ntimes)
          }
       //
       // phi
-      // if(!ipa)
-      //       phi=fEvPlane; //align first particle of each event with event plane
-      // else{
+         //      if(!ipa)
+         //phi=fEvPlane; //align first particle of each event with event plane
+         //else{
        double v2 = fV2Para->Eval(pt);
        fdNdPhi->SetParameter(0,v2);
        fdNdPhi->SetParameter(1,fEvPlane);
        phi=fdNdPhi->GetRandom();
-      // }
+       //     }
          
          pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
          theta=TMath::ATan2(pt,pl);
@@ -422,6 +540,14 @@ void AliGenParam::GenerateN(Int_t ntimes)
        // select decay particles
              Int_t np=fDecayer->ImportParticles(particles);
 
+       if(fForceConv) np=ForceGammaConversion(particles,np);
+
+       // for(int iPart=0; iPart<np; iPart++){
+       //   TParticle *gamma = (TParticle *) particles->At(iPart);
+       //   printf("%i %i:", iPart, gamma->GetPdgCode());
+       //   printf("%i %i %i|",gamma->GetFirstMother(),gamma->GetFirstDaughter(),gamma->GetLastDaughter());
+       // }
+
              //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
              if (fGeometryAcceptance) 
                if (!CheckAcceptanceGeometry(np,particles)) continue;
index 591ed79..4b20a26 100644 (file)
@@ -46,11 +46,22 @@ public:
     virtual void SetWeighting(Weighting_t flag = kAnalog) {fAnalog = flag;}    
     virtual void SetDeltaPt(Float_t delta=0.01) {fDeltaPt = delta;}
     virtual void SetDecayer(AliDecayer* decayer) {fDecayer = decayer;}
+    virtual void SetForceGammaConversion(Bool_t force=kTRUE) {fForceConv = force;}
     virtual void Draw(const char * opt);
     TF1 *  GetPt() { return fPtPara;}
     TF1 *  GetY() {return fYPara;}
     Float_t GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax);
 
+    static double ScreenVar(double Z, double e0, double eps){ return 136/pow(Z,0.333333)*e0/eps/(1-eps); }
+    static double ScreenFunc1(double d);
+    static double ScreenFunc2(double d);
+    static double AuxScreenFunc1(double d, double Fz){ return 3*ScreenFunc1(d)-ScreenFunc2(d)-Fz; }
+    static double AuxScreenFunc2(double d, double Fz){ return 3*0.5*ScreenFunc1(d)-0.5*ScreenFunc2(d)-Fz; }
+    static TVector3 OrthogonalVector(TVector3 &inVec);
+    double EnergyFraction(double Z, double E);
+    double PolarAngle(double E);
+    Int_t ForceGammaConversion(TClonesArray *particles, Int_t nPart);
+  
 protected:
     Double_t (*fPtParaFunc)(const Double_t*, const Double_t*); //! Pointer to Pt parametrisation function
     Double_t (*fYParaFunc )(const Double_t*, const Double_t*); //! Pointer to Y parametrisation function
@@ -69,6 +80,7 @@ protected:
     Float_t     fDeltaPt;      // pT sampling in steps of fDeltaPt
     Bool_t      fSelectAll;    // Flag for transportation of Background while using SetForceDecay()
     AliDecayer  *fDecayer;     // ! Pointer to pythia object for decays
+    Bool_t      fForceConv;    //
 
 private:
     AliGenParam(const AliGenParam &Param);
index ae90bae..891f9f8 100644 (file)
@@ -429,6 +429,7 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  223,11,1); // omega     
        ForceParticleDecay(  331,11,1); // etaprime     
        ForceParticleDecay(  333,11,1); // phi     
+       ForceParticleDecay(  443,11,1); // jpsi     
        break;
     case kDiElectronEM:
        ForceParticleDecay(  111,11,2); // pi^0     
@@ -437,6 +438,7 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  223,11,2); // omega     
        ForceParticleDecay(  331,11,2); // etaprime     
        ForceParticleDecay(  333,11,2); // phi     
+       ForceParticleDecay(  443,11,2); // jpsi     
        break;
     case kGammaEM:
        ForceParticleDecay(  111,22,1); // pi^0     
@@ -445,6 +447,7 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  223,22,1); // omega     
        ForceParticleDecay(  331,22,1); // etaprime     
        ForceParticleDecay(  333,22,1); // phi     
+       ForceParticleDecay(  443,22,1); // jpsi     
        break;
     case kBeautyUpgrade:
         ForceBeautyUpgrade();
@@ -687,7 +690,10 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t m
     Int_t ifirst=fPythia->GetMDCY(kc,2);
     Int_t ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
     Double_t norm = 0.;
-    for (Int_t channel=ifirst; channel<=ilast;channel++) norm+=fPythia->GetBRAT(channel);
+    for (Int_t channel=ifirst; channel<=ilast;channel++){
+      norm+=fPythia->GetBRAT(channel);
+      printf("%f ",fPythia->GetBRAT(channel));
+    }printf("\n");
     if (norm < 1.-1.e-12 || norm > 1.+1.e-12) {
         char pName[16];
         fPythia->Pyname(particle,pName);