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   fSelectedParticles(kGenHadrons)
  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(fSelectedParticles&kGenPizero){
 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[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
 110 #else
 111     fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 112 #endif
 113   }
 114
 115   // eta
 116   if(fSelectedParticles&kGenEta){
 117   AliGenParam *geneta=0;
 118   Char_t nameEta[10];
 119   snprintf(nameEta,10,"Eta");
 120   geneta = new AliGenParam(fNPart/0.825, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
 121   geneta->SetYRange(fYMin/0.825, fYMax/0.825);
 122   AddSource2Generator(nameEta,geneta);
 123   TF1 *fPtEta = geneta->GetPt();
 124 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 125     fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
 126 #else
 127     fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 128 #endif
 129   }
 130
 131   // rho
 132   if(fSelectedParticles&kGenRho){
 133   AliGenParam *genrho=0;
 134   Char_t nameRho[10];
 135   snprintf(nameRho,10,"Rho");
 136   genrho = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
 137   genrho->SetYRange(fYMin/0.775, fYMax/0.775);
 138   AddSource2Generator(nameRho,genrho);
 139   TF1 *fPtRho = genrho->GetPt();
 140 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 141     fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
 142 #else
 143     fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 144 #endif
 145   }
 146
 147   // omega
 148   if(fSelectedParticles&kGenOmega){
 149   AliGenParam *genomega=0;
 150   Char_t nameOmega[10];
 151   snprintf(nameOmega,10,"Omega");
 152   genomega = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
 153   genomega->SetYRange(fYMin/0.775, fYMax/0.775);
 154   AddSource2Generator(nameOmega,genomega);
 155   TF1 *fPtOmega = genomega->GetPt();
 156 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 157     fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
 158 #else
 159     fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 160 #endif
 161   }
 162
 163   // etaprime
 164   if(fSelectedParticles&kGenEtaprime){
 165   AliGenParam *genetaprime=0;
 166   Char_t nameEtaprime[10];
 167   snprintf(nameEtaprime,10,"Etaprime");
 168   genetaprime = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
 169   genetaprime->SetYRange(fYMin/0.725, fYMax/0.725);
 170   AddSource2Generator(nameEtaprime,genetaprime);
 171   TF1 *fPtEtaprime = genetaprime->GetPt();
 172 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 173     fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
 174 #else
 175     fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 176 #endif
 177   }
 178
 179   // phi
 180   if(fSelectedParticles&kGenPhi){
 181   AliGenParam *genphi=0;
 182   Char_t namePhi[10];
 183   snprintf(namePhi,10,"Phi");
 184   genphi = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
 185   genphi->SetYRange(fYMin/0.725, fYMax/0.725);
 186   AddSource2Generator(namePhi,genphi);
 187   TF1 *fPtPhi = genphi->GetPt();
 188 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 189     fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
 190 #else
 191     fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 192 #endif
 193   }
 194
 195   // jpsi
 196   if(fSelectedParticles&kGenJpsi){
 197   AliGenParam *genjpsi=0;
 198   Char_t nameJpsi[10];
 199   snprintf(nameJpsi,10,"Jpsi");
 200   genjpsi = new AliGenParam(fNPart/0.525, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
 201   genjpsi->SetYRange(fYMin/0.525, fYMax/0.525);
 202   AddSource2Generator(nameJpsi,genjpsi);
 203   TF1 *fPtJpsi = genjpsi->GetPt();
 204 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 205     fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
 206 #else
 207     fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 208 #endif
 209   }
 210
 211   // direct gamma
 212   if(fDecayMode!=kGammaEM) return;
 213   if(fSelectedParticles&kGenDirectRealGamma){
 214     AliGenParam *genDirectRealG=0;
 215     Char_t nameDirectRealG[10];
 216     snprintf(nameDirectRealG,10,"DirectRealGamma");
 217     genDirectRealG = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDirectRealGamma, "DUMMY");
 218     genDirectRealG->SetYRange(fYMin, fYMax);
 219     AddSource2Generator(nameDirectRealG,genDirectRealG);
 220     TF1 *fPtDirectRealG = genDirectRealG->GetPt();
 221 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 222     fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
 223 #else
 224     fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 225 #endif
 226   }
 227
 228   if(fSelectedParticles&kGenDirectVirtGamma){
 229     TDatabasePDG::Instance()->AddParticle("DirectVirtGamma","DirectVirtGamma",0,true,0,0,"GaugeBoson",220000);
 230     AliGenParam *genDirectVirtG=0;
 231     Char_t nameDirectVirtG[10];
 232     snprintf(nameDirectVirtG,10,"DirectVirtGamma");
 233     genDirectVirtG = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kDirectVirtGamma, "DUMMY");
 234     genDirectVirtG->SetYRange(fYMin/0.775, fYMax/0.775);
 235     AddSource2Generator(nameDirectVirtG,genDirectVirtG);
 236     TF1 *fPtDirectVirtG = genDirectVirtG->GetPt();
 237 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 238     fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,1.e-6);
 239 #else
 240     fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 241 #endif
 242   }
 243 }
 244
 245 //-------------------------------------------------------------------
 246 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
 247                                          AliGenParam* const genSource)
 248 {
 249 // add sources to the cocktail
 250   Double_t phiMin = fPhiMin*180./TMath::Pi();
 251   Double_t phiMax = fPhiMax*180./TMath::Pi();
 252
 253   genSource->SetPtRange(fPtMin, fPtMax);
 254   genSource->SetPhiRange(phiMin, phiMax);
 255   genSource->SetWeighting(fWeightingMode);
 256   genSource->SetForceGammaConversion(fForceConv);
 257   if (!TVirtualMC::GetMC()) genSource->SetDecayer(fDecayer);
 258   genSource->Init();
 259
 260   AddGenerator(genSource,nameSource,1.); // Adding Generator
 261 }
 262
 263 //-------------------------------------------------------------------
 264 void AliGenEMCocktail::Init()
 265 {
 266   // Initialisation
 267   TIter next(fEntries);
 268   AliGenCocktailEntry *entry;
 269   if (fStack) {
 270     while((entry = (AliGenCocktailEntry*)next())) {
 271       entry->Generator()->SetStack(fStack);
 272     }
 273   }
 274 }
 275
 276 //_________________________________________________________________________
 277 void AliGenEMCocktail::Generate()
 278 {
 279   // Generate event
 280   TIter next(fEntries);
 281   AliGenCocktailEntry *entry = 0;
 282   AliGenCocktailEntry *preventry = 0;
 283   AliGenerator* gen = 0;
 284
 285   if (fHeader) delete fHeader;
 286   fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
 287
 288   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 289
 290   // Generate the vertex position used by all generators
 291   if(fVertexSmear == kPerEvent) Vertex();
 292
 293   //Reseting stack
 294   AliRunLoader * runloader = AliRunLoader::Instance();
 295   if (runloader)
 296     if (runloader->Stack())
 297       runloader->Stack()->Clean();
 298
 299   // Loop over generators and generate events
 300   Int_t igen = 0;
 301   Float_t evPlane;
 302   Rndm(&evPlane,1);
 303   evPlane*=TMath::Pi()*2;
 304   while((entry = (AliGenCocktailEntry*)next())) {
 305     gen = entry->Generator();
 306     gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 307
 308     if (fNPart > 0) {
 309       igen++;
 310       if (igen == 1) entry->SetFirst(0);
 311       else  entry->SetFirst((partArray->GetEntriesFast())+1);
 312       gen->SetEventPlane(evPlane);
 313       gen->Generate();
 314       entry->SetLast(partArray->GetEntriesFast());
 315       preventry = entry;
 316     }
 317   }
 318   next.Reset();
 319
 320   // Setting weights for proper absolute normalization
 321   Int_t iPart, iMother;
 322   Int_t pdgMother = 0;
 323   Double_t weight = 0.;
 324   Double_t dNdy = 0.;
 325   Int_t maxPart = partArray->GetEntriesFast();
 326   for(iPart=0; iPart<maxPart; iPart++){
 327     TParticle *part = gAlice->GetMCApp()->Particle(iPart);
 328     iMother = part->GetFirstMother();
 329     TParticle *mother = 0;
 330     if (iMother>=0){
 331       mother = gAlice->GetMCApp()->Particle(iMother);
 332       pdgMother = mother->GetPdgCode();
 333     }
 334     else
 335       pdgMother = part->GetPdgCode();
 336
 337     switch (pdgMother){
 338     case 111:
 339       dNdy = fYieldArray[kPizero];
 340       break;
 341     case 221:
 342       dNdy = fYieldArray[kEta];
 343       break;
 344     case 113:
 345       dNdy = fYieldArray[kRho];
 346       break;
 347     case 223:
 348       dNdy = fYieldArray[kOmega];
 349       break;
 350     case 331:
 351       dNdy = fYieldArray[kEtaprime];
 352       break;
 353     case 333:
 354       dNdy = fYieldArray[kPhi];
 355       break;
 356     case 443:
 357       dNdy = fYieldArray[kJpsi];
 358       break;
 359     case 22:
 360       dNdy = fYieldArray[kDirectRealGamma];
 361       break;
 362     case 220000:
 363       dNdy = fYieldArray[kDirectVirtGamma];
 364       break;
 365
 366     default:
 367       dNdy = 0.;
 368     }
 369
 370     weight = dNdy*part->GetWeight();
 371     part->SetWeight(weight);
 372   }
 373
 374   fHeader->SetNProduced(maxPart);
 375
 376
 377   TArrayF eventVertex;
 378   eventVertex.Set(3);
 379   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
 380
 381   fHeader->SetPrimaryVertex(eventVertex);
 382
 383   gAlice->SetGenEventHeader(fHeader);
 384 }