Updated macros for L* analysis (Neelima)
[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
18// Class to create cocktails 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#include <TObjArray.h>
29#include <TParticle.h>
30#include <TF1.h>
31#include <TVirtualMC.h>
32#include <TPDGCode.h>
33#include "AliGenCocktailEventHeader.h"
34
35#include "AliGenCocktailEntry.h"
36#include "AliGenEMCocktail.h"
37#include "AliGenEMlib.h"
38#include "AliGenBox.h"
39#include "AliGenParam.h"
40#include "AliMC.h"
41#include "AliRun.h"
42#include "AliStack.h"
43#include "AliDecayer.h"
44#include "AliDecayerPythia.h"
45#include "AliLog.h"
46#include "AliGenCorrHF.h"
47
48ClassImp(AliGenEMCocktail)
49
50//________________________________________________________________________
51AliGenEMCocktail::AliGenEMCocktail()
52 :AliGenCocktail(),
53 fDecayer(0),
54 fDecayMode(kAll),
55 fWeightingMode(kNonAnalog),
56 fNPart(1000),
57 fYieldArray()
58{
59// Constructor
60
61}
62
63//_________________________________________________________________________
64AliGenEMCocktail::~AliGenEMCocktail()
65{
66// Destructor
67
68}
69
70//_________________________________________________________________________
71void AliGenEMCocktail::CreateCocktail()
72{
73// create and add sources to the cocktail
74
75 fDecayer->SetForceDecay(fDecayMode);
76 fDecayer->ForceDecay();
77
78// Set kinematic limits
79 Double_t ptMin = fPtMin;
80 Double_t ptMax = fPtMax;
81 Double_t yMin = fYMin;;
82 Double_t yMax = fYMax;;
83 Double_t phiMin = fPhiMin*180./TMath::Pi();
84 Double_t phiMax = fPhiMax*180./TMath::Pi();
85 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));
86 AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
87
88// Create and add electron sources to the generator
89
90// pizero
91 AliGenParam * genpizero=0;
92 Char_t namePizero[10];
2e13c33b 93 snprintf(namePizero,10, "Pizero");
e40b9538 94 genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
95 AddSource2Generator(namePizero,genpizero);
96 TF1 *fPtPizero = genpizero->GetPt();
faf00b0d 97#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
98 fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
99#else
100 fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
101#endif
e40b9538 102// eta
103 AliGenParam * geneta=0;
104 Char_t nameEta[10];
2e13c33b 105 snprintf(nameEta,10, "Eta");
e40b9538 106 geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
107 AddSource2Generator(nameEta,geneta);
108 TF1 *fPtEta = geneta->GetPt();
faf00b0d 109#if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
110 fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
111#else
112 fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
113#endif
e40b9538 114// rho
115 AliGenParam * genrho=0;
116 Char_t nameRho[10];
2e13c33b 117 snprintf(nameRho,10, "Rho");
e40b9538 118 genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
119 AddSource2Generator(nameRho,genrho);
120 TF1 *fPtRho = genrho->GetPt();
faf00b0d 121#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
122 fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
123#else
124 fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
125#endif
e40b9538 126// omega
127 AliGenParam * genomega=0;
128 Char_t nameOmega[10];
2e13c33b 129 snprintf(nameOmega,10, "Omega");
e40b9538 130 genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
131 AddSource2Generator(nameOmega,genomega);
132 TF1 *fPtOmega = genomega->GetPt();
faf00b0d 133#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
134 fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
135#else
136 fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(double*)0,1.e-6);
137#endif
e40b9538 138// etaprime
139 AliGenParam * genetaprime=0;
140 Char_t nameEtaprime[10];
2e13c33b 141 snprintf(nameEtaprime,10, "Etaprime");
e40b9538 142 genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
143 AddSource2Generator(nameEtaprime,genetaprime);
144 TF1 *fPtEtaprime = genetaprime->GetPt();
faf00b0d 145#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
146 fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
147#else
148 fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
149#endif
e40b9538 150// phi
151 AliGenParam * genphi=0;
152 Char_t namePhi[10];
2e13c33b 153 snprintf(namePhi, 10, "Phi");
e40b9538 154 genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
155 AddSource2Generator(namePhi,genphi);
156 TF1 *fPtPhi = genphi->GetPt();
faf00b0d 157#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
158 fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
159#else
160 fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t*)0,1.e-6);
161#endif
e40b9538 162}
163
164//-------------------------------------------------------------------
165void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
166 AliGenParam* const genSource)
167{
168// add sources to the cocktail
169 Double_t phiMin = fPhiMin*180./TMath::Pi();
170 Double_t phiMax = fPhiMax*180./TMath::Pi();
171
172 genSource->SetPtRange(fPtMin, fPtMax);
173 genSource->SetYRange(fYMin, fYMax);
174 genSource->SetPhiRange(phiMin, phiMax);
175 genSource->SetWeighting(fWeightingMode);
176 if (!gMC) genSource->SetDecayer(fDecayer);
177 genSource->Init();
178
179 AddGenerator(genSource,nameSource,1.); // Adding Generator
180}
181
182//-------------------------------------------------------------------
183void AliGenEMCocktail::Init()
184{
185// Initialisation
186 TIter next(fEntries);
187 AliGenCocktailEntry *entry;
188 if (fStack) {
189 while((entry = (AliGenCocktailEntry*)next())) {
190 entry->Generator()->SetStack(fStack);
191 }
192 }
193}
194
195//_________________________________________________________________________
196void AliGenEMCocktail::Generate()
197{
198// Generate event
199 TIter next(fEntries);
200 AliGenCocktailEntry *entry = 0;
201 AliGenCocktailEntry *preventry = 0;
202 AliGenerator* gen = 0;
203
204 if (fHeader) delete fHeader;
205 fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
206
207 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
208
209// Generate the vertex position used by all generators
210 if(fVertexSmear == kPerEvent) Vertex();
211
212//Reseting stack
213 AliRunLoader * runloader = AliRunLoader::Instance();
214 if (runloader)
215 if (runloader->Stack())
216 runloader->Stack()->Clean();
217
218// Loop over generators and generate events
219 Int_t igen = 0;
6078e216 220 Float_t evPlane;
221 Rndm(&evPlane,1);
222 evPlane*=TMath::Pi()*2;
e40b9538 223 while((entry = (AliGenCocktailEntry*)next())) {
224 gen = entry->Generator();
e40b9538 225 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
21391258 226 gen->SetTime(fTime);
e40b9538 227
228 if (fNPart > 0) {
229 igen++;
230 if (igen == 1) entry->SetFirst(0);
231 else entry->SetFirst((partArray->GetEntriesFast())+1);
6078e216 232 gen->SetEventPlane(evPlane);
e40b9538 233 gen->SetNumberParticles(fNPart);
234 gen->Generate();
235 entry->SetLast(partArray->GetEntriesFast());
236 preventry = entry;
237 }
238 }
239 next.Reset();
240
241// Setting weights for proper absolute normalization
242 Int_t iPart, iMother;
243 Int_t pdgMother = 0;
244 Double_t weight = 0.;
245 Double_t dNdy = 0.;
246 Int_t maxPart = partArray->GetEntriesFast();
247 for(iPart=0; iPart<maxPart; iPart++){
248 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
249 iMother = part->GetFirstMother();
250 TParticle *mother = 0;
251 if (iMother>=0){
252 mother = gAlice->GetMCApp()->Particle(iMother);
253 pdgMother = mother->GetPdgCode();
254 }
255 else
256 pdgMother = part->GetPdgCode();
257 switch (pdgMother){
258 case 111:
259 dNdy = fYieldArray[kGenPizero];
260 break;
261 case 221:
262 dNdy = fYieldArray[kGenEta];
263 break;
264 case 113:
265 dNdy = fYieldArray[kGenRho];
266 break;
267 case 223:
268 dNdy = fYieldArray[kGenOmega];
269 break;
270 case 331:
271 dNdy = fYieldArray[kGenEtaprime];
272 break;
273 case 333:
274 dNdy = fYieldArray[kGenPhi];
275 break;
276
277 default:
278 dNdy = 0.;
279 }
280
281 weight = dNdy*part->GetWeight();
282 part->SetWeight(weight);
283 }
284
285 fHeader->SetNProduced(maxPart);
286
287
288 TArrayF eventVertex;
289 eventVertex.Set(3);
290 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
291
292 fHeader->SetPrimaryVertex(eventVertex);
21391258 293 fHeader->SetInteractionTime(fTime);
e40b9538 294
295 gAlice->SetGenEventHeader(fHeader);
296}
297
298