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