]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenEMCocktail.cxx
- update to newest fit parameters
[u/mrichter/AliRoot.git] / EVGEN / AliGenEMCocktail.cxx
index 5bb28dad48174469e781feb355f0c7d84187eb5c..2db7bfdd04ed041d11bd02fcf6cfca9b4a52705a 100644 (file)
@@ -31,6 +31,7 @@
 #include <TF1.h>
 #include <TVirtualMC.h>
 #include <TPDGCode.h>
+#include <TDatabasePDG.h>
 #include "AliGenCocktailEventHeader.h"
 
 #include "AliGenCocktailEntry.h"
@@ -56,14 +57,14 @@ AliGenEMCocktail::AliGenEMCocktail()
    fWeightingMode(kNonAnalog),
    fNPart(1000),
   fYieldArray(),
-  fPtSelect(0),
-  fCentrality(0),
-  fV2Systematic(0),
+  fPtSelect(AliGenEMlib::kPizero7TeVpp),
+  fCentrality(AliGenEMlib::kpp),
+  fV2Systematic(AliGenEMlib::kNoV2Sys),
   fForceConv(kFALSE),
-  fHeaviestParticle(kGENs)
+  fSelectedParticles(kGenPizero)
 {
   // Constructor
-
+  SetHeaviestHadron(kGenPhi);
 }
 
 //_________________________________________________________________________
@@ -73,6 +74,15 @@ AliGenEMCocktail::~AliGenEMCocktail()
 
 }
 
+//_________________________________________________________________________
+void AliGenEMCocktail::SetHeaviestHadron(ParticeGenerator_t part)
+{
+  Int_t val=kGenPizero;
+  while(val<part) val|=val<<1;
+
+  fSelectedParticles=val;
+}
+
 //_________________________________________________________________________
 void AliGenEMCocktail::CreateCocktail()
 {
@@ -95,104 +105,152 @@ void AliGenEMCocktail::CreateCocktail()
   AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
 
   // Create and add electron sources to the generator
-
   // pizero
+  if(fSelectedParticles&kGenPizero){
   AliGenParam *genpizero=0;
   Char_t namePizero[10];    
   snprintf(namePizero,10,"Pizero");    
-  genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
+    //fNPart/0.925: increase number of particles so that we have the chosen number of particles in the chosen eta range
+  genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
+    //fYMin/0.925: increase eta range, so that the electron yield is constant (<5% change) over the chosen eta range
+  genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
   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");    
-  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)
-  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");    
-  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)
-  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");    
-  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)
-  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");    
-  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)
-  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");    
-  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)
-  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");    
-  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)
-  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);
+    fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 #endif
+  }
 
+  // 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[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
+#else
+    fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+#endif
+  }
+
+  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[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,1.e-6);
+#else
+    fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
+#endif
+  }
 }
 
 //-------------------------------------------------------------------
@@ -204,11 +262,10 @@ 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);
-  if (!gMC) genSource->SetDecayer(fDecayer);  
+  if (!TVirtualMC::GetMC()) genSource->SetDecayer(fDecayer);  
   genSource->Init();
     
   AddGenerator(genSource,nameSource,1.); // Adding Generator    
@@ -264,7 +321,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,32 +344,40 @@ void AliGenEMCocktail::Generate()
     }
     else
       pdgMother = part->GetPdgCode();
+
     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.;
     }
+
     weight = dNdy*part->GetWeight();
     part->SetWeight(weight);
   }