EVGEN/AliGenEMCocktail.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliGenEMCocktail.cxx 40702 2010-04-26 13:09:52Z morsch $ */
  17
  18 // Class to create the cocktail for physics with electrons, di-electrons,
  19 // and photons from the decay of the following sources:
  20 // pizero, eta, rho, omega, etaprime, phi
  21 // Kinematic distributions of the sources are taken from AliGenEMlib.
  22 // Decay channels can be selected via the method SetDecayMode.
  23 // Particles can be generated flat in pT with weights according to the
  24 // chosen pT distributions from AliGenEMlib (weighting mode: kNonAnalog),
  25 // or they are generated according to the pT distributions themselves
  26 // (weighting mode: kAnalog)
  27
  28
  29 #include <TObjArray.h>
  30 #include <TParticle.h>
  31 #include <TF1.h>
  32 #include <TVirtualMC.h>
  33 #include <TPDGCode.h>
  34 #include <TDatabasePDG.h>
  35 #include "AliGenCocktailEventHeader.h"
  36
  37 #include "AliGenCocktailEntry.h"
  38 #include "AliGenEMCocktail.h"
  39 #include "AliGenEMlib.h"
  40 #include "AliGenBox.h"
  41 #include "AliGenParam.h"
  42 #include "AliMC.h"
  43 #include "AliRun.h"
  44 #include "AliStack.h"
  45 #include "AliDecayer.h"
  46 #include "AliDecayerPythia.h"
  47 #include "AliLog.h"
  48 #include "AliGenCorrHF.h"
  49
  50 ClassImp(AliGenEMCocktail)
  51
  52 //________________________________________________________________________
  53 AliGenEMCocktail::AliGenEMCocktail()
  54 :AliGenCocktail(),
  55    fDecayer(0),
  56    fDecayMode(kAll),
  57    fWeightingMode(kNonAnalog),
  58    fNPart(1000),
  59   fYieldArray(),
  60   fPtSelect(AliGenEMlib::kPizero7TeVpp),
  61   fCentrality(AliGenEMlib::kpp),
  62   fV2Systematic(AliGenEMlib::kNoV2Sys),
  63   fForceConv(kFALSE),
  64   fSelectedParticles(kGenPizero)
  65 {
  66   // Constructor
  67   SetHeaviestHadron(kGenPhi);
  68 }
  69
  70 //_________________________________________________________________________
  71 AliGenEMCocktail::~AliGenEMCocktail()
  72 {
  73   // Destructor
  74
  75 }
  76
  77 //_________________________________________________________________________
  78 void AliGenEMCocktail::SetHeaviestHadron(ParticeGenerator_t part)
  79 {
  80   Int_t val=kGenPizero;
  81   while(val<part) val|=val<<1;
  82
  83   fSelectedParticles=val;
  84 }
  85
  86 //_________________________________________________________________________
  87 void AliGenEMCocktail::CreateCocktail()
  88 {
  89   // create and add sources to the cocktail
  90
  91   fDecayer->SetForceDecay(fDecayMode);
  92   fDecayer->ForceDecay();
  93
  94   // Set kinematic limits
  95   Double_t ptMin  = fPtMin;
  96   Double_t ptMax  = fPtMax;
  97   Double_t yMin   = fYMin;;
  98   Double_t yMax   = fYMax;;
  99   Double_t phiMin = fPhiMin*180./TMath::Pi();
 100   Double_t phiMax = fPhiMax*180./TMath::Pi();
 101   AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax));
 102   AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
 103
 104   //Initialize user selection for Pt Parameterization and centrality:
 105   AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
 106
 107   // Create and add electron sources to the generator
 108   // pizero
 109   if(fSelectedParticles&kGenPizero){
 110   AliGenParam *genpizero=0;
 111   Char_t namePizero[10];
 112   snprintf(namePizero,10,"Pizero");
 113     //fNPart/0.925: increase number of particles so that we have the chosen number of particles in the chosen eta range
 114   genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
 115     //fYMin/0.925: increase eta range, so that the electron yield is constant (<5% change) over the chosen eta range
 116   genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
 117   AddSource2Generator(namePizero,genpizero);
 118   TF1 *fPtPizero = genpizero->GetPt();
 119 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 120     fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
 121 #else
 122     fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 123 #endif
 124   }
 125
 126   // eta
 127   if(fSelectedParticles&kGenEta){
 128   AliGenParam *geneta=0;
 129   Char_t nameEta[10];
 130   snprintf(nameEta,10,"Eta");
 131   geneta = new AliGenParam(fNPart/0.825, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
 132   geneta->SetYRange(fYMin/0.825, fYMax/0.825);
 133   AddSource2Generator(nameEta,geneta);
 134   TF1 *fPtEta = geneta->GetPt();
 135 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 136     fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
 137 #else
 138     fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 139 #endif
 140   }
 141
 142   // rho
 143   if(fSelectedParticles&kGenRho){
 144   AliGenParam *genrho=0;
 145   Char_t nameRho[10];
 146   snprintf(nameRho,10,"Rho");
 147   genrho = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
 148   genrho->SetYRange(fYMin/0.775, fYMax/0.775);
 149   AddSource2Generator(nameRho,genrho);
 150   TF1 *fPtRho = genrho->GetPt();
 151 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 152     fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
 153 #else
 154     fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 155 #endif
 156   }
 157
 158   // omega
 159   if(fSelectedParticles&kGenOmega){
 160   AliGenParam *genomega=0;
 161   Char_t nameOmega[10];
 162   snprintf(nameOmega,10,"Omega");
 163   genomega = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
 164   genomega->SetYRange(fYMin/0.775, fYMax/0.775);
 165   AddSource2Generator(nameOmega,genomega);
 166   TF1 *fPtOmega = genomega->GetPt();
 167 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 168     fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
 169 #else
 170     fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 171 #endif
 172   }
 173
 174   // etaprime
 175   if(fSelectedParticles&kGenEtaprime){
 176   AliGenParam *genetaprime=0;
 177   Char_t nameEtaprime[10];
 178   snprintf(nameEtaprime,10,"Etaprime");
 179   genetaprime = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
 180   genetaprime->SetYRange(fYMin/0.725, fYMax/0.725);
 181   AddSource2Generator(nameEtaprime,genetaprime);
 182   TF1 *fPtEtaprime = genetaprime->GetPt();
 183 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 184     fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
 185 #else
 186     fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 187 #endif
 188   }
 189
 190   // phi
 191   if(fSelectedParticles&kGenPhi){
 192   AliGenParam *genphi=0;
 193   Char_t namePhi[10];
 194   snprintf(namePhi,10,"Phi");
 195   genphi = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
 196   genphi->SetYRange(fYMin/0.725, fYMax/0.725);
 197   AddSource2Generator(namePhi,genphi);
 198   TF1 *fPtPhi = genphi->GetPt();
 199 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 200     fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
 201 #else
 202     fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 203 #endif
 204   }
 205
 206   // jpsi
 207   if(fSelectedParticles&kGenJpsi){
 208   AliGenParam *genjpsi=0;
 209   Char_t nameJpsi[10];
 210   snprintf(nameJpsi,10,"Jpsi");
 211   genjpsi = new AliGenParam(fNPart/0.525, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
 212   genjpsi->SetYRange(fYMin/0.525, fYMax/0.525);
 213   AddSource2Generator(nameJpsi,genjpsi);
 214   TF1 *fPtJpsi = genjpsi->GetPt();
 215 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 216     fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
 217 #else
 218     fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 219 #endif
 220   }
 221
 222   // direct gamma
 223   if(fDecayMode!=kGammaEM) return;
 224   if(fSelectedParticles&kGenDirectRealGamma){
 225     AliGenParam *genDirectRealG=0;
 226     Char_t nameDirectRealG[10];
 227     snprintf(nameDirectRealG,10,"DirectRealGamma");
 228     genDirectRealG = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDirectRealGamma, "DUMMY");
 229     genDirectRealG->SetYRange(fYMin, fYMax);
 230     AddSource2Generator(nameDirectRealG,genDirectRealG);
 231     TF1 *fPtDirectRealG = genDirectRealG->GetPt();
 232 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 233     fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
 234 #else
 235     fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 236 #endif
 237   }
 238
 239   if(fSelectedParticles&kGenDirectVirtGamma){
 240     TDatabasePDG::Instance()->AddParticle("DirectVirtGamma","DirectVirtGamma",0,true,0,0,"GaugeBoson",220000);
 241     AliGenParam *genDirectVirtG=0;
 242     Char_t nameDirectVirtG[10];
 243     snprintf(nameDirectVirtG,10,"DirectVirtGamma");
 244     genDirectVirtG = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kDirectVirtGamma, "DUMMY");
 245     genDirectVirtG->SetYRange(fYMin/0.775, fYMax/0.775);
 246     AddSource2Generator(nameDirectVirtG,genDirectVirtG);
 247     TF1 *fPtDirectVirtG = genDirectVirtG->GetPt();
 248 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 249     fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,1.e-6);
 250 #else
 251     fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 252 #endif
 253   }
 254 }
 255
 256 //-------------------------------------------------------------------
 257 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
 258                                          AliGenParam* const genSource)
 259 {
 260 // add sources to the cocktail
 261   Double_t phiMin = fPhiMin*180./TMath::Pi();
 262   Double_t phiMax = fPhiMax*180./TMath::Pi();
 263
 264   genSource->SetPtRange(fPtMin, fPtMax);
 265   genSource->SetPhiRange(phiMin, phiMax);
 266   genSource->SetWeighting(fWeightingMode);
 267   genSource->SetForceGammaConversion(fForceConv);
 268   if (!TVirtualMC::GetMC()) genSource->SetDecayer(fDecayer);
 269   genSource->Init();
 270
 271   AddGenerator(genSource,nameSource,1.); // Adding Generator
 272 }
 273
 274 //-------------------------------------------------------------------
 275 void AliGenEMCocktail::Init()
 276 {
 277   // Initialisation
 278   TIter next(fEntries);
 279   AliGenCocktailEntry *entry;
 280   if (fStack) {
 281     while((entry = (AliGenCocktailEntry*)next())) {
 282       entry->Generator()->SetStack(fStack);
 283     }
 284   }
 285 }
 286
 287 //_________________________________________________________________________
 288 void AliGenEMCocktail::Generate()
 289 {
 290   // Generate event
 291   TIter next(fEntries);
 292   AliGenCocktailEntry *entry = 0;
 293   AliGenCocktailEntry *preventry = 0;
 294   AliGenerator* gen = 0;
 295
 296   if (fHeader) delete fHeader;
 297   fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
 298
 299   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 300
 301   // Generate the vertex position used by all generators
 302   if(fVertexSmear == kPerEvent) Vertex();
 303
 304   //Reseting stack
 305   AliRunLoader * runloader = AliRunLoader::Instance();
 306   if (runloader)
 307     if (runloader->Stack())
 308       runloader->Stack()->Clean();
 309
 310   // Loop over generators and generate events
 311   Int_t igen = 0;
 312   Float_t evPlane;
 313   Rndm(&evPlane,1);
 314   evPlane*=TMath::Pi()*2;
 315   while((entry = (AliGenCocktailEntry*)next())) {
 316     gen = entry->Generator();
 317     gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 318
 319     if (fNPart > 0) {
 320       igen++;
 321       if (igen == 1) entry->SetFirst(0);
 322       else  entry->SetFirst((partArray->GetEntriesFast())+1);
 323       gen->SetEventPlane(evPlane);
 324       gen->Generate();
 325       entry->SetLast(partArray->GetEntriesFast());
 326       preventry = entry;
 327     }
 328   }
 329   next.Reset();
 330
 331   // Setting weights for proper absolute normalization
 332   Int_t iPart, iMother;
 333   Int_t pdgMother = 0;
 334   Double_t weight = 0.;
 335   Double_t dNdy = 0.;
 336   Int_t maxPart = partArray->GetEntriesFast();
 337   for(iPart=0; iPart<maxPart; iPart++){
 338     TParticle *part = gAlice->GetMCApp()->Particle(iPart);
 339     iMother = part->GetFirstMother();
 340     TParticle *mother = 0;
 341     if (iMother>=0){
 342       mother = gAlice->GetMCApp()->Particle(iMother);
 343       pdgMother = mother->GetPdgCode();
 344     }
 345     else
 346       pdgMother = part->GetPdgCode();
 347
 348     switch (pdgMother){
 349     case 111:
 350       dNdy = fYieldArray[kPizero];
 351       break;
 352     case 221:
 353       dNdy = fYieldArray[kEta];
 354       break;
 355     case 113:
 356       dNdy = fYieldArray[kRho];
 357       break;
 358     case 223:
 359       dNdy = fYieldArray[kOmega];
 360       break;
 361     case 331:
 362       dNdy = fYieldArray[kEtaprime];
 363       break;
 364     case 333:
 365       dNdy = fYieldArray[kPhi];
 366       break;
 367     case 443:
 368       dNdy = fYieldArray[kJpsi];
 369       break;
 370     case 22:
 371       dNdy = fYieldArray[kDirectRealGamma];
 372       break;
 373     case 220000:
 374       dNdy = fYieldArray[kDirectVirtGamma];
 375       break;
 376
 377     default:
 378       dNdy = 0.;
 379     }
 380
 381     weight = dNdy*part->GetWeight();
 382     part->SetWeight(weight);
 383   }
 384
 385   fHeader->SetNProduced(maxPart);
 386
 387
 388   TArrayF eventVertex;
 389   eventVertex.Set(3);
 390   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
 391
 392   fHeader->SetPrimaryVertex(eventVertex);
 393
 394   gAlice->SetGenEventHeader(fHeader);
 395 }