SetSeed implementation
[u/mrichter/AliRoot.git] / 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 }