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         fCollisionSystem(AliGenEMlib::kpp7TeV),
  61         fPtSelectPi0(AliGenEMlib::kPizeroParam),
  62         fPtSelectEta(AliGenEMlib::kEtaParampp),
  63         fPtSelectOmega(AliGenEMlib::kOmegaParampp),
  64         fPtSelectPhi(AliGenEMlib::kPhiParampp),
  65         fCentrality(AliGenEMlib::kpp),
  66         fV2Systematic(AliGenEMlib::kNoV2Sys),
  67         fForceConv(kFALSE),
  68         fSelectedParticles(kGenHadrons)
  69 {
  70         // Constructor
  71 }
  72
  73 //_________________________________________________________________________
  74 AliGenEMCocktail::~AliGenEMCocktail()
  75 {
  76         // Destructor
  77 }
  78
  79 //_________________________________________________________________________
  80 void AliGenEMCocktail::SetHeaviestHadron(ParticleGenerator_t part)
  81 {
  82         Int_t val=kGenPizero;
  83         while(val<part) val|=val<<1;
  84
  85         fSelectedParticles=val;
  86         return;
  87 }
  88
  89 //_________________________________________________________________________
  90 void AliGenEMCocktail::CreateCocktail()
  91 {
  92         // create and add sources to the cocktail
  93
  94         fDecayer->SetForceDecay(fDecayMode);
  95         fDecayer->ForceDecay();
  96
  97         // Set kinematic limits
  98         Double_t ptMin  = fPtMin;
  99         Double_t ptMax  = fPtMax;
 100         Double_t yMin   = fYMin;;
 101         Double_t yMax   = fYMax;;
 102         Double_t phiMin = fPhiMin*180./TMath::Pi();
 103         Double_t phiMax = fPhiMax*180./TMath::Pi();
 104         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));
 105         AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
 106         AliInfo(Form("Selected Params:collision system - %d , centrality - %d, pi0 param - %d, eta param - %d, omega param - %d, phi param - %d",fCollisionSystem, fCentrality, fPtSelectPi0, fPtSelectEta, fPtSelectOmega, fPtSelectPhi));
 107         //Initialize user selection for Pt Parameterization and centrality:
 108         AliGenEMlib::SelectParams(fCollisionSystem, fPtSelectPi0, fPtSelectEta, fPtSelectOmega, fPtSelectPhi, fCentrality,fV2Systematic);
 109
 110         // Create and add electron sources to the generator
 111         // pizero
 112         if(fSelectedParticles&kGenPizero){
 113                 AliGenParam *genpizero=0;
 114                 Char_t namePizero[10];
 115                 snprintf(namePizero,10,"Pizero");
 116                         //fNPart/0.925: increase number of particles so that we have the chosen number of particles in the chosen eta range
 117         //      genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
 118                         //fYMin/0.925: increase eta range, so that the electron yield is constant (<5% change) over the chosen eta range
 119                 genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
 120
 121                 // NOTE Theo: fNPart/0.925: increase number of particles so that we have the chosen number of particles in the chosen eta range
 122                 // NOTE Theo: fYMin/0.925: increase eta range, so that the electron yield is constant (<5% change) over the chosen eta range
 123                 // NOTE Friederike: the additional factors here cannot be fixed numbers, if you need them
 124                 //                                      generate a setting which puts them for you but never do it hardcoded - electrons are not the only ones
 125                 //                                      using the cocktail
 126                 genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
 127                 genpizero->SetYRange(fYMin, fYMax);
 128
 129                 AddSource2Generator(namePizero,genpizero);
 130                 TF1 *fPtPizero = genpizero->GetPt();
 131                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 132                         fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
 133                 #else
 134                         fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 135                 #endif
 136         }
 137
 138         // eta
 139         if(fSelectedParticles&kGenEta){
 140                 AliGenParam *geneta=0;
 141                 Char_t nameEta[10];
 142                 snprintf(nameEta,10,"Eta");
 143                 // NOTE: the additional factors are set back to one as they are not the same for photons and electrons
 144                 geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
 145                 geneta->SetYRange(fYMin, fYMax);
 146
 147                 AddSource2Generator(nameEta,geneta);
 148                 TF1 *fPtEta = geneta->GetPt();
 149                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 150                         fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
 151                 #else
 152                         fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 153                 #endif
 154         }
 155
 156         // rho
 157         if(fSelectedParticles&kGenRho0){
 158                 AliGenParam *genrho=0;
 159                 Char_t nameRho[10];
 160                 snprintf(nameRho,10,"Rho");
 161                 // NOTE: the additional factors are set back to one as they are not the same for photons and electrons
 162                 genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho0, "DUMMY");
 163                 genrho->SetYRange(fYMin, fYMax);
 164                 AddSource2Generator(nameRho,genrho);
 165                 TF1 *fPtRho = genrho->GetPt();
 166                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 167                         fYieldArray[kRho0] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
 168                 #else
 169                         fYieldArray[kRho0] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 170                 #endif
 171         }
 172
 173         // omega
 174         if(fSelectedParticles&kGenOmega){
 175                 AliGenParam *genomega=0;
 176                 Char_t nameOmega[10];
 177                 snprintf(nameOmega,10,"Omega");
 178                 // NOTE: the additional factors are set back to one as they are not the same for photons and electrons
 179                 genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
 180                 genomega->SetYRange(fYMin, fYMax);
 181                 AddSource2Generator(nameOmega,genomega);
 182                 TF1 *fPtOmega = genomega->GetPt();
 183                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 184                         fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
 185                 #else
 186                         fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 187                 #endif
 188         }
 189
 190         // etaprime
 191         if(fSelectedParticles&kGenEtaprime){
 192                 AliGenParam *genetaprime=0;
 193                 Char_t nameEtaprime[10];
 194                 snprintf(nameEtaprime,10,"Etaprime");
 195                 // NOTE: the additional factors are set back to one as they are not the same for photons and electrons
 196                 genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
 197                 genetaprime->SetYRange(fYMin, fYMax);
 198                 AddSource2Generator(nameEtaprime,genetaprime);
 199                 TF1 *fPtEtaprime = genetaprime->GetPt();
 200                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 201                         fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
 202                 #else
 203                         fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 204                 #endif
 205         }
 206
 207         // phi
 208         if(fSelectedParticles&kGenPhi){
 209                 AliGenParam *genphi=0;
 210                 Char_t namePhi[10];
 211                 snprintf(namePhi,10,"Phi");
 212                 // NOTE: the additional factors are set back to one as they are not the same for photons and electrons
 213                 genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
 214                 genphi->SetYRange(fYMin, fYMax);
 215                 AddSource2Generator(namePhi,genphi);
 216                 TF1 *fPtPhi = genphi->GetPt();
 217                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 218                         fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
 219                 #else
 220                         fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 221                 #endif
 222         }
 223
 224         // jpsi
 225         if(fSelectedParticles&kGenJpsi){
 226                 AliGenParam *genjpsi=0;
 227                 Char_t nameJpsi[10];
 228                 snprintf(nameJpsi,10,"Jpsi");
 229                 // NOTE: the additional factors are set back to one as they are not the same for photons and electrons
 230                 genjpsi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
 231                 genjpsi->SetYRange(fYMin, fYMax);
 232                 AddSource2Generator(nameJpsi,genjpsi);
 233                 TF1 *fPtJpsi = genjpsi->GetPt();
 234                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 235                         fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
 236                 #else
 237                         fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 238                 #endif
 239         }
 240
 241         // sigma
 242         if(fSelectedParticles&kGenSigma0){
 243                 AliGenParam * gensigma=0;
 244                 Char_t nameSigma[10];
 245                 snprintf(nameSigma,10, "Sigma0");
 246                 gensigma = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kSigma0, "DUMMY");
 247                 gensigma->SetYRange(fYMin, fYMax);
 248
 249                 AddSource2Generator(nameSigma,gensigma);
 250                 TF1 *fPtSigma = gensigma->GetPt();
 251                 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
 252                         fYieldArray[kSigma0] = fPtSigma->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
 253                 #else
 254                         fYieldArray[kSigma0] = fPtSigma->Integral(fPtMin,fPtMax,1.e-6);
 255                 #endif
 256         }
 257
 258         // k0short
 259         if(fSelectedParticles&kGenK0s){
 260                 AliGenParam * genkzeroshort=0;
 261                 Char_t nameK0short[10];
 262                 snprintf(nameK0short, 10, "K0short");
 263                 genkzeroshort = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kK0s, "DUMMY");
 264                 genkzeroshort->SetYRange(fYMin, fYMax);
 265                 AddSource2Generator(nameK0short,genkzeroshort);
 266                 TF1 *fPtK0short = genkzeroshort->GetPt();
 267                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 268                         fYieldArray[kK0s] = fPtK0short->Integral(fPtMin,fPtMax,1.e-6);
 269                 #else
 270                         fYieldArray[kK0s] = fPtK0short->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
 271                 #endif
 272         }
 273
 274         // Delta++
 275         if(fSelectedParticles&kGenDeltaPlPl){
 276                 AliGenParam * genkdeltaPlPl=0;
 277                 Char_t nameDeltaPlPl[10];
 278                 snprintf(nameDeltaPlPl, 10, "DeltaPlPl");
 279                 genkdeltaPlPl = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDeltaPlPl, "DUMMY");
 280                 genkdeltaPlPl->SetYRange(fYMin, fYMax);
 281                 AddSource2Generator(nameDeltaPlPl,genkdeltaPlPl);
 282                 TF1 *fPtDeltaPlPl = genkdeltaPlPl->GetPt();
 283                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 284                         fYieldArray[kDeltaPlPl] = fPtDeltaPlPl->Integral(fPtMin,fPtMax,1.e-6);
 285                 #else
 286                         fYieldArray[kDeltaPlPl] = fPtDeltaPlPl->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
 287                 #endif
 288         }
 289
 290         // Delta+
 291         if(fSelectedParticles&kGenDeltaPl){
 292                 AliGenParam * genkdeltaPl=0;
 293                 Char_t nameDeltaPl[10];
 294                 snprintf(nameDeltaPl, 10, "DeltaPl");
 295                 genkdeltaPl = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDeltaPl, "DUMMY");
 296                 genkdeltaPl->SetYRange(fYMin, fYMax);
 297                 AddSource2Generator(nameDeltaPl,genkdeltaPl);
 298                 TF1 *fPtDeltaPl = genkdeltaPl->GetPt();
 299                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 300                         fYieldArray[kDeltaPl] = fPtDeltaPl->Integral(fPtMin,fPtMax,1.e-6);
 301                 #else
 302                         fYieldArray[kDeltaPl] = fPtDeltaPl->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
 303                 #endif
 304         }
 305
 306         // Delta-
 307         if(fSelectedParticles&kGenDeltaMi){
 308                 AliGenParam * genkdeltaMi=0;
 309                 Char_t nameDeltaMi[10];
 310                 snprintf(nameDeltaMi, 10, "DeltaMi");
 311                 genkdeltaMi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDeltaMi, "DUMMY");
 312                 genkdeltaMi->SetYRange(fYMin, fYMax);
 313                 AddSource2Generator(nameDeltaMi,genkdeltaMi);
 314                 TF1 *fPtDeltaMi = genkdeltaMi->GetPt();
 315                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 316                         fYieldArray[kDeltaMi] = fPtDeltaMi->Integral(fPtMin,fPtMax,1.e-6);
 317                 #else
 318                         fYieldArray[kDeltaMi] = fPtDeltaMi->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
 319                 #endif
 320         }
 321
 322         // Delta0
 323         if(fSelectedParticles&kGenDeltaZero){
 324                 AliGenParam * genkdeltaZero=0;
 325                 Char_t nameDeltaZero[10];
 326                 snprintf(nameDeltaZero, 10, "DeltaZero");
 327                 genkdeltaZero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDeltaZero, "DUMMY");
 328                 genkdeltaZero->SetYRange(fYMin, fYMax);
 329                 AddSource2Generator(nameDeltaZero,genkdeltaZero);
 330                 TF1 *fPtDeltaZero = genkdeltaZero->GetPt();
 331                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 332                         fYieldArray[kDeltaZero] = fPtDeltaZero->Integral(fPtMin,fPtMax,1.e-6);
 333                 #else
 334                         fYieldArray[kDeltaZero] = fPtDeltaZero->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
 335                 #endif
 336         }
 337
 338         // rho+
 339         if(fSelectedParticles&kGenRhoPl){
 340                 AliGenParam * genkrhoPl=0;
 341                 Char_t nameRhoPl[10];
 342                 snprintf(nameRhoPl, 10, "RhoPl");
 343                 genkrhoPl = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRhoPl, "DUMMY");
 344                 genkrhoPl->SetYRange(fYMin, fYMax);
 345                 AddSource2Generator(nameRhoPl,genkrhoPl);
 346                 TF1 *fPtRhoPl = genkrhoPl->GetPt();
 347                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 348                         fYieldArray[kRhoPl] = fPtRhoPl->Integral(fPtMin,fPtMax,1.e-6);
 349                 #else
 350                         fYieldArray[kRhoPl] = fPtRhoPl->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
 351                 #endif
 352         }
 353
 354         // rho-
 355         if(fSelectedParticles&kGenRhoMi){
 356                 AliGenParam * genkrhoMi=0;
 357                 Char_t nameRhoMi[10];
 358                 snprintf(nameRhoMi, 10, "RhoMi");
 359                 genkrhoMi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRhoMi, "DUMMY");
 360                 genkrhoMi->SetYRange(fYMin, fYMax);
 361                 AddSource2Generator(nameRhoMi,genkrhoMi);
 362                 TF1 *fPtRhoMi = genkrhoMi->GetPt();
 363                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 364                         fYieldArray[kRhoMi] = fPtRhoMi->Integral(fPtMin,fPtMax,1.e-6);
 365                 #else
 366                         fYieldArray[kRhoMi] = fPtRhoMi->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
 367                 #endif
 368         }
 369
 370         // K0*
 371         if(fSelectedParticles&kGenK0star){
 372                 AliGenParam * genkK0star=0;
 373                 Char_t nameK0star[10];
 374                 snprintf(nameK0star, 10, "K0star");
 375                 genkK0star = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kK0star, "DUMMY");
 376                 genkK0star->SetYRange(fYMin, fYMax);
 377                 AddSource2Generator(nameK0star,genkK0star);
 378                 TF1 *fPtK0star = genkK0star->GetPt();
 379                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 380                         fYieldArray[kK0star] = fPtK0star->Integral(fPtMin,fPtMax,1.e-6);
 381                 #else
 382                         fYieldArray[kK0star] = fPtK0star->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
 383                 #endif
 384         }
 385
 386         // direct gamma
 387         if(fDecayMode!=kGammaEM) return;
 388
 389         if(fSelectedParticles&kGenDirectRealGamma){
 390                 AliGenParam *genDirectRealG=0;
 391                 Char_t nameDirectRealG[10];
 392                 snprintf(nameDirectRealG,10,"DirectRealGamma");
 393                 // NOTE: the additional factors are set back to one as they are not the same for photons and electrons
 394                 genDirectRealG = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDirectRealGamma, "DUMMY");
 395                 genDirectRealG->SetYRange(fYMin, fYMax);
 396                 AddSource2Generator(nameDirectRealG,genDirectRealG);
 397                 TF1 *fPtDirectRealG = genDirectRealG->GetPt();
 398                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 399                         fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
 400                 #else
 401                         fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 402                 #endif
 403         }
 404
 405         if(fSelectedParticles&kGenDirectVirtGamma){
 406                 TDatabasePDG::Instance()->AddParticle("DirectVirtGamma","DirectVirtGamma",0,true,0,0,"GaugeBoson",220000);
 407                 AliGenParam *genDirectVirtG=0;
 408                 Char_t nameDirectVirtG[10];
 409                 snprintf(nameDirectVirtG,10,"DirectVirtGamma");
 410                 // NOTE: the additional factors are set back to one as they are not the same for photons and electrons
 411                 genDirectVirtG = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDirectVirtGamma, "DUMMY");
 412                 genDirectVirtG->SetYRange(fYMin, fYMax);
 413                 AddSource2Generator(nameDirectVirtG,genDirectVirtG);
 414                 TF1 *fPtDirectVirtG = genDirectVirtG->GetPt();
 415                 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
 416                         fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,1.e-6);
 417                 #else
 418                         fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
 419                 #endif
 420         }
 421 }
 422
 423 //-------------------------------------------------------------------
 424 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
 425                                          AliGenParam* const genSource)
 426 {
 427         // add sources to the cocktail
 428         Double_t phiMin = fPhiMin*180./TMath::Pi();
 429         Double_t phiMax = fPhiMax*180./TMath::Pi();
 430
 431         genSource->SetPtRange(fPtMin, fPtMax);
 432         genSource->SetPhiRange(phiMin, phiMax);
 433         genSource->SetWeighting(fWeightingMode);
 434         genSource->SetForceGammaConversion(fForceConv);
 435         if (!TVirtualMC::GetMC()) genSource->SetDecayer(fDecayer);
 436         genSource->Init();
 437
 438         AddGenerator(genSource,nameSource,1.); // Adding Generator
 439 }
 440
 441 //-------------------------------------------------------------------
 442 void AliGenEMCocktail::Init()
 443 {
 444         // Initialisation
 445         TIter next(fEntries);
 446         AliGenCocktailEntry *entry;
 447         if (fStack) {
 448                 while((entry = (AliGenCocktailEntry*)next())) {
 449                 entry->Generator()->SetStack(fStack);
 450                 }
 451         }
 452 }
 453
 454 //_________________________________________________________________________
 455 void AliGenEMCocktail::Generate()
 456 {
 457         // Generate event
 458         TIter next(fEntries);
 459         AliGenCocktailEntry *entry = 0;
 460         AliGenerator* gen = 0;
 461
 462         if (fHeader) delete fHeader;
 463         fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
 464
 465         const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 466
 467         // Generate the vertex position used by all generators
 468         if(fVertexSmear == kPerEvent) Vertex();
 469
 470         //Reseting stack
 471         AliRunLoader * runloader = AliRunLoader::Instance();
 472         if (runloader)
 473                 if (runloader->Stack())
 474                 runloader->Stack()->Clean();
 475
 476         // Loop over generators and generate events
 477         Int_t igen = 0;
 478         Float_t evPlane;
 479         Rndm(&evPlane,1);
 480         evPlane*=TMath::Pi()*2;
 481         while((entry = (AliGenCocktailEntry*)next())) {
 482                 gen = entry->Generator();
 483                 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 484
 485                 if (fNPart > 0) {
 486                 igen++;
 487                 if (igen == 1) entry->SetFirst(0);
 488                 else  entry->SetFirst((partArray->GetEntriesFast())+1);
 489                 gen->SetEventPlane(evPlane);
 490                 gen->Generate();
 491                 entry->SetLast(partArray->GetEntriesFast());
 492                 }
 493         }
 494         next.Reset();
 495
 496         // Setting weights for proper absolute normalization
 497         Int_t iPart, iMother;
 498         Int_t pdgMother = 0;
 499         Double_t weight = 0.;
 500         Double_t dNdy = 0.;
 501         Int_t maxPart = partArray->GetEntriesFast();
 502         for(iPart=0; iPart<maxPart; iPart++){
 503                 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
 504                 iMother = part->GetFirstMother();
 505                 TParticle *mother = 0;
 506                 if (iMother>=0){
 507                         mother = gAlice->GetMCApp()->Particle(iMother);
 508                         pdgMother = mother->GetPdgCode();
 509                 } else pdgMother = part->GetPdgCode();
 510
 511                 switch (pdgMother){
 512                         case 111:
 513                                 dNdy = fYieldArray[kPizero];
 514                                 break;
 515                         case 221:
 516                                 dNdy = fYieldArray[kEta];
 517                                 break;
 518                         case 113:
 519                                 dNdy = fYieldArray[kRho0];
 520                                 break;
 521                         case 223:
 522                                 dNdy = fYieldArray[kOmega];
 523                                 break;
 524                         case 331:
 525                                 dNdy = fYieldArray[kEtaprime];
 526                                 break;
 527                         case 333:
 528                                 dNdy = fYieldArray[kPhi];
 529                                 break;
 530                         case 443:
 531                                 dNdy = fYieldArray[kJpsi];
 532                                 break;
 533                         case 22:
 534                                 dNdy = fYieldArray[kDirectRealGamma];
 535                                 break;
 536                         case 220000:
 537                                 dNdy = fYieldArray[kDirectVirtGamma];
 538                                 break;
 539                         case 3212:
 540                                 dNdy = fYieldArray[kSigma0];
 541                                 break;
 542                         case 310:
 543                                 dNdy = fYieldArray[kK0s];
 544                                 break;
 545                         case 2224:
 546                                 dNdy = fYieldArray[kDeltaPlPl];
 547                                 break;
 548                         case 2214:
 549                                 dNdy = fYieldArray[kDeltaPl];
 550                                 break;
 551                         case 1114:
 552                                 dNdy = fYieldArray[kDeltaMi];
 553                                 break;
 554                         case 2114:
 555                                 dNdy = fYieldArray[kDeltaZero];
 556                                 break;
 557                         case 213:
 558                                 dNdy = fYieldArray[kRhoPl];
 559                                 break;
 560                         case -213:
 561                                 dNdy = fYieldArray[kRhoMi];
 562                                 break;
 563                         case 313:
 564                                 dNdy = fYieldArray[kK0star];
 565                                 break;
 566                         default:
 567                                 dNdy = 0.;
 568                 }
 569
 570                 weight = dNdy*part->GetWeight();
 571                 part->SetWeight(weight);
 572         }
 573
 574         fHeader->SetNProduced(maxPart);
 575
 576
 577         TArrayF eventVertex;
 578         eventVertex.Set(3);
 579         for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
 580
 581         fHeader->SetPrimaryVertex(eventVertex);
 582
 583         gAlice->SetGenEventHeader(fHeader);
 584 }