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(0),
  61   fCentrality(0),
  62   fV2Systematic(0),
  63   fForceConv(kFALSE),
  64   fHeaviestParticle(kGenJpsi)
  65 {
  66   // Constructor
  67
  68 }
  69
  70 //_________________________________________________________________________
  71 AliGenEMCocktail::~AliGenEMCocktail()
  72 {
  73   // Destructor
  74
  75 }
  76
  77 //_________________________________________________________________________
  78 void AliGenEMCocktail::CreateCocktail()
  79 {
  80   // create and add sources to the cocktail
  81
  82   fDecayer->SetForceDecay(fDecayMode);
  83   fDecayer->ForceDecay();
  84
  85   // Set kinematic limits
  86   Double_t ptMin  = fPtMin;
  87   Double_t ptMax  = fPtMax;
  88   Double_t yMin   = fYMin;;
  89   Double_t yMax   = fYMax;;
  90   Double_t phiMin = fPhiMin*180./TMath::Pi();
  91   Double_t phiMax = fPhiMax*180./TMath::Pi();
  92   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));
  93   AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
  94
  95   //Initialize user selection for Pt Parameterization and centrality:
  96   AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
  97
  98   // Create and add electron sources to the generator
  99   // pizero
 100   if(fHeaviestParticle<kGenPizero)return;
 101   AliGenParam *genpizero=0;
 102   Char_t namePizero[10];
 103   snprintf(namePizero,10,"Pizero");
 104   genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
 105   genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
 106   AddSource2Generator(namePizero,genpizero);
 107   TF1 *fPtPizero = genpizero->GetPt();
 108 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 109   fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
 110 #else
 111   fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 112 #endif
 113
 114   // eta
 115   if(fHeaviestParticle<kGenEta)return;
 116   AliGenParam *geneta=0;
 117   Char_t nameEta[10];
 118   snprintf(nameEta,10,"Eta");
 119   geneta = new AliGenParam(fNPart/0.825, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
 120   geneta->SetYRange(fYMin/0.825, fYMax/0.825);
 121   AddSource2Generator(nameEta,geneta);
 122   TF1 *fPtEta = geneta->GetPt();
 123 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 124   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
 125 #else
 126   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 127 #endif
 128
 129   // rho
 130   if(fHeaviestParticle<kGenRho)return;
 131   AliGenParam *genrho=0;
 132   Char_t nameRho[10];
 133   snprintf(nameRho,10,"Rho");
 134   genrho = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
 135   genrho->SetYRange(fYMin/0.775, fYMax/0.775);
 136   AddSource2Generator(nameRho,genrho);
 137   TF1 *fPtRho = genrho->GetPt();
 138 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 139   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
 140 #else
 141   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 142 #endif
 143
 144   // omega
 145   if(fHeaviestParticle<kGenOmega)return;
 146   AliGenParam *genomega=0;
 147   Char_t nameOmega[10];
 148   snprintf(nameOmega,10,"Omega");
 149   genomega = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
 150   genomega->SetYRange(fYMin/0.775, fYMax/0.775);
 151   AddSource2Generator(nameOmega,genomega);
 152   TF1 *fPtOmega = genomega->GetPt();
 153 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 154   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
 155 #else
 156   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 157 #endif
 158
 159   // etaprime
 160   if(fHeaviestParticle<kGenEtaprime)return;
 161   AliGenParam *genetaprime=0;
 162   Char_t nameEtaprime[10];
 163   snprintf(nameEtaprime,10,"Etaprime");
 164   genetaprime = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
 165   genetaprime->SetYRange(fYMin/0.725, fYMax/0.725);
 166   AddSource2Generator(nameEtaprime,genetaprime);
 167   TF1 *fPtEtaprime = genetaprime->GetPt();
 168 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 169   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
 170 #else
 171   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 172 #endif
 173
 174   // phi
 175   if(fHeaviestParticle<kGenPhi)return;
 176   AliGenParam *genphi=0;
 177   Char_t namePhi[10];
 178   snprintf(namePhi,10,"Phi");
 179   genphi = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
 180   genphi->SetYRange(fYMin/0.725, fYMax/0.725);
 181   AddSource2Generator(namePhi,genphi);
 182   TF1 *fPtPhi = genphi->GetPt();
 183 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 184   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
 185 #else
 186   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 187 #endif
 188
 189   // jpsi
 190   if(fHeaviestParticle<kGenJpsi)return;
 191   AliGenParam *genjpsi=0;
 192   Char_t nameJpsi[10];
 193   snprintf(nameJpsi,10,"Jpsi");
 194   genjpsi = new AliGenParam(fNPart/0.525, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
 195   genjpsi->SetYRange(fYMin/0.525, fYMax/0.525);
 196   AddSource2Generator(nameJpsi,genjpsi);
 197   TF1 *fPtJpsi = genjpsi->GetPt();
 198 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 199   fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
 200 #else
 201   fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 202 #endif
 203
 204   // prompt gamma
 205   if(fDecayMode==kGammaEM){
 206     if(fHeaviestParticle<kGenPromptRealGamma)return;
 207     TDatabasePDG::Instance()->AddParticle("PromptRealGamma","PromptRealGamma",0,true,0,0,"GaugeBoson",221000);
 208     //gMC->DefineParticle(221000, "PromptGamma", kPTGamma, 0, 0, 0,"Gamma", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kFALSE);
 209     AliGenParam *genPromptRealG=0;
 210     Char_t namePromptRealG[10];
 211     snprintf(namePromptRealG,10,"PromptRealGamma");
 212     genPromptRealG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kPromptRealGamma, "DUMMY");
 213     genPromptRealG->SetYRange(fYMin, fYMax);
 214     AddSource2Generator(namePromptRealG,genPromptRealG);
 215     TF1 *fPtPromptRealG = genPromptRealG->GetPt();
 216 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 217     fYieldArray[kGenPromptRealGamma] = fPtPromptRealG->Integral(fPtMin,fPtMax,1.e-6);
 218 #else
 219     fYieldArray[kGenPromptRealGamma] = fPtPromptRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 220 #endif
 221
 222     if(fHeaviestParticle<kGenThermRealGamma)return;
 223     TDatabasePDG::Instance()->AddParticle("ThermRealGamma","ThermRealGamma",0,true,0,0,"GaugeBoson",222000);
 224     //gMC->DefineParticle(221000, "ThermGamma", kPTGamma, 0, 0, 0,"Gamma", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kFALSE);
 225     AliGenParam *genThermRealG=0;
 226     Char_t nameThermRealG[10];
 227     snprintf(nameThermRealG,10,"ThermRealGamma");
 228     genThermRealG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kThermRealGamma, "DUMMY");
 229     genThermRealG->SetYRange(fYMin, fYMax);
 230     AddSource2Generator(nameThermRealG,genThermRealG);
 231     TF1 *fPtThermRealG = genThermRealG->GetPt();
 232 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 233     fYieldArray[kGenThermRealGamma] = fPtThermRealG->Integral(fPtMin,fPtMax,1.e-6);
 234 #else
 235     fYieldArray[kGenThermRealGamma] = fPtThermRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 236 #endif
 237
 238     if(fHeaviestParticle<kGenPromptVirtGamma)return;
 239     TDatabasePDG::Instance()->AddParticle("PromptVirtGamma","PromptVirtGamma",0,true,0,0,"GaugeBoson",223000);
 240     AliGenParam *genPromptVirtG=0;
 241     Char_t namePromptVirtG[10];
 242     snprintf(namePromptVirtG,10,"PromptVirtGamma");
 243     genPromptVirtG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kPromptVirtGamma, "DUMMY");
 244     genPromptVirtG->SetYRange(fYMin, fYMax);
 245     AddSource2Generator(namePromptVirtG,genPromptVirtG);
 246     TF1 *fPtPromptVirtG = genPromptVirtG->GetPt();
 247 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 248     fYieldArray[kGenPromptVirtGamma] = fPtPromptVirtG->Integral(fPtMin,fPtMax,1.e-6);
 249 #else
 250     fYieldArray[kGenPromptVirtGamma] = fPtPromptVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 251 #endif
 252
 253     if(fHeaviestParticle<kGenThermVirtGamma)return;
 254     TDatabasePDG::Instance()->AddParticle("ThermVirtGamma","ThermVirtGamma",0,true,0,0,"GaugeBoson",224000);
 255     AliGenParam *genThermVirtG=0;
 256     Char_t nameThermVirtG[10];
 257     snprintf(nameThermVirtG,10,"ThermVirtGamma");
 258     genThermVirtG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kThermVirtGamma, "DUMMY");
 259     genThermVirtG->SetYRange(fYMin, fYMax);
 260     AddSource2Generator(nameThermVirtG,genThermVirtG);
 261     TF1 *fPtThermVirtG = genThermVirtG->GetPt();
 262 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 263     fYieldArray[kGenThermVirtGamma] = fPtThermVirtG->Integral(fPtMin,fPtMax,1.e-6);
 264 #else
 265     fYieldArray[kGenThermVirtGamma] = fPtThermVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 266 #endif
 267   }
 268
 269 }
 270
 271 //-------------------------------------------------------------------
 272 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
 273                                          AliGenParam* const genSource)
 274 {
 275 // add sources to the cocktail
 276   Double_t phiMin = fPhiMin*180./TMath::Pi();
 277   Double_t phiMax = fPhiMax*180./TMath::Pi();
 278
 279   genSource->SetPtRange(fPtMin, fPtMax);
 280   genSource->SetPhiRange(phiMin, phiMax);
 281   genSource->SetWeighting(fWeightingMode);
 282   genSource->SetForceGammaConversion(fForceConv);
 283   if (!gMC) genSource->SetDecayer(fDecayer);
 284   genSource->Init();
 285
 286   AddGenerator(genSource,nameSource,1.); // Adding Generator
 287 }
 288
 289 //-------------------------------------------------------------------
 290 void AliGenEMCocktail::Init()
 291 {
 292   // Initialisation
 293   TIter next(fEntries);
 294   AliGenCocktailEntry *entry;
 295   if (fStack) {
 296     while((entry = (AliGenCocktailEntry*)next())) {
 297       entry->Generator()->SetStack(fStack);
 298     }
 299   }
 300 }
 301
 302 //_________________________________________________________________________
 303 void AliGenEMCocktail::Generate()
 304 {
 305   // Generate event
 306   TIter next(fEntries);
 307   AliGenCocktailEntry *entry = 0;
 308   AliGenCocktailEntry *preventry = 0;
 309   AliGenerator* gen = 0;
 310
 311   if (fHeader) delete fHeader;
 312   fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
 313
 314   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 315
 316   // Generate the vertex position used by all generators
 317   if(fVertexSmear == kPerEvent) Vertex();
 318
 319   //Reseting stack
 320   AliRunLoader * runloader = AliRunLoader::Instance();
 321   if (runloader)
 322     if (runloader->Stack())
 323       runloader->Stack()->Clean();
 324
 325   // Loop over generators and generate events
 326   Int_t igen = 0;
 327   Float_t evPlane;
 328   Rndm(&evPlane,1);
 329   evPlane*=TMath::Pi()*2;
 330   while((entry = (AliGenCocktailEntry*)next())) {
 331     gen = entry->Generator();
 332     gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 333
 334     if (fNPart > 0) {
 335       igen++;
 336       if (igen == 1) entry->SetFirst(0);
 337       else  entry->SetFirst((partArray->GetEntriesFast())+1);
 338       gen->SetEventPlane(evPlane);
 339       gen->Generate();
 340       entry->SetLast(partArray->GetEntriesFast());
 341       preventry = entry;
 342     }
 343   }
 344   next.Reset();
 345
 346   // Setting weights for proper absolute normalization
 347   Int_t iPart, iMother;
 348   Int_t pdgMother = 0;
 349   Double_t weight = 0.;
 350   Double_t dNdy = 0.;
 351   Int_t maxPart = partArray->GetEntriesFast();
 352   for(iPart=0; iPart<maxPart; iPart++){
 353     TParticle *part = gAlice->GetMCApp()->Particle(iPart);
 354     iMother = part->GetFirstMother();
 355     TParticle *mother = 0;
 356     if (iMother>=0){
 357       mother = gAlice->GetMCApp()->Particle(iMother);
 358       pdgMother = mother->GetPdgCode();
 359     }
 360     else
 361       pdgMother = part->GetPdgCode();
 362
 363     switch (pdgMother){
 364     case 111:
 365       dNdy = fYieldArray[kGenPizero];
 366       break;
 367     case 221:
 368       dNdy = fYieldArray[kGenEta];
 369       break;
 370     case 113:
 371       dNdy = fYieldArray[kGenRho];
 372       break;
 373     case 223:
 374       dNdy = fYieldArray[kGenOmega];
 375       break;
 376     case 331:
 377       dNdy = fYieldArray[kGenEtaprime];
 378       break;
 379     case 333:
 380       dNdy = fYieldArray[kGenPhi];
 381       break;
 382     case 443:
 383       dNdy = fYieldArray[kGenJpsi];
 384       break;
 385     default:
 386       dNdy = 0.;
 387     }
 388
 389     switch (pdgMother){
 390     case 221000:
 391       dNdy = fYieldArray[kGenPromptRealGamma];
 392       break;
 393     case 223000:
 394       dNdy = fYieldArray[kGenPromptVirtGamma];
 395       break;
 396     case 222000:
 397       dNdy = fYieldArray[kGenThermRealGamma];
 398       break;
 399     case 224000:
 400       dNdy = fYieldArray[kGenThermVirtGamma];
 401       break;
 402     }
 403     weight = dNdy*part->GetWeight();
 404     part->SetWeight(weight);
 405   }
 406
 407   fHeader->SetNProduced(maxPart);
 408
 409
 410   TArrayF eventVertex;
 411   eventVertex.Set(3);
 412   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
 413
 414   fHeader->SetPrimaryVertex(eventVertex);
 415
 416   gAlice->SetGenEventHeader(fHeader);
 417 }