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