- finish implementation of direct gammas in the EM cocktail (prompt, thermal, real...
authormorsch <andreas.morsch@cern.ch>
Mon, 14 Apr 2014 13:42:40 +0000 (15:42 +0200)
committermorsch <andreas.morsch@cern.ch>
Mon, 14 Apr 2014 13:42:40 +0000 (15:42 +0200)
 Theodor Rascanu <trascanu@stud.uni-frankfurt.de>

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

index 2d73011..bb6f8d3 100644 (file)
@@ -61,7 +61,7 @@ AliGenEMCocktail::AliGenEMCocktail()
   fCentrality(0),
   fV2Systematic(0),
   fForceConv(kFALSE),
-  fHeaviestParticle(kGenJpsi)
+  fSelectedParticles(kGenHadrons)
 {
   // Constructor
 
@@ -97,7 +97,7 @@ void AliGenEMCocktail::CreateCocktail()
 
   // Create and add electron sources to the generator
   // pizero
-  if(fHeaviestParticle<kGenPizero)return;
+  if(fSelectedParticles&kGenPizero){
   AliGenParam *genpizero=0;
   Char_t namePizero[10];    
   snprintf(namePizero,10,"Pizero");    
@@ -106,13 +106,14 @@ void AliGenEMCocktail::CreateCocktail()
   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);
+    fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+    fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+  }
 
   // eta  
-  if(fHeaviestParticle<kGenEta)return;
+  if(fSelectedParticles&kGenEta){
   AliGenParam *geneta=0;
   Char_t nameEta[10];    
   snprintf(nameEta,10,"Eta");    
@@ -121,13 +122,14 @@ void AliGenEMCocktail::CreateCocktail()
   AddSource2Generator(nameEta,geneta);
   TF1 *fPtEta = geneta->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
-  fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
+    fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+    fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+  }
 
   // rho  
-  if(fHeaviestParticle<kGenRho)return;
+  if(fSelectedParticles&kGenRho){
   AliGenParam *genrho=0;
   Char_t nameRho[10];    
   snprintf(nameRho,10,"Rho");    
@@ -136,13 +138,14 @@ void AliGenEMCocktail::CreateCocktail()
   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);
+    fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+    fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+  }
   
   // omega
-  if(fHeaviestParticle<kGenOmega)return;
+  if(fSelectedParticles&kGenOmega){
   AliGenParam *genomega=0;
   Char_t nameOmega[10];    
   snprintf(nameOmega,10,"Omega");    
@@ -151,13 +154,14 @@ void AliGenEMCocktail::CreateCocktail()
   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);
+    fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+    fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+  }
 
   // etaprime
-  if(fHeaviestParticle<kGenEtaprime)return;
+  if(fSelectedParticles&kGenEtaprime){
   AliGenParam *genetaprime=0;
   Char_t nameEtaprime[10];    
   snprintf(nameEtaprime,10,"Etaprime");    
@@ -166,13 +170,14 @@ void AliGenEMCocktail::CreateCocktail()
   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);
+    fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+    fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+  }
 
   // phi  
-  if(fHeaviestParticle<kGenPhi)return;
+  if(fSelectedParticles&kGenPhi){
   AliGenParam *genphi=0;
   Char_t namePhi[10];    
   snprintf(namePhi,10,"Phi");    
@@ -181,13 +186,14 @@ void AliGenEMCocktail::CreateCocktail()
   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);
+    fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+    fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+  }
 
   // jpsi  
-  if(fHeaviestParticle<kGenJpsi)return;
+  if(fSelectedParticles&kGenJpsi){
   AliGenParam *genjpsi=0;
   Char_t nameJpsi[10];    
   snprintf(nameJpsi,10,"Jpsi");    
@@ -196,76 +202,44 @@ void AliGenEMCocktail::CreateCocktail()
   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);
+    fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
 #else
-  fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
-#endif
-
-  // prompt gamma
-  if(fDecayMode==kGammaEM){
-    if(fHeaviestParticle<kGenPromptRealGamma)return;
-    TDatabasePDG::Instance()->AddParticle("PromptRealGamma","PromptRealGamma",0,true,0,0,"GaugeBoson",221000);
-    //gMC->DefineParticle(221000, "PromptGamma", kPTGamma, 0, 0, 0,"Gamma", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kFALSE);
-    AliGenParam *genPromptRealG=0;
-    Char_t namePromptRealG[10];    
-    snprintf(namePromptRealG,10,"PromptRealGamma");    
-    genPromptRealG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kPromptRealGamma, "DUMMY");
-    genPromptRealG->SetYRange(fYMin, fYMax);
-    AddSource2Generator(namePromptRealG,genPromptRealG);
-    TF1 *fPtPromptRealG = genPromptRealG->GetPt();
-#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
-    fYieldArray[kGenPromptRealGamma] = fPtPromptRealG->Integral(fPtMin,fPtMax,1.e-6);
-#else
-    fYieldArray[kGenPromptRealGamma] = fPtPromptRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
-#endif
-
-    if(fHeaviestParticle<kGenThermRealGamma)return;
-    TDatabasePDG::Instance()->AddParticle("ThermRealGamma","ThermRealGamma",0,true,0,0,"GaugeBoson",222000);
-    //gMC->DefineParticle(221000, "ThermGamma", kPTGamma, 0, 0, 0,"Gamma", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kFALSE);
-    AliGenParam *genThermRealG=0;
-    Char_t nameThermRealG[10];    
-    snprintf(nameThermRealG,10,"ThermRealGamma");    
-    genThermRealG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kThermRealGamma, "DUMMY");
-    genThermRealG->SetYRange(fYMin, fYMax);
-    AddSource2Generator(nameThermRealG,genThermRealG);
-    TF1 *fPtThermRealG = genThermRealG->GetPt();
-#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
-    fYieldArray[kGenThermRealGamma] = fPtThermRealG->Integral(fPtMin,fPtMax,1.e-6);
-#else
-    fYieldArray[kGenThermRealGamma] = fPtThermRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+    fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+  }
 
-    if(fHeaviestParticle<kGenPromptVirtGamma)return;
-    TDatabasePDG::Instance()->AddParticle("PromptVirtGamma","PromptVirtGamma",0,true,0,0,"GaugeBoson",223000);
-    AliGenParam *genPromptVirtG=0;
-    Char_t namePromptVirtG[10];    
-    snprintf(namePromptVirtG,10,"PromptVirtGamma");    
-    genPromptVirtG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kPromptVirtGamma, "DUMMY");
-    genPromptVirtG->SetYRange(fYMin, fYMax);
-    AddSource2Generator(namePromptVirtG,genPromptVirtG);
-    TF1 *fPtPromptVirtG = genPromptVirtG->GetPt();
+  // direct gamma
+  if(fDecayMode!=kGammaEM) return;
+  if(fSelectedParticles&kGenDirectRealGamma){
+    AliGenParam *genDirectRealG=0;
+    Char_t nameDirectRealG[10];    
+    snprintf(nameDirectRealG,10,"DirectRealGamma");    
+    genDirectRealG = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDirectRealGamma, "DUMMY");
+    genDirectRealG->SetYRange(fYMin, fYMax);
+    AddSource2Generator(nameDirectRealG,genDirectRealG);
+    TF1 *fPtDirectRealG = genDirectRealG->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
-    fYieldArray[kGenPromptVirtGamma] = fPtPromptVirtG->Integral(fPtMin,fPtMax,1.e-6);
+    fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
 #else
-    fYieldArray[kGenPromptVirtGamma] = fPtPromptVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+    fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+  }
 
-    if(fHeaviestParticle<kGenThermVirtGamma)return;
-    TDatabasePDG::Instance()->AddParticle("ThermVirtGamma","ThermVirtGamma",0,true,0,0,"GaugeBoson",224000);
-    AliGenParam *genThermVirtG=0;
-    Char_t nameThermVirtG[10];    
-    snprintf(nameThermVirtG,10,"ThermVirtGamma");    
-    genThermVirtG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kThermVirtGamma, "DUMMY");
-    genThermVirtG->SetYRange(fYMin, fYMax);
-    AddSource2Generator(nameThermVirtG,genThermVirtG);
-    TF1 *fPtThermVirtG = genThermVirtG->GetPt();
+  if(fSelectedParticles&kGenDirectVirtGamma){
+    TDatabasePDG::Instance()->AddParticle("DirectVirtGamma","DirectVirtGamma",0,true,0,0,"GaugeBoson",220000);
+    AliGenParam *genDirectVirtG=0;
+    Char_t nameDirectVirtG[10];    
+    snprintf(nameDirectVirtG,10,"DirectVirtGamma");    
+    genDirectVirtG = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kDirectVirtGamma, "DUMMY");
+    genDirectVirtG->SetYRange(fYMin/0.775, fYMax/0.775);
+    AddSource2Generator(nameDirectVirtG,genDirectVirtG);
+    TF1 *fPtDirectVirtG = genDirectVirtG->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
-    fYieldArray[kGenThermVirtGamma] = fPtThermVirtG->Integral(fPtMin,fPtMax,1.e-6);
+    fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,1.e-6);
 #else
-    fYieldArray[kGenThermVirtGamma] = fPtThermVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+    fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
   }
-  
 }
 
 //-------------------------------------------------------------------
@@ -362,44 +336,37 @@ void AliGenEMCocktail::Generate()
 
     switch (pdgMother){
     case 111:
-      dNdy = fYieldArray[kGenPizero];
+      dNdy = fYieldArray[kPizero];
       break;
     case 221:
-      dNdy = fYieldArray[kGenEta];
+      dNdy = fYieldArray[kEta];
       break;
     case 113:
-      dNdy = fYieldArray[kGenRho];
+      dNdy = fYieldArray[kRho];
       break;
     case 223:
-      dNdy = fYieldArray[kGenOmega];
+      dNdy = fYieldArray[kOmega];
       break;
     case 331:
-      dNdy = fYieldArray[kGenEtaprime];
+      dNdy = fYieldArray[kEtaprime];
       break;
     case 333:
-      dNdy = fYieldArray[kGenPhi];
+      dNdy = fYieldArray[kPhi];
       break;
     case 443:
-      dNdy = fYieldArray[kGenJpsi];
+      dNdy = fYieldArray[kJpsi];
+      break;
+    case 22:
+      dNdy = fYieldArray[kDirectRealGamma];
       break;
+    case 220000:
+      dNdy = fYieldArray[kDirectVirtGamma];
+      break;
+      
     default:
       dNdy = 0.;
     }
 
-    switch (pdgMother){
-    case 221000:
-      dNdy = fYieldArray[kGenPromptRealGamma];
-      break;
-    case 223000:
-      dNdy = fYieldArray[kGenPromptVirtGamma];
-      break;
-    case 222000:
-      dNdy = fYieldArray[kGenThermRealGamma];
-      break;      
-    case 224000:
-      dNdy = fYieldArray[kGenThermVirtGamma];
-      break;   
-    }
     weight = dNdy*part->GetWeight();
     part->SetWeight(weight);
   }    
index 96ed249..89b7a50 100644 (file)
@@ -23,7 +23,8 @@ class AliGenEMCocktail : public AliGenCocktail
 public:
 
     AliGenEMCocktail();
-  enum GeneratorCode { kGenPizero=0, kGenEta, kGenRho, kGenOmega, kGenEtaprime, kGenPhi, kGenJpsi, kGenPromptRealGamma, kGenThermRealGamma, kGenPromptVirtGamma, kGenThermVirtGamma, kGENs };
+  enum GeneratorCode { kPizero=0, kEta, kRho, kOmega, kEtaprime, kPhi, kJpsi, kDirectRealGamma, kDirectVirtGamma, kGENs };
+  enum SelectPartice { kGenPizero=0x001, kGenEta=0x002, kGenRho=0x004, kGenOmega=0x008, kGenEtaprime=0x010, kGenPhi=0x020, kGenJpsi=0x040, kGenDirectRealGamma=0x100, kGenDirectVirtGamma=0x200, kGenHadrons=0x7f, kGenGammas=0x300 };
 
     virtual ~AliGenEMCocktail();    
     virtual void Init();
@@ -39,7 +40,7 @@ public:
   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; }
+  void    SelectMotherParticles(Int_t part){ fSelectedParticles=part; }
     
 private:
     AliGenEMCocktail(const AliGenEMCocktail &cocktail); 
@@ -59,7 +60,7 @@ private:
   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
+  Int_t fSelectedParticles; //which particles to simulate
 
     ClassDef(AliGenEMCocktail,1)  //  cocktail for EM physics
 };
index 9c80a67..ff46597 100644 (file)
@@ -50,15 +50,15 @@ Double_t AliGenEMlib::CrossOverRc(const double a, const double b, const double x
 }
 
 const Double_t AliGenEMlib::fgkV2param[16][15] = {
-  // charged pion                                                                                                                        cent, based on
+  // charged pion                                                                                                                        cent, based on: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/FlowPAGQM2012talkIdentified
   {  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.551541e-02, 1.438274e+00, 4.626379e-02, 2.512477e+00, 1.371824e+00, 2.964543e-02, 4.630670e+00, 4.228889e+00, 6.037970e-02, 1.425269e-03, 1.144124e+00, 0, 1, 9.154016e-04, 1.288285e+00 }  // 0-5,  Francesco
-  ,{ 1.171360e-01, 1.333046e+00, 4.536752e-02, 3.046448e+00, 3.903714e+00, 4.407124e-02, 9.122534e-01, 4.834519e+00, 1.186237e-01, 2.179274e-03, 8.968478e-01, 0, 1, 1.501201e-03, 9.902785e-01 }  // 5-10, Francesco
-  ,{ 1.748423e-01, 1.285211e+00, 4.219624e-02, 4.019148e+00, 4.255047e+00, 7.956751e-03, 1.184731e-01,-9.211391e+00, 5.768716e-01, 3.127110e-03, 6.808650e-01, 0, 1, 2.786807e-03, 6.159338e-01 }  // 10-20,Francesco
-  ,{ 2.152937e-01, 1.405391e+00, 5.037925e-02, 3.214458e+00, 3.991894e+00, 3.655882e-02, 1.968766e-01,-1.637650e+01, 7.023397e+00, 4.573453e-03, 6.031381e-01, 0, 1, 3.564348e-03, 5.748053e-01 }  // 20-30,Francesco
-  ,{ 2.409800e-01, 1.476557e+00, 5.759362e-02, 3.339713e+00, 3.642386e+00,-1.544366e-02, 1.098611e-01,-1.373154e+01, 1.471955e+00, 5.200180e-03, 6.315474e-01, 0, 1, 3.776112e-03, 6.298605e-01 }  // 30-40,Francesco
-  ,{ 2.495087e-01, 1.543711e+00, 6.217817e-02, 3.517101e+00, 4.558221e+00, 6.021316e-02, 1.486822e-01,-5.769155e+00, 5.576843e-01, 5.348029e-03, 7.255976e-01, 0, 1, 3.531350e-03, 7.661694e-01 }  // 40-50,Francesco
-  ,{ 2.166449e-01, 1.931014e+00, 8.195656e-02, 2.226742e+00, 3.106472e+00, 1.058786e-01, 8.558786e-01, 4.006680e+00, 2.476313e-01, 5.137623e-03, 9.104401e-01, 0, 1, 2.477450e-03, 1.109649e+00 }  // 50-60,Francesco
+  ,{ 6.551541e-02, 1.438274e+00, 4.626379e-02, 2.512477e+00, 1.371824e+00, 2.964543e-02, 4.630670e+00, 4.228889e+00, 6.037970e-02, 1.425269e-03, 1.144124e+00, 0, 1, 9.154016e-04, 1.288285e+00 }  // 0-5
+  ,{ 1.171360e-01, 1.333046e+00, 4.536752e-02, 3.046448e+00, 3.903714e+00, 4.407124e-02, 9.122534e-01, 4.834519e+00, 1.186237e-01, 2.179274e-03, 8.968478e-01, 0, 1, 1.501201e-03, 9.902785e-01 }  // 5-10
+  ,{ 1.748423e-01, 1.285211e+00, 4.219624e-02, 4.019148e+00, 4.255047e+00, 7.956751e-03, 1.184731e-01,-9.211391e+00, 5.768716e-01, 3.127110e-03, 6.808650e-01, 0, 1, 2.786807e-03, 6.159338e-01 }  // 10-20
+  ,{ 2.152937e-01, 1.405391e+00, 5.037925e-02, 3.214458e+00, 3.991894e+00, 3.655882e-02, 1.968766e-01,-1.637650e+01, 7.023397e+00, 4.573453e-03, 6.031381e-01, 0, 1, 3.564348e-03, 5.748053e-01 }  // 20-30
+  ,{ 2.409800e-01, 1.476557e+00, 5.759362e-02, 3.339713e+00, 3.642386e+00,-1.544366e-02, 1.098611e-01,-1.373154e+01, 1.471955e+00, 5.200180e-03, 6.315474e-01, 0, 1, 3.776112e-03, 6.298605e-01 }  // 30-40
+  ,{ 2.495087e-01, 1.543711e+00, 6.217817e-02, 3.517101e+00, 4.558221e+00, 6.021316e-02, 1.486822e-01,-5.769155e+00, 5.576843e-01, 5.348029e-03, 7.255976e-01, 0, 1, 3.531350e-03, 7.661694e-01 }  // 40-50
+  ,{ 2.166449e-01, 1.931014e+00, 8.195656e-02, 2.226742e+00, 3.106472e+00, 1.058786e-01, 8.558786e-01, 4.006680e+00, 2.476313e-01, 5.137623e-03, 9.104401e-01, 0, 1, 2.477450e-03, 1.109649e+00 }  // 50-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 }   // 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
@@ -92,17 +92,17 @@ const Double_t AliGenEMlib::fgkThermPtParam[16][2] = {
   {  0.0000000000, 0.0000000000 } // pp no V2
   ,{ 0.0000000000, 0.0000000000 } // 0-5
   ,{ 0.0000000000, 0.0000000000 } // 5-10
-  ,{ 2.581823e+01, 3.187900e+00 } // 10-20 //from: https://aliceinfo.cern.ch/Notes/node/249
+  ,{ 3.447105e+01, 3.416818e+00 } // 10-20 //based on: https://aliceinfo.cern.ch/Notes/node/249
   ,{ 0.0000000000, 0.0000000000 } // 20-30
   ,{ 0.0000000000, 0.0000000000 } // 30-40
   ,{ 0.0000000000, 0.0000000000 } // 40-50
   ,{ 0.0000000000, 0.0000000000 } // 50-60
-  ,{ 7.177551e+02, 4.946179e+00 } // 0-10  //from: https://aliceinfo.cern.ch/Notes/node/249
-  ,{ 2.328661e+00, 2.635257e+00 } // 20-40 //from: https://twiki.cern.ch/twiki/pub/ALICE/ALICEDirectPhotonSpectrumPaper/directPbPb.pdf
+  ,{ 3.888847e+02, 4.502683e+00 } // 0-10  //based on: https://aliceinfo.cern.ch/Notes/node/249
+  ,{ 1.766210e+00, 2.473812e+00 } // 20-40 //based on: https://twiki.cern.ch/twiki/pub/ALICE/ALICEDirectPhotonSpectrumPaper/directPbPb.pdf
   ,{ 0.0000000000, 0.0000000000 } // 40-60
   ,{ 0.0000000000, 0.0000000000 } // 60-80
-  ,{ 1.919280e+01, 2.946472e+00 } // 0-20  //from: https://twiki.cern.ch/twiki/pub/ALICE/ALICEDirectPhotonSpectrumPaper/directPbPb.pdf
-  ,{ 0.0000000000, 0.0000000000 } // 0-40 
+  ,{ 1.576151e+01, 2.841202e+00 } // 0-20  //based on: https://twiki.cern.ch/twiki/pub/ALICE/ALICEDirectPhotonSpectrumPaper/directPbPb.pdf
+  ,{ 4.263499e+01, 3.249843e+00 } // 0-40  //based on: https://aliceinfo.cern.ch/Figure/node/2866
   ,{ 0.0000000000, 0.0000000000 } // 20-80
   ,{ 0.0000000000, 0.0000000000 } // 40-80
 };
@@ -119,11 +119,10 @@ const Double_t AliGenEMlib::fgkMtFactor[2][8] = {
   //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
   //https://aliceinfo.cern.ch/Notes/node/87
   //best guess:
-  {1., 0.48, 1.0, 0.9, 0.4, 0.25, 0.004, 0.}, //pp
-  {1., 0.48, 1.0, 0.9, 0.4, 0.25, 0.0195, 0.}  //PbPb
+  {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
 };
 
 //==========================================================================
@@ -211,120 +210,108 @@ Double_t AliGenEMlib::PtExponential(const Double_t *px, const Double_t *c){
 
 // Hagedorn with additional Powerlaw
 Double_t AliGenEMlib::PtModifiedHagedornPowerlaw(const Double_t *px, const Double_t *c){
-  double pt=px[0]+0.0001;
-  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,-c[9]);
+  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);
 }
 
-Double_t AliGenEMlib::IntegratedKrollWada(Double_t mh){
-  if(mh<0.003) return 0;
-  const double me=0.000511;
-  return 2*log(mh/me/exp(7.0/4.0))/411.0/TMath::Pi();
+// 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;
+}
+
+// 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();
 }
 
 //--------------------------------------------------------------------------
 //
-//                             PromptRealGamma
+//                             DirectRealGamma
 //
 //--------------------------------------------------------------------------
-Int_t AliGenEMlib::IpPromptRealGamma(TRandom *)
-{
-  return 221000;
-}
-
 Double_t AliGenEMlib::PtPromptRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  //if(*px<0.001) return 0;
-  const static Double_t promptGammaPtParam[10] = { 7.019259e-02, 6.771695e-01, 8.249346e-01, 5.720419e+00, 1.848869e+01, 2.629075e+01, 1.061126e+01, 3.699205e+01, 5.253572e-02, 5.167275e+00 };
-  //{ 5.449971e-02, 3.843241e-01, 9.469766e-01, 4.993039e+00, 5.342451e+00, 4.457944e+00, 5.555146e+00, 4.458580e+00, 6.035177e-02, 5.102109e+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);
 }
 
-Double_t AliGenEMlib::YPromptRealGamma( const Double_t *px, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::PtThermalRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return YFlat(*px);
-}
-
-Double_t AliGenEMlib::V2PromptRealGamma( const Double_t */*px*/, const Double_t */*dummy*/ )
-{
-  return 0.0;
+  return PtExponential(px,fgkThermPtParam[fgSelectedCentrality]);
 }
 
-//--------------------------------------------------------------------------
-//
-//                             PromptVirtGamma
-//
-//--------------------------------------------------------------------------
-Int_t AliGenEMlib::IpPromptVirtGamma(TRandom *)
+Double_t AliGenEMlib::PtDirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return 223000;
+  return PtPromptRealGamma(px,px)+PtThermalRealGamma(px,px);
 }
 
-Double_t AliGenEMlib::PtPromptVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+Int_t AliGenEMlib::IpDirectRealGamma(TRandom *)
 {
-  return IntegratedKrollWada(*px)*PtPromptRealGamma(px,px);
+  return 22;
 }
 
-Double_t AliGenEMlib::YPromptVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::YDirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
   return YFlat(*px);
 }
 
-Double_t AliGenEMlib::V2PromptVirtGamma( const Double_t */*px*/, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::V2DirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return 0.0;
+  const static Double_t v2Param[3][15] = {
+    { 1.211795e-01, 9.813671e-01, 0.000000e+00, 3.056960e+00, 2.380183e+00, -7.833212e-02, 5.000000e-01, 3.056960e+00, 1.195000e-01, 1.183293e-02, 1.252249e+00, 0, 1, 4.876263e-03, 1.518526e+00 }  // 00-20, based on: https://aliceinfo.cern.ch/Notes/node/249
+    ,{ 1.619000e-01, 2.185695e+00, 0.000000e+00, 1.637681e+00, 1.000000e+00, -1.226399e-06, 3.092027e+00, 3.064692e+00, 1.619000e-01, 2.264320e-02, 1.028641e+00, 0, 1, 8.172203e-03, 1.271637e+00 } // 20-40
+    ,{ 1.335000e-01, 1.331963e+00, 0.000000e+00, 2.252315e+00, 1.198383e+00, -5.861987e-02, 7.132859e-01, 2.252315e+00, 2.934249e-01, 1.571589e-02, 1.001131e+00, 0, 1, 5.179715e-03, 1.329344e+00 } // 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;
+  }
+  return 0;
 }
 
+
 //--------------------------------------------------------------------------
 //
-//                             ThermRealGamma
+//                             DirectVirtGamma
 //
 //--------------------------------------------------------------------------
-Int_t AliGenEMlib::IpThermRealGamma(TRandom *)
-{
-  return 222000;
-}
-
-Double_t AliGenEMlib::PtThermRealGamma( const Double_t *px, const Double_t */*dummy*/ )
-{
-  return PtExponential(px,fgkThermPtParam[fgSelectedCentrality]);
-}
-
-Double_t AliGenEMlib::YThermRealGamma( const Double_t *px, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::PtPromptVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return YFlat(*px);
+  return IntegratedKrollWada(px,px)*PtPromptRealGamma(px,px);
 }
 
-Double_t AliGenEMlib::V2ThermRealGamma( const Double_t *px, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::PtThermalVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return KEtScal(*px,8);
+  return IntegratedKrollWada(px,px)*PtThermalRealGamma(px,px);
 }
 
-//--------------------------------------------------------------------------
-//
-//                             ThermVirtGamma
-//
-//--------------------------------------------------------------------------
-Int_t AliGenEMlib::IpThermVirtGamma(TRandom *)
+Double_t AliGenEMlib::PtDirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return 224000;
+  return IntegratedKrollWada(px,px)*(PtPromptRealGamma(px,px)+PtThermalRealGamma(px,px));
 }
 
-Double_t AliGenEMlib::PtThermVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+Int_t AliGenEMlib::IpDirectVirtGamma(TRandom *)
 {
-  return IntegratedKrollWada(*px)*PtThermRealGamma(px,px);
+  return 220000;
 }
 
-Double_t AliGenEMlib::YThermVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::YDirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
   return YFlat(*px);
 }
 
-Double_t AliGenEMlib::V2ThermVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+Double_t AliGenEMlib::V2DirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 {
-  return KEtScal(*px,8);
+  return V2DirectRealGamma(px,px);
 }
 
 //--------------------------------------------------------------------------
@@ -340,9 +327,12 @@ Int_t AliGenEMlib::IpPizero(TRandom *)
 
 Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
 {
-  // double pigammacorr=2.385389e-01*log(*px+0.001)+1.557687e+00;
-  // pigammacorr*=9.513666e-03*log(*px+0.001)+9.509347e-01;
-  // return pigammacorr*PtPromptRealGamma(px,px);  //misuse pion for direct gammas
+  // 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 
 
@@ -358,6 +348,9 @@ Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
   Double_t kb=0.;
   Double_t kd=0.;
 
+  double n1,n2,n3;
+  int oldCent;
+
   switch(fgSelectedPtParam|fgSelectedCentrality) {
     // fit to pi charged v1
     // charged pion from ToF, unidentified hadrons scaled with pion from TPC
@@ -407,7 +400,26 @@ Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
     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 kPichargedPbPb|k0020:
+    oldCent=fgSelectedCentrality;
+    fgSelectedCentrality=k0010;
+    n1=PtPizero(px,px);
+    fgSelectedCentrality=k1020;
+    n2=PtPizero(px,px);
+    fgSelectedCentrality=oldCent;
+    return (n1+n2)/2;
+    break;
+  case kPichargedPbPb|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
@@ -502,7 +514,8 @@ Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
 
 Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 {
-  double n1,n2,v1,v2;
+  double n1,n2,n3,n4,n5;
+  double v1,v2,v3,v4,v5;
   switch(fgSelectedCentrality) {
   case k0010:
     n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
@@ -511,6 +524,15 @@ Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
     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]);
@@ -518,6 +540,19 @@ Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
     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]);
@@ -717,7 +752,18 @@ Int_t AliGenEMlib::IpJpsi(TRandom *)
 Double_t AliGenEMlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/ )
 {
   // Jpsi pT
-  return MtScal(*px,6);
+  // 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;
 }
 
 Double_t AliGenEMlib::YJpsi( const Double_t *py, const Double_t */*dummy*/ )
@@ -786,7 +832,7 @@ Double_t AliGenEMlib::KEtScal(const Double_t pt, Int_t np)
   // 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);
+  // return std::max(val+sys,0.0);
   return V2Pizero(&scaledPt, (Double_t*) 0);
 }
 
@@ -797,7 +843,7 @@ Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
   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);
+  return std::max(val+sys,0.0);
 }
 
 Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
@@ -845,17 +891,11 @@ GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
 
    switch (param) 
     {
-    case kPromptRealGamma:
-      func=PtPromptRealGamma;
-      break;
-    case kPromptVirtGamma:
-      func=PtPromptVirtGamma;
+     case kDirectRealGamma:
+       func=PtDirectRealGamma;
       break;
-    case kThermRealGamma:
-      func=PtThermRealGamma;
-      break;
-    case kThermVirtGamma:
-      func=PtThermVirtGamma;
+     case kDirectVirtGamma:
+       func=PtDirectVirtGamma;
       break;
     case kPizero:
       func=PtPizero;
@@ -894,17 +934,11 @@ GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
 
    switch (param) 
     {
-    case kPromptRealGamma:
-      func=YPromptRealGamma;
-      break;
-    case kPromptVirtGamma:
-      func=YPromptVirtGamma;
+    case kDirectRealGamma:
+      func=YDirectRealGamma;
       break;
-    case kThermRealGamma:
-      func=YThermRealGamma;
-      break;
-    case kThermVirtGamma:
-      func=YThermVirtGamma;
+    case kDirectVirtGamma:
+      func=YDirectVirtGamma;
       break;
     case kPizero:
          func=YPizero;
@@ -943,17 +977,11 @@ GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
 
    switch (param) 
     {
-    case kPromptRealGamma:
-      func=IpPromptRealGamma;
-      break;
-    case kPromptVirtGamma:
-      func=IpPromptVirtGamma;
+    case kDirectRealGamma:
+      func=IpDirectRealGamma;
       break;
-    case kThermRealGamma:
-      func=IpThermRealGamma;
-      break;
-    case kThermVirtGamma:
-      func=IpThermVirtGamma;
+    case kDirectVirtGamma:
+      func=IpDirectVirtGamma;
       break;
     case kPizero:
          func=IpPizero;
@@ -992,17 +1020,11 @@ GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
 
   switch (param) 
     {
-    case kPromptRealGamma:
-      func=V2PromptRealGamma;
-      break;
-    case kPromptVirtGamma:
-      func=V2PromptVirtGamma;
-      break;
-    case kThermRealGamma:
-      func=V2ThermRealGamma;
+    case kDirectRealGamma:
+      func=V2DirectRealGamma;
       break;
-    case kThermVirtGamma:
-      func=V2ThermVirtGamma;
+    case kDirectVirtGamma:
+      func=V2DirectVirtGamma;
       break;
     case kPizero:
       func=V2Pizero;
index 6372507..95af158 100644 (file)
@@ -21,7 +21,7 @@ class TRandom;
 class AliGenEMlib :public AliGenLib {
 public:
     
-  enum Particle_t{kPromptRealGamma, kPromptVirtGamma, kThermRealGamma, kThermVirtGamma, kPizero, kEta, kRho, kOmega, kEtaprime, kPhi, kJpsi};
+  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=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_t{kLoV2Sys=-1, kNoV2Sys=0, kUpV2Sys=+1};
@@ -81,29 +81,24 @@ public:
 
   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 IntegratedKrollWada(Double_t mh);
+  static Double_t PtDoublePowerlaw(const Double_t *pt, const Double_t *param);
+  static Double_t IntegratedKrollWada(const Double_t *mh, const Double_t *);
 
-  // prompt gamma
-  static Int_t    IpPromptRealGamma(TRandom *ran);
+  // direct gamma
   static Double_t PtPromptRealGamma(const Double_t *px, const Double_t *dummy);
-  static Double_t YPromptRealGamma(const Double_t *py, const Double_t *dummy);
-  static Double_t V2PromptRealGamma(const Double_t *px, const Double_t *dummy);
-
-  static Int_t    IpPromptVirtGamma(TRandom *ran);
   static Double_t PtPromptVirtGamma(const Double_t *px, const Double_t *dummy);
-  static Double_t YPromptVirtGamma(const Double_t *py, const Double_t *dummy);
-  static Double_t V2PromptVirtGamma(const Double_t *px, const Double_t *dummy);
-
-  // thermal gamma
-  static Int_t    IpThermRealGamma(TRandom *ran);
-  static Double_t PtThermRealGamma(const Double_t *px, const Double_t *dummy);
-  static Double_t YThermRealGamma(const Double_t *py, const Double_t *dummy);
-  static Double_t V2ThermRealGamma(const Double_t *px, const Double_t *dummy);
-
-  static Int_t    IpThermVirtGamma(TRandom *ran);
-  static Double_t PtThermVirtGamma(const Double_t *px, const Double_t *dummy);
-  static Double_t YThermVirtGamma(const Double_t *py, const Double_t *dummy);
-  static Double_t V2ThermVirtGamma(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);
index 3e92b14..a2dafd2 100644 (file)
@@ -225,11 +225,6 @@ void AliGenParam::RotateVector(Double_t *pin, Double_t *pout, Double_t costheta,
   return;
 }
 
-Double_t AliGenParam::IntegratedKrollWada(Double_t mh){
-  if(mh<0.003) return 0;
-  return 2*log(mh/0.000511/exp(7.0/4.0))/411.0/TMath::Pi();
-}
-
 double AliGenParam::ScreenFunction1(double screenVariable){
   if(screenVariable>1)
     return 42.24 - 8.368 * log(screenVariable + 0.952);
@@ -321,8 +316,8 @@ Double_t AliGenParam::RandomMass(Double_t mh){
   while(true){
     double y=fRandom->Rndm();
     double mee=2*0.000511*TMath::Power(2*0.000511/mh,-y); //inverse of the enveloping cumulative distribution
-    double apxkw=2.0/3.0/137.036/TMath::Pi()/mee;
-    double val=fRandom->Uniform(0,apxkw);  //enveloping probability density0
+    double apxkw=2.0/3.0/137.036/TMath::Pi()/mee; //enveloping probability density
+    double val=fRandom->Uniform(0,apxkw);
     double kw=apxkw*sqrt(1-4*0.000511*0.000511/mee/mee)*(1+2*0.000511*0.000511/mee/mee)*1*1*TMath::Power(1-mee*mee/mh/mh,3);
     if(val<kw)
       return mee;
@@ -334,8 +329,8 @@ Int_t AliGenParam::VirtualGammaPairProduction(TClonesArray *particles, Int_t nPa
   Int_t nPartNew=nPart;
   for(int iPart=0; iPart<nPart; iPart++){      
     TParticle *gamma = (TParticle *) particles->At(iPart);
-    if(gamma->GetPdgCode()!=223000 || gamma->GetPdgCode()!=224000) continue;
-    //if(gamma->Energy()<0.001022) continue; //can never be
+    if(gamma->GetPdgCode()!=220000) continue;
+    if(gamma->Pt()<0.002941) continue;  //approximation of kw in AliGenEMlib is 0 below 0.002941
     double mass=RandomMass(gamma->Pt());
 
     // lepton pair kinematics in virtual photon rest frame 
@@ -378,6 +373,7 @@ Int_t AliGenParam::VirtualGammaPairProduction(TClonesArray *particles, Int_t nPa
     new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, e2V4, vtx);
     nPartNew++;
   }
+  return nPartNew;
 }
 
 Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
@@ -390,7 +386,6 @@ Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
     TParticle *gamma = (TParticle *) particles->At(iPart);
     if(gamma->GetPdgCode()!=22) continue;
     if(gamma->Energy()<0.001022) continue;
-
     TVector3 gammaV3(gamma->Px(),gamma->Py(),gamma->Pz());
     double frac=RandomEnergyFraction(1,gamma->Energy());
     double Ee1=frac*gamma->Energy();
@@ -424,8 +419,7 @@ Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
     new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e2V3,Ee2), vtx);
     nPartNew++;
   }
-  // particles->Compress();
-  return particles->GetEntriesFast();
+  return nPartNew;
 }
 
 //____________________________________________________________
@@ -579,7 +573,7 @@ void AliGenParam::GenerateN(Int_t ntimes)
       Int_t iTemp = iPart;
 
       // custom pdg codes for to destinguish direct photons
-      if((iPart>=221000) && (iPart<=229000)) iPart=22;
+      if(iPart==220000) iPart=22;
 
          fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;          
          TParticlePDG *particle = pDataBase->GetParticle(iPart);
@@ -664,13 +658,12 @@ void AliGenParam::GenerateN(Int_t ntimes)
              Int_t np=fDecayer->ImportParticles(particles);
 
        iPart=iTemp;
-       if(fForceConv) np=ForceGammaConversion(particles,np);
-       if(iPart==223000 || iPart==224000){
-         // wgtp*=IntegratedKrollWada(pt);
-         // wgtch*=IntegratedKrollWada(pt);
-         // np=VirtualGammaPairProduction(particles,np)
+       if(iPart==220000){
+         TParticle *gamma = (TParticle *)particles->At(0);
+         gamma->SetPdgCode(iPart);
+         np=VirtualGammaPairProduction(particles,np);
        }
-
+       if(fForceConv) np=ForceGammaConversion(particles,np);
 
              //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
              if (fGeometryAcceptance) 
index 522096a..b00789c 100644 (file)
@@ -55,7 +55,6 @@ public:
     static TVector3 OrthogonalVector(TVector3 &inVec);
     static void RotateVector(Double_t *pin, Double_t *pout, Double_t costheta, Double_t sintheta,
                           Double_t cosphi, Double_t sinphi);
-    static double IntegratedKrollWada(Double_t mh);
     static double ScreenFunction1(double d);
     static double ScreenFunction2(double d);
     double RandomEnergyFraction(double Z, double E);