]> 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
Afterburner to generate light nuclei
[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"
45#include "AliDecayer.h"
46#include "AliDecayerPythia.h"
47#include "AliLog.h"
48#include "AliGenCorrHF.h"
49
50ClassImp(AliGenEMCocktail)
51
52//________________________________________________________________________
53AliGenEMCocktail::AliGenEMCocktail()
4ae1c9f0 54:AliGenCocktail(),
41b25ae9 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)
e40b9538 69{
41b25ae9 70 // Constructor
e40b9538 71}
72
73//_________________________________________________________________________
74AliGenEMCocktail::~AliGenEMCocktail()
75{
41b25ae9 76 // Destructor
e40b9538 77}
78
1b589c8a 79//_________________________________________________________________________
41b25ae9 80void AliGenEMCocktail::SetHeaviestHadron(ParticleGenerator_t part)
1b589c8a 81{
41b25ae9 82 Int_t val=kGenPizero;
83 while(val<part) val|=val<<1;
1b589c8a 84
41b25ae9 85 fSelectedParticles=val;
86 return;
1b589c8a 87}
88
e40b9538 89//_________________________________________________________________________
90void AliGenEMCocktail::CreateCocktail()
91{
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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();
41b25ae9 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 }
e40b9538 421}
422
423//-------------------------------------------------------------------
424void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
425 AliGenParam* const genSource)
426{
41b25ae9 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
e40b9538 439}
440
441//-------------------------------------------------------------------
442void AliGenEMCocktail::Init()
443{
41b25ae9 444 // Initialisation
445 TIter next(fEntries);
446 AliGenCocktailEntry *entry;
447 if (fStack) {
448 while((entry = (AliGenCocktailEntry*)next())) {
449 entry->Generator()->SetStack(fStack);
450 }
451 }
e40b9538 452}
453
454//_________________________________________________________________________
455void AliGenEMCocktail::Generate()
456{
41b25ae9 457 // Generate event
458 TIter next(fEntries);
459 AliGenCocktailEntry *entry = 0;
41b25ae9 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());
41b25ae9 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);
e40b9538 584}