]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliGenEMCocktail.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
CMake: Removing all lhapdf link dependencies
[u/mrichter/AliRoot.git] / EVGEN / AliGenEMCocktail.cxx
CommitLineData
e40b9538 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
4ae1c9f0 18// Class to create the cocktail for physics with electrons, di-electrons,
e40b9538 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
4ae1c9f0 28
e40b9538 29#include <TObjArray.h>
30#include <TParticle.h>
31#include <TF1.h>
32#include <TVirtualMC.h>
33#include <TPDGCode.h>
71443190 34#include <TDatabasePDG.h>
e40b9538 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"
e40b9538 45#include "AliLog.h"
46#include "AliGenCorrHF.h"
47
48ClassImp(AliGenEMCocktail)
49
50//________________________________________________________________________
51AliGenEMCocktail::AliGenEMCocktail()
4ae1c9f0 52:AliGenCocktail(),
41b25ae9 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)
e40b9538 67{
41b25ae9 68 // Constructor
e40b9538 69}
70
71//_________________________________________________________________________
72AliGenEMCocktail::~AliGenEMCocktail()
73{
41b25ae9 74 // Destructor
e40b9538 75}
76
1b589c8a 77//_________________________________________________________________________
41b25ae9 78void AliGenEMCocktail::SetHeaviestHadron(ParticleGenerator_t part)
1b589c8a 79{
41b25ae9 80 Int_t val=kGenPizero;
81 while(val<part) val|=val<<1;
1b589c8a 82
41b25ae9 83 fSelectedParticles=val;
84 return;
1b589c8a 85}
86
e40b9538 87//_________________________________________________________________________
88void AliGenEMCocktail::CreateCocktail()
89{
41b25ae9 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
9f9ba3fd 117 // genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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 }
e40b9538 419}
420
421//-------------------------------------------------------------------
422void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
423 AliGenParam* const genSource)
424{
41b25ae9 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
e40b9538 437}
438
439//-------------------------------------------------------------------
440void AliGenEMCocktail::Init()
441{
41b25ae9 442 // Initialisation
443 TIter next(fEntries);
444 AliGenCocktailEntry *entry;
445 if (fStack) {
446 while((entry = (AliGenCocktailEntry*)next())) {
447 entry->Generator()->SetStack(fStack);
448 }
449 }
e40b9538 450}
451
452//_________________________________________________________________________
453void AliGenEMCocktail::Generate()
454{
41b25ae9 455 // Generate event
456 TIter next(fEntries);
457 AliGenCocktailEntry *entry = 0;
41b25ae9 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());
41b25ae9 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);
e40b9538 582}