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