- correct minor issues with forcred gamma decays, and speed up
authormorsch <andreas.morsch@cern.ch>
Mon, 31 Mar 2014 16:07:47 +0000 (18:07 +0200)
committermorsch <andreas.morsch@cern.ch>
Mon, 31 Mar 2014 16:07:47 +0000 (18:07 +0200)
- correct mt scaling factors
- first implementation for direct gammas
Theodor Rascanu

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

index 5bb28da..2d73011 100644 (file)
@@ -31,6 +31,7 @@
 #include <TF1.h>
 #include <TVirtualMC.h>
 #include <TPDGCode.h>
+#include <TDatabasePDG.h>
 #include "AliGenCocktailEventHeader.h"
 
 #include "AliGenCocktailEntry.h"
@@ -60,7 +61,7 @@ AliGenEMCocktail::AliGenEMCocktail()
   fCentrality(0),
   fV2Systematic(0),
   fForceConv(kFALSE),
-  fHeaviestParticle(kGENs)
+  fHeaviestParticle(kGenJpsi)
 {
   // Constructor
 
@@ -95,12 +96,13 @@ void AliGenEMCocktail::CreateCocktail()
   AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
 
   // Create and add electron sources to the generator
-
   // pizero
+  if(fHeaviestParticle<kGenPizero)return;
   AliGenParam *genpizero=0;
   Char_t namePizero[10];    
   snprintf(namePizero,10,"Pizero");    
-  genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
+  genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
+  genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
   AddSource2Generator(namePizero,genpizero);
   TF1 *fPtPizero = genpizero->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
@@ -114,7 +116,8 @@ void AliGenEMCocktail::CreateCocktail()
   AliGenParam *geneta=0;
   Char_t nameEta[10];    
   snprintf(nameEta,10,"Eta");    
-  geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
+  geneta = new AliGenParam(fNPart/0.825, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
+  geneta->SetYRange(fYMin/0.825, fYMax/0.825);
   AddSource2Generator(nameEta,geneta);
   TF1 *fPtEta = geneta->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
@@ -128,7 +131,8 @@ void AliGenEMCocktail::CreateCocktail()
   AliGenParam *genrho=0;
   Char_t nameRho[10];    
   snprintf(nameRho,10,"Rho");    
-  genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
+  genrho = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
+  genrho->SetYRange(fYMin/0.775, fYMax/0.775);
   AddSource2Generator(nameRho,genrho);
   TF1 *fPtRho = genrho->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
@@ -142,7 +146,8 @@ void AliGenEMCocktail::CreateCocktail()
   AliGenParam *genomega=0;
   Char_t nameOmega[10];    
   snprintf(nameOmega,10,"Omega");    
-  genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
+  genomega = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
+  genomega->SetYRange(fYMin/0.775, fYMax/0.775);
   AddSource2Generator(nameOmega,genomega);
   TF1 *fPtOmega = genomega->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
@@ -156,7 +161,8 @@ void AliGenEMCocktail::CreateCocktail()
   AliGenParam *genetaprime=0;
   Char_t nameEtaprime[10];    
   snprintf(nameEtaprime,10,"Etaprime");    
-  genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
+  genetaprime = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
+  genetaprime->SetYRange(fYMin/0.725, fYMax/0.725);
   AddSource2Generator(nameEtaprime,genetaprime);
   TF1 *fPtEtaprime = genetaprime->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
@@ -170,7 +176,8 @@ void AliGenEMCocktail::CreateCocktail()
   AliGenParam *genphi=0;
   Char_t namePhi[10];    
   snprintf(namePhi,10,"Phi");    
-  genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
+  genphi = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
+  genphi->SetYRange(fYMin/0.725, fYMax/0.725);
   AddSource2Generator(namePhi,genphi);
   TF1 *fPtPhi = genphi->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
@@ -184,7 +191,8 @@ void AliGenEMCocktail::CreateCocktail()
   AliGenParam *genjpsi=0;
   Char_t nameJpsi[10];    
   snprintf(nameJpsi,10,"Jpsi");    
-  genjpsi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
+  genjpsi = new AliGenParam(fNPart/0.525, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
+  genjpsi->SetYRange(fYMin/0.525, fYMax/0.525);
   AddSource2Generator(nameJpsi,genjpsi);
   TF1 *fPtJpsi = genjpsi->GetPt();
 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
@@ -193,6 +201,71 @@ void AliGenEMCocktail::CreateCocktail()
   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);
+#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();
+#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
+    fYieldArray[kGenPromptVirtGamma] = fPtPromptVirtG->Integral(fPtMin,fPtMax,1.e-6);
+#else
+    fYieldArray[kGenPromptVirtGamma] = fPtPromptVirtG->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 ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
+    fYieldArray[kGenThermVirtGamma] = fPtThermVirtG->Integral(fPtMin,fPtMax,1.e-6);
+#else
+    fYieldArray[kGenThermVirtGamma] = fPtThermVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+#endif
+  }
+  
 }
 
 //-------------------------------------------------------------------
@@ -204,7 +277,6 @@ void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
   Double_t phiMax = fPhiMax*180./TMath::Pi();
 
   genSource->SetPtRange(fPtMin, fPtMax);  
-  genSource->SetYRange(fYMin, fYMax);
   genSource->SetPhiRange(phiMin, phiMax);
   genSource->SetWeighting(fWeightingMode);
   genSource->SetForceGammaConversion(fForceConv);
@@ -264,7 +336,6 @@ void AliGenEMCocktail::Generate()
       if (igen == 1) entry->SetFirst(0);               
       else  entry->SetFirst((partArray->GetEntriesFast())+1);
       gen->SetEventPlane(evPlane);
-      gen->SetNumberParticles(fNPart);         
       gen->Generate();
       entry->SetLast(partArray->GetEntriesFast());
       preventry = entry;
@@ -288,6 +359,7 @@ void AliGenEMCocktail::Generate()
     }
     else
       pdgMother = part->GetPdgCode();
+
     switch (pdgMother){
     case 111:
       dNdy = fYieldArray[kGenPizero];
@@ -310,10 +382,24 @@ void AliGenEMCocktail::Generate()
     case 443:
       dNdy = fYieldArray[kGenJpsi];
       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 947506c..96ed249 100644 (file)
@@ -23,7 +23,7 @@ class AliGenEMCocktail : public AliGenCocktail
 public:
 
     AliGenEMCocktail();
-  enum GeneratorCode { kGenPizero=0, kGenEta, kGenRho, kGenOmega, kGenEtaprime, kGenPhi, kGenJpsi, kGENs };    
+  enum GeneratorCode { kGenPizero=0, kGenEta, kGenRho, kGenOmega, kGenEtaprime, kGenPhi, kGenJpsi, kGenPromptRealGamma, kGenThermRealGamma, kGenPromptVirtGamma, kGenThermVirtGamma, kGENs };
 
     virtual ~AliGenEMCocktail();    
     virtual void Init();
@@ -41,10 +41,6 @@ public:
   void    SetForceGammaConversion(Bool_t force=kTRUE){ fForceConv=force; }
   void    SetHeaviestParticle(Int_t part){ fHeaviestParticle=part; }
     
-protected:
-
-  TVector3& OrthogonalVector(TVector3 &inVec);
-    
 private:
     AliGenEMCocktail(const AliGenEMCocktail &cocktail); 
     AliGenEMCocktail & operator=(const AliGenEMCocktail &cocktail); 
index 5be02c9..9c80a67 100644 (file)
@@ -52,13 +52,13 @@ 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
   {  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
+  ,{ 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
   ,{ 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
@@ -69,10 +69,48 @@ const Double_t AliGenEMlib::fgkV2param[16][15] = {
   ,{ 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
 };
 
+const Double_t AliGenEMlib::fgkRawPtOfV2Param[16][10] = {
+  {  0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // pp no V2
+  ,{ 2.181446e+08, 9.412925e-01, 1.158774e-01, 3.020303e+01, 6.790828e+00, 9.999996e+01, 2.616827e+00, 3.980492e+00, 1.225169e+07, 5.575243e+00 } // 0-5
+  ,{ 3.006215e+08, 9.511881e-01, 1.192788e-01, 2.981931e+01, 5.068175e+01, 9.999993e+01, 2.650635e+00, 4.073982e+00, 2.508045e+07, 5.621039e+00 } // 5-10
+  ,{ 1.643438e+09, 9.604242e-01, 1.218512e-01, 2.912684e+01, 1.164242e+00, 9.999709e+01, 2.662326e+00, 4.027795e+00, 7.020810e+07, 5.696860e+00 } // 10-20
+  ,{ 8.109985e+08, 9.421935e-01, 1.328020e-01, 2.655910e+01, 1.053677e+00, 9.999812e+01, 2.722949e+00, 3.964547e+00, 6.104096e+07, 5.694703e+00 } // 20-30
+  ,{ 5.219789e+08, 9.417339e-01, 1.417541e-01, 2.518080e+01, 7.430803e-02, 9.303295e+01, 2.780227e+00, 3.909570e+00, 4.723116e+07, 5.778375e+00 } // 30-40
+  ,{ 2.547159e+08, 9.481459e-01, 2.364858e-01, 1.689288e+01, 3.858883e+00, 6.352619e+00, 2.742270e+00, 3.855226e+00, 3.120535e+07, 5.878677e+00 } // 40-50
+  ,{ 9.396097e+07, 9.304491e-01, 3.244940e-01, 1.291696e+01, 2.854367e+00, 6.325908e+00, 2.828258e+00, 4.133699e+00, 1.302739e+07, 5.977896e+00 } // 50-60
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-10 
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 20-40
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-60
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 60-80
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-20 
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-40 
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 20-80
+  ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-80
+};
+
+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
+  ,{ 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
+  ,{ 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 
+  ,{ 0.0000000000, 0.0000000000 } // 20-80
+  ,{ 0.0000000000, 0.0000000000 } // 40-80
+};
+
 // MASS   0=>PIZERO, 1=>ETA, 2=>RHO, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI, 6=>JPSI
-const Double_t AliGenEMlib::fgkHM[7] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946, 3.0969};
+const Double_t AliGenEMlib::fgkHM[8] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946, 3.0969, 0.0};
 
-const Double_t AliGenEMlib::fgkMtFactor[2][7] = { 
+const Double_t AliGenEMlib::fgkMtFactor[2][8] = { 
   // {1.0, 0.5, 1.0, 0.9, 0.4, 0.23, 0.054},  // factor for pp from arXiv:1110.3929
   // {1.0, 0.55, 1.0, 0.9, 0.4, 0.25, 0.004}    // factor for PbPb from arXiv:1110.3929
   //{1., 0.48, 1.0, 0.9, 0.25, 0.4}, (old values)
@@ -82,9 +120,10 @@ const Double_t AliGenEMlib::fgkMtFactor[2][7] = {
   //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.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
+  {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
 };
 
 //==========================================================================
@@ -162,7 +201,131 @@ Double_t AliGenEMlib::PtTsallis(const Double_t pt,
   return invYield*(2*TMath::Pi()*pt);
 }
 
+// Exponential
+Double_t AliGenEMlib::PtExponential(const Double_t *px, const Double_t *c){
+  const double &pt=px[0];
+  Double_t invYield = c[0]*exp(-pt*c[1]);
+  
+  return invYield*(2*TMath::Pi()*pt);
+}
+
+// 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]);
+  
+  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();
+}
+
+//--------------------------------------------------------------------------
+//
+//                             PromptRealGamma
+//
+//--------------------------------------------------------------------------
+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 };
+  return PtModifiedHagedornPowerlaw(px,promptGammaPtParam)*GetTAA(fgSelectedCentrality);
+}
+
+Double_t AliGenEMlib::YPromptRealGamma( 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;
+}
+
+//--------------------------------------------------------------------------
+//
+//                             PromptVirtGamma
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpPromptVirtGamma(TRandom *)
+{
+  return 223000;
+}
+
+Double_t AliGenEMlib::PtPromptVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+{
+  return IntegratedKrollWada(*px)*PtPromptRealGamma(px,px);
+}
+
+Double_t AliGenEMlib::YPromptVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+{
+  return YFlat(*px);
+}
+
+Double_t AliGenEMlib::V2PromptVirtGamma( const Double_t */*px*/, const Double_t */*dummy*/ )
+{
+  return 0.0;
+}
+
+//--------------------------------------------------------------------------
+//
+//                             ThermRealGamma
+//
+//--------------------------------------------------------------------------
+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*/ )
+{
+  return YFlat(*px);
+}
+
+Double_t AliGenEMlib::V2ThermRealGamma( const Double_t *px, const Double_t */*dummy*/ )
+{
+  return KEtScal(*px,8);
+}
+
+//--------------------------------------------------------------------------
+//
+//                             ThermVirtGamma
+//
+//--------------------------------------------------------------------------
+Int_t AliGenEMlib::IpThermVirtGamma(TRandom *)
+{
+  return 224000;
+}
+
+Double_t AliGenEMlib::PtThermVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+{
+  return IntegratedKrollWada(*px)*PtThermRealGamma(px,px);
+}
+
+Double_t AliGenEMlib::YThermVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+{
+  return YFlat(*px);
+}
 
+Double_t AliGenEMlib::V2ThermVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
+{
+  return KEtScal(*px,8);
+}
 
 //--------------------------------------------------------------------------
 //
@@ -177,6 +340,10 @@ 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
+  
   // 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.;
@@ -335,29 +502,20 @@ Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
 
 Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 {
-  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;
+    n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
+    v1=V2Param(px,fgkV2param[k0005]);
+    n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
+    v2=V2Param(px,fgkV2param[k0510]);
     return (n1*v1+n2*v2)/(n1+n2);
     break;
   case 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;
+    n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k2030]);
+    v1=V2Param(px,fgkV2param[k2030]);
+    n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k3040]);
+    v2=V2Param(px,fgkV2param[k3040]);
     return (n1*v1+n2*v2)/(n1+n2);
     break;
 
@@ -366,8 +524,6 @@ Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
   }
 }
 
-
-
 //--------------------------------------------------------------------------
 //
 //                              Eta
@@ -626,11 +782,12 @@ Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
 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 scaledPt = sqrt(pow(2.0/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);
+  return V2Pizero(&scaledPt, (Double_t*) 0);
 }
 
 Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
@@ -647,10 +804,28 @@ Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
 {
   // Flat v2
 
-  return 1.0;
+  return 0.0;
 }
 
-
+Double_t AliGenEMlib::GetTAA(Int_t cent){
+  const static Double_t taa[16] = { 1.0,    // pp
+                                   26.32,  // 0-5
+                                   20.56,  // 5-10
+                                   14.39,  // 10-20
+                                   8.70,   // 20-30
+                                   5.001,  // 30-40
+                                   2.675,  // 40-50
+                                   1.317,  // 50-60
+                                   23.44,  // 0-10
+                                   6.85,   // 20-40
+                                   1.996,  // 40-60
+                                   0.4174, // 60-80
+                                   18.91,  // 0-20
+                                   12.88,  // 0-40
+                                   3.088,  // 20-80
+                                   1.207}; // 40-80
+  return taa[cent];  
+}
 
 //==========================================================================
 //
@@ -670,6 +845,18 @@ GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
 
    switch (param) 
     {
+    case kPromptRealGamma:
+      func=PtPromptRealGamma;
+      break;
+    case kPromptVirtGamma:
+      func=PtPromptVirtGamma;
+      break;
+    case kThermRealGamma:
+      func=PtThermRealGamma;
+      break;
+    case kThermVirtGamma:
+      func=PtThermVirtGamma;
+      break;
     case kPizero:
       func=PtPizero;
       break;
@@ -707,6 +894,18 @@ GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
 
    switch (param) 
     {
+    case kPromptRealGamma:
+      func=YPromptRealGamma;
+      break;
+    case kPromptVirtGamma:
+      func=YPromptVirtGamma;
+      break;
+    case kThermRealGamma:
+      func=YThermRealGamma;
+      break;
+    case kThermVirtGamma:
+      func=YThermVirtGamma;
+      break;
     case kPizero:
          func=YPizero;
          break;
@@ -744,6 +943,18 @@ GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
 
    switch (param) 
     {
+    case kPromptRealGamma:
+      func=IpPromptRealGamma;
+      break;
+    case kPromptVirtGamma:
+      func=IpPromptVirtGamma;
+      break;
+    case kThermRealGamma:
+      func=IpThermRealGamma;
+      break;
+    case kThermVirtGamma:
+      func=IpThermVirtGamma;
+      break;
     case kPizero:
          func=IpPizero;
          break;
@@ -781,6 +992,18 @@ 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;
+      break;
+    case kThermVirtGamma:
+      func=V2ThermVirtGamma;
+      break;
     case kPizero:
       func=V2Pizero;
       break;
index 81c0864..6372507 100644 (file)
@@ -21,10 +21,10 @@ 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};
+  enum Particle_t{kPromptRealGamma, kPromptVirtGamma, kThermRealGamma, kThermVirtGamma, kPizero, kEta, kRho, kOmega, kEtaprime, kPhi, kJpsi};
+  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};
   
   AliGenEMlib() { } ;
 
@@ -79,8 +79,31 @@ public:
                            const Double_t T,
                            const Double_t n);
 
+  static Double_t PtExponential(const Double_t *pt, const Double_t *param);
+  static Double_t PtModifiedHagedornPowerlaw(const Double_t *pt, const Double_t *param);
+  static Double_t IntegratedKrollWada(Double_t mh);
 
+  // prompt gamma
+  static Int_t    IpPromptRealGamma(TRandom *ran);
+  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);
 
   // Pizero
     static Int_t    IpPizero(TRandom *ran);
@@ -131,13 +154,16 @@ public:
   static Double_t V2Param(const Double_t *px, const Double_t *param);
   static Double_t V2Flat(const Double_t *px, const Double_t *param);
   static Double_t KEtScal(Double_t pt, Int_t np);
+  static Double_t GetTAA(Int_t cent);
 
   static Double_t CrossOverLc(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 fgkV2param[16][15]; // array of v2 parameter
-  static const Double_t fgkHM[7];
-  static const Double_t fgkMtFactor[2][7];
+  static const Double_t fgkV2param[16][15];          // parameters of pi v2
+  static const Double_t fgkRawPtOfV2Param[16][10];   // parameters of the raw pt spectrum of v2 analysys
+  static const Double_t fgkThermPtParam[16][2];      // parameters of thermal gamma pt
+  static const Double_t fgkHM[8];                    // particle masses
+  static const Double_t fgkMtFactor[2][8];           // mt scaling factor
 
   ClassDef(AliGenEMlib,0)
 };
index be599f9..3e92b14 100644 (file)
@@ -215,83 +215,199 @@ TVector3 AliGenParam::OrthogonalVector(TVector3 &inVec){
   return res;
 }
 
-double AliGenParam::ScreenFunc1(double d){
-  if(d>1)
-    return 21.12-4.184*log(d+0.952);
+void AliGenParam::RotateVector(Double_t *pin, Double_t *pout, Double_t costheta, Double_t sintheta,
+    Double_t cosphi, Double_t sinphi)
+{
+  // Perform rotation
+  pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
+  pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
+  pout[2] = -1.0  * pin[0] * sintheta + pin[2] * 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);
   else
-    return 20.867-3.242*d+0.652*d*d;
+    return 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
 }
 
-double AliGenParam::ScreenFunc2(double d){
-  if(d>1)
-    return 21.12-4.184*log(d+0.952);
+double AliGenParam::ScreenFunction2(double screenVariable){
+  if(screenVariable>1)
+    return 42.24 - 8.368 * log(screenVariable + 0.952);
   else
-    return 20.209-1.93*d-0.086*d*d;
+    return 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
 }
 
-double AliGenParam::EnergyFraction(double Z, double E){
-  double e0=0.000511/E;
+double AliGenParam::RandomEnergyFraction(double Z, double photonEnergy){
   double aZ=Z/137.036;
-  double dmin=ScreenVar(Z,e0,0.5);
+  double epsilon ;
+  double epsilon0Local = 0.000511 / photonEnergy ;
+
+  // Do it fast if photon energy < 2. MeV
+  if (photonEnergy < 0.002 )
+    {
+      epsilon = epsilon0Local + (0.5 - epsilon0Local) * fRandom->Rndm();
+    }
+  else
+    {
+      double fZ = 8*log(Z)/3;
   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);
+      if (photonEnergy > 0.050) fZ += 8*fcZ;
+      
+      // Limits of the screening variable
+      double screenFactor = 136. * epsilon0Local / pow (Z,1/3);
+      double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
+      double screenMin = std::min(4.*screenFactor,screenMax) ;
+
+      // Limits of the energy sampling
+      double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
+      double epsilonMin = std::max(epsilon0Local,epsilon1);
+      double epsilonRange = 0.5 - epsilonMin ;
+
+      // Sample the energy rate of the created electron (or positron)
+      double screen;
+      double gReject ;
+
+      double f10 = ScreenFunction1(screenMin) - fZ;
+      double f20 = ScreenFunction2(screenMin) - fZ;
+      double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
+      double normF2 = std::max(1.5 * f20,0.);
+
+      do 
+       {
+         if (normF1 / (normF1 + normF2) > fRandom->Rndm() )
+           {
+             epsilon = 0.5 - epsilonRange * pow(fRandom->Rndm(), 0.333333) ;
+             screen = screenFactor / (epsilon * (1. - epsilon));
+             gReject = (ScreenFunction1(screen) - fZ) / f10 ;
+           }
+         else
+           {
+             epsilon = epsilonMin + epsilonRange * fRandom->Rndm();
+             screen = screenFactor / (epsilon * (1 - epsilon));
+             gReject = (ScreenFunction2(screen) - fZ) / f20 ;
+           }
+       } while ( gReject < fRandom->Rndm() );
+    
+    }   //  End of epsilon sampling
+  return epsilon;
+}
+
+double AliGenParam::RandomPolarAngle(){
+  double u;
+  const double a1 = 0.625;
+  double a2 = 3. * a1;
+  //  double d = 27. ;
+
+  //  if (9. / (9. + d) > fRandom->Rndm())
+  if (0.25 > fRandom->Rndm())
+    {
+      u = - log(fRandom->Rndm() * fRandom->Rndm()) / a1 ;
+    }
+  else
+    {
+      u = - log(fRandom->Rndm() * fRandom->Rndm()) / a2 ;
+    }
+  return u*0.000511;
+}
   
-  double normval=1/(0.5*(ScreenFunc1(dmin)-0.5*Fz)+0.1666667*(ScreenFunc2(dmin)-0.5*Fz));
+Double_t AliGenParam::RandomMass(Double_t mh){
   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 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 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;
   }
 }
 
-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::VirtualGammaPairProduction(TClonesArray *particles, Int_t nPart)
+{
+  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
+    double mass=RandomMass(gamma->Pt());
+
+    // lepton pair kinematics in virtual photon rest frame 
+    double Ee=mass/2;
+    double Pe=TMath::Sqrt((Ee+0.000511)*(Ee-0.000511));
+
+    double costheta = (2.0 * gRandom->Rndm()) - 1.;
+    double sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
+    double phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
+    double sinphi   = TMath::Sin(phi);
+    double cosphi   = TMath::Cos(phi);
+    
+    // momentum vectors of leptons in virtual photon rest frame
+    Double_t pProd1[3] = {Pe * sintheta * cosphi,
+                         Pe * sintheta * sinphi,
+                         Pe * costheta};
+
+    Double_t pProd2[3] = {-1.0 * Pe * sintheta * cosphi,
+                         -1.0 * Pe * sintheta * sinphi,
+                         -1.0 * Pe * costheta}; 
+
+    // lepton 4-vectors in properly rotated virtual photon rest frame
+    Double_t pRot1[3] = {0.};
+    RotateVector(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
+    Double_t pRot2[3] = {0.};
+    RotateVector(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi); 
+
+    TLorentzVector e1V4(pRot1[0],pRot1[1],pRot1[2],Ee);
+    TLorentzVector e2V4(pRot2[0],pRot2[1],pRot2[2],Ee);
+  
+    TVector3 boost(gamma->Px(),gamma->Py(),gamma->Pz());
+    boost*=1/sqrt(gamma->P()*gamma->P()+mass*mass);
+    e1V4.Boost(boost);
+    e2V4.Boost(boost);
+
+    TLorentzVector vtx;
+    gamma->ProductionVertex(vtx);
+    new((*particles)[nPartNew]) TParticle(11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, e1V4, vtx);
+    nPartNew++;
+    new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, e2V4, vtx);
+    nPartNew++;
+  }
 }
 
 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
+  //     and: G4LivermoreGammaConversionModel.cc
   Int_t nPartNew=nPart;
   for(int iPart=0; iPart<nPart; iPart++){      
     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());
-    Float_t az=fRandom->Uniform(TMath::Pi()*2);
-    double frac=EnergyFraction(1,gamma->Energy());
+    double frac=RandomEnergyFraction(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);
+    double Pe1=sqrt((Ee1+0.000511)*(Ee1-0.000511));
+    double Pe2=sqrt((Ee2+0.000511)*(Ee2-0.000511));
     
     TVector3 rotAxis(OrthogonalVector(gammaV3));
+    Float_t az=fRandom->Uniform(TMath::Pi()*2);
     rotAxis.Rotate(az,gammaV3);
     TVector3 e1V3(gammaV3);
-    e1V3.Rotate(PolarAngle(Ee1),rotAxis);
+    double u=RandomPolarAngle();
+    e1V3.Rotate(u/Ee1,rotAxis);
     e1V3=e1V3.Unit();
     e1V3*=Pe1;
     TVector3 e2V3(gammaV3);
-    e2V3.Rotate(-PolarAngle(Ee2),rotAxis);
+    e2V3.Rotate(-u/Ee2,rotAxis);
     e2V3=e2V3.Unit();
     e2V3*=Pe2;
     // gamma = new TParticle(*gamma);
@@ -460,6 +576,11 @@ void AliGenParam::GenerateN(Int_t ntimes)
       //
       // particle type 
          Int_t iPart = fIpParaFunc(fRandom);
+      Int_t iTemp = iPart;
+
+      // custom pdg codes for to destinguish direct photons
+      if((iPart>=221000) && (iPart<=229000)) iPart=22;
+
          fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;          
          TParticlePDG *particle = pDataBase->GetParticle(iPart);
          Float_t am = particle->Mass();
@@ -468,6 +589,7 @@ void AliGenParam::GenerateN(Int_t ntimes)
       //
       // y
          ty = TMath::TanH(fYPara->GetRandom());
+
       //
       // pT
          if (fAnalog == kAnalog) {
@@ -508,6 +630,7 @@ void AliGenParam::GenerateN(Int_t ntimes)
          p[0]=pt*TMath::Cos(phi);
          p[1]=pt*TMath::Sin(phi);
          p[2]=pl;
+
          if(fVertexSmear==kPerTrack) {
              Rndm(random,6);
              for (j=0;j<3;j++) {
@@ -540,13 +663,14 @@ void AliGenParam::GenerateN(Int_t ntimes)
        // select decay particles
              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)
+       }
 
-       // 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) 
index 4b20a26..522096a 100644 (file)
@@ -52,14 +52,16 @@ public:
     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);
+    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);
+    double RandomPolarAngle();
+    double RandomMass(Double_t mh);
+    Int_t VirtualGammaPairProduction(TClonesArray *particles, Int_t nPart);
     Int_t ForceGammaConversion(TClonesArray *particles, Int_t nPart);
   
 protected: