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 "AliGenCocktailEventHeader.h"
  35
  36 #include "AliGenCocktailEntry.h"
  37 #include "AliGenEMCocktail.h"
  38 #include "AliGenEMlib.h"
  39 #include "AliGenBox.h"
  40 #include "AliGenParam.h"
  41 #include "AliMC.h"
  42 #include "AliRun.h"
  43 #include "AliStack.h"
  44 #include "AliDecayer.h"
  45 #include "AliDecayerPythia.h"
  46 #include "AliLog.h"
  47 #include "AliGenCorrHF.h"
  48
  49 ClassImp(AliGenEMCocktail)
  50
  51 //________________________________________________________________________
  52 AliGenEMCocktail::AliGenEMCocktail()
  53 :AliGenCocktail(),
  54    fDecayer(0),
  55    fDecayMode(kAll),
  56    fWeightingMode(kNonAnalog),
  57    fNPart(1000),
  58   fYieldArray(),
  59   fPtSelect(0),
  60   fCentrality(0),
  61   fV2Systematic(0),
  62   fForceConv(kFALSE),
  63   fHeaviestParticle(kGENs)
  64 {
  65   // Constructor
  66
  67 }
  68
  69 //_________________________________________________________________________
  70 AliGenEMCocktail::~AliGenEMCocktail()
  71 {
  72   // Destructor
  73
  74 }
  75
  76 //_________________________________________________________________________
  77 void AliGenEMCocktail::CreateCocktail()
  78 {
  79   // create and add sources to the cocktail
  80
  81   fDecayer->SetForceDecay(fDecayMode);
  82   fDecayer->ForceDecay();
  83
  84   // Set kinematic limits
  85   Double_t ptMin  = fPtMin;
  86   Double_t ptMax  = fPtMax;
  87   Double_t yMin   = fYMin;;
  88   Double_t yMax   = fYMax;;
  89   Double_t phiMin = fPhiMin*180./TMath::Pi();
  90   Double_t phiMax = fPhiMax*180./TMath::Pi();
  91   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));
  92   AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
  93
  94   //Initialize user selection for Pt Parameterization and centrality:
  95   AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
  96
  97   // Create and add electron sources to the generator
  98
  99   // pizero
 100   AliGenParam *genpizero=0;
 101   Char_t namePizero[10];
 102   snprintf(namePizero,10,"Pizero");
 103   genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
 104   AddSource2Generator(namePizero,genpizero);
 105   TF1 *fPtPizero = genpizero->GetPt();
 106 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 107   fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
 108 #else
 109   fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 110 #endif
 111
 112   // eta
 113   if(fHeaviestParticle<kGenEta)return;
 114   AliGenParam *geneta=0;
 115   Char_t nameEta[10];
 116   snprintf(nameEta,10,"Eta");
 117   geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
 118   AddSource2Generator(nameEta,geneta);
 119   TF1 *fPtEta = geneta->GetPt();
 120 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 121   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
 122 #else
 123   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 124 #endif
 125
 126   // rho
 127   if(fHeaviestParticle<kGenRho)return;
 128   AliGenParam *genrho=0;
 129   Char_t nameRho[10];
 130   snprintf(nameRho,10,"Rho");
 131   genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
 132   AddSource2Generator(nameRho,genrho);
 133   TF1 *fPtRho = genrho->GetPt();
 134 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 135   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
 136 #else
 137   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 138 #endif
 139
 140   // omega
 141   if(fHeaviestParticle<kGenOmega)return;
 142   AliGenParam *genomega=0;
 143   Char_t nameOmega[10];
 144   snprintf(nameOmega,10,"Omega");
 145   genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
 146   AddSource2Generator(nameOmega,genomega);
 147   TF1 *fPtOmega = genomega->GetPt();
 148 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 149   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
 150 #else
 151   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 152 #endif
 153
 154   // etaprime
 155   if(fHeaviestParticle<kGenEtaprime)return;
 156   AliGenParam *genetaprime=0;
 157   Char_t nameEtaprime[10];
 158   snprintf(nameEtaprime,10,"Etaprime");
 159   genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
 160   AddSource2Generator(nameEtaprime,genetaprime);
 161   TF1 *fPtEtaprime = genetaprime->GetPt();
 162 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 163   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
 164 #else
 165   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 166 #endif
 167
 168   // phi
 169   if(fHeaviestParticle<kGenPhi)return;
 170   AliGenParam *genphi=0;
 171   Char_t namePhi[10];
 172   snprintf(namePhi,10,"Phi");
 173   genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
 174   AddSource2Generator(namePhi,genphi);
 175   TF1 *fPtPhi = genphi->GetPt();
 176 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 177   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
 178 #else
 179   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 180 #endif
 181
 182   // jpsi
 183   if(fHeaviestParticle<kGenJpsi)return;
 184   AliGenParam *genjpsi=0;
 185   Char_t nameJpsi[10];
 186   snprintf(nameJpsi,10,"Jpsi");
 187   genjpsi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
 188   AddSource2Generator(nameJpsi,genjpsi);
 189   TF1 *fPtJpsi = genjpsi->GetPt();
 190 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 191   fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
 192 #else
 193   fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 194 #endif
 195
 196 }
 197
 198 //-------------------------------------------------------------------
 199 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
 200                                          AliGenParam* const genSource)
 201 {
 202 // add sources to the cocktail
 203   Double_t phiMin = fPhiMin*180./TMath::Pi();
 204   Double_t phiMax = fPhiMax*180./TMath::Pi();
 205
 206   genSource->SetPtRange(fPtMin, fPtMax);
 207   genSource->SetYRange(fYMin, fYMax);
 208   genSource->SetPhiRange(phiMin, phiMax);
 209   genSource->SetWeighting(fWeightingMode);
 210   genSource->SetForceGammaConversion(fForceConv);
 211   if (!gMC) genSource->SetDecayer(fDecayer);
 212   genSource->Init();
 213
 214   AddGenerator(genSource,nameSource,1.); // Adding Generator
 215 }
 216
 217 //-------------------------------------------------------------------
 218 void AliGenEMCocktail::Init()
 219 {
 220   // Initialisation
 221   TIter next(fEntries);
 222   AliGenCocktailEntry *entry;
 223   if (fStack) {
 224     while((entry = (AliGenCocktailEntry*)next())) {
 225       entry->Generator()->SetStack(fStack);
 226     }
 227   }
 228 }
 229
 230 //_________________________________________________________________________
 231 void AliGenEMCocktail::Generate()
 232 {
 233   // Generate event
 234   TIter next(fEntries);
 235   AliGenCocktailEntry *entry = 0;
 236   AliGenCocktailEntry *preventry = 0;
 237   AliGenerator* gen = 0;
 238
 239   if (fHeader) delete fHeader;
 240   fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
 241
 242   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 243
 244   // Generate the vertex position used by all generators
 245   if(fVertexSmear == kPerEvent) Vertex();
 246
 247   //Reseting stack
 248   AliRunLoader * runloader = AliRunLoader::Instance();
 249   if (runloader)
 250     if (runloader->Stack())
 251       runloader->Stack()->Clean();
 252
 253   // Loop over generators and generate events
 254   Int_t igen = 0;
 255   Float_t evPlane;
 256   Rndm(&evPlane,1);
 257   evPlane*=TMath::Pi()*2;
 258   while((entry = (AliGenCocktailEntry*)next())) {
 259     gen = entry->Generator();
 260     gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 261
 262     if (fNPart > 0) {
 263       igen++;
 264       if (igen == 1) entry->SetFirst(0);
 265       else  entry->SetFirst((partArray->GetEntriesFast())+1);
 266       gen->SetEventPlane(evPlane);
 267       gen->SetNumberParticles(fNPart);
 268       gen->Generate();
 269       entry->SetLast(partArray->GetEntriesFast());
 270       preventry = entry;
 271     }
 272   }
 273   next.Reset();
 274
 275   // Setting weights for proper absolute normalization
 276   Int_t iPart, iMother;
 277   Int_t pdgMother = 0;
 278   Double_t weight = 0.;
 279   Double_t dNdy = 0.;
 280   Int_t maxPart = partArray->GetEntriesFast();
 281   for(iPart=0; iPart<maxPart; iPart++){
 282     TParticle *part = gAlice->GetMCApp()->Particle(iPart);
 283     iMother = part->GetFirstMother();
 284     TParticle *mother = 0;
 285     if (iMother>=0){
 286       mother = gAlice->GetMCApp()->Particle(iMother);
 287       pdgMother = mother->GetPdgCode();
 288     }
 289     else
 290       pdgMother = part->GetPdgCode();
 291     switch (pdgMother){
 292     case 111:
 293       dNdy = fYieldArray[kGenPizero];
 294       break;
 295     case 221:
 296       dNdy = fYieldArray[kGenEta];
 297       break;
 298     case 113:
 299       dNdy = fYieldArray[kGenRho];
 300       break;
 301     case 223:
 302       dNdy = fYieldArray[kGenOmega];
 303       break;
 304     case 331:
 305       dNdy = fYieldArray[kGenEtaprime];
 306       break;
 307     case 333:
 308       dNdy = fYieldArray[kGenPhi];
 309       break;
 310     case 443:
 311       dNdy = fYieldArray[kGenJpsi];
 312       break;
 313
 314     default:
 315       dNdy = 0.;
 316     }
 317     weight = dNdy*part->GetWeight();
 318     part->SetWeight(weight);
 319   }
 320
 321   fHeader->SetNProduced(maxPart);
 322
 323
 324   TArrayF eventVertex;
 325   eventVertex.Set(3);
 326   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
 327
 328   fHeader->SetPrimaryVertex(eventVertex);
 329
 330   gAlice->SetGenEventHeader(fHeader);
 331 }