1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliGenEMCocktail.cxx 40702 2010-04-26 13:09:52Z morsch $ */
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)
29 #include <TObjArray.h>
30 #include <TParticle.h>
32 #include <TVirtualMC.h>
34 #include <TDatabasePDG.h>
35 #include "AliGenCocktailEventHeader.h"
37 #include "AliGenCocktailEntry.h"
38 #include "AliGenEMCocktail.h"
39 #include "AliGenEMlib.h"
40 #include "AliGenBox.h"
41 #include "AliGenParam.h"
45 #include "AliDecayer.h"
46 #include "AliDecayerPythia.h"
48 #include "AliGenCorrHF.h"
50 ClassImp(AliGenEMCocktail)
52 //________________________________________________________________________
53 AliGenEMCocktail::AliGenEMCocktail()
57 fWeightingMode(kNonAnalog),
64 fSelectedParticles(kGenHadrons)
70 //_________________________________________________________________________
71 AliGenEMCocktail::~AliGenEMCocktail()
77 //_________________________________________________________________________
78 void AliGenEMCocktail::CreateCocktail()
80 // create and add sources to the cocktail
82 fDecayer->SetForceDecay(fDecayMode);
83 fDecayer->ForceDecay();
85 // Set kinematic limits
86 Double_t ptMin = fPtMin;
87 Double_t ptMax = fPtMax;
88 Double_t yMin = fYMin;;
89 Double_t yMax = fYMax;;
90 Double_t phiMin = fPhiMin*180./TMath::Pi();
91 Double_t phiMax = fPhiMax*180./TMath::Pi();
92 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));
93 AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
95 //Initialize user selection for Pt Parameterization and centrality:
96 AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
98 // Create and add electron sources to the generator
100 if(fSelectedParticles&kGenPizero){
101 AliGenParam *genpizero=0;
102 Char_t namePizero[10];
103 snprintf(namePizero,10,"Pizero");
104 genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
105 genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
106 AddSource2Generator(namePizero,genpizero);
107 TF1 *fPtPizero = genpizero->GetPt();
108 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
109 fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
111 fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
116 if(fSelectedParticles&kGenEta){
117 AliGenParam *geneta=0;
119 snprintf(nameEta,10,"Eta");
120 geneta = new AliGenParam(fNPart/0.825, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
121 geneta->SetYRange(fYMin/0.825, fYMax/0.825);
122 AddSource2Generator(nameEta,geneta);
123 TF1 *fPtEta = geneta->GetPt();
124 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
125 fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
127 fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
132 if(fSelectedParticles&kGenRho){
133 AliGenParam *genrho=0;
135 snprintf(nameRho,10,"Rho");
136 genrho = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
137 genrho->SetYRange(fYMin/0.775, fYMax/0.775);
138 AddSource2Generator(nameRho,genrho);
139 TF1 *fPtRho = genrho->GetPt();
140 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
141 fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
143 fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
148 if(fSelectedParticles&kGenOmega){
149 AliGenParam *genomega=0;
150 Char_t nameOmega[10];
151 snprintf(nameOmega,10,"Omega");
152 genomega = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
153 genomega->SetYRange(fYMin/0.775, fYMax/0.775);
154 AddSource2Generator(nameOmega,genomega);
155 TF1 *fPtOmega = genomega->GetPt();
156 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
157 fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
159 fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
164 if(fSelectedParticles&kGenEtaprime){
165 AliGenParam *genetaprime=0;
166 Char_t nameEtaprime[10];
167 snprintf(nameEtaprime,10,"Etaprime");
168 genetaprime = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
169 genetaprime->SetYRange(fYMin/0.725, fYMax/0.725);
170 AddSource2Generator(nameEtaprime,genetaprime);
171 TF1 *fPtEtaprime = genetaprime->GetPt();
172 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
173 fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
175 fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
180 if(fSelectedParticles&kGenPhi){
181 AliGenParam *genphi=0;
183 snprintf(namePhi,10,"Phi");
184 genphi = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
185 genphi->SetYRange(fYMin/0.725, fYMax/0.725);
186 AddSource2Generator(namePhi,genphi);
187 TF1 *fPtPhi = genphi->GetPt();
188 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
189 fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
191 fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
196 if(fSelectedParticles&kGenJpsi){
197 AliGenParam *genjpsi=0;
199 snprintf(nameJpsi,10,"Jpsi");
200 genjpsi = new AliGenParam(fNPart/0.525, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
201 genjpsi->SetYRange(fYMin/0.525, fYMax/0.525);
202 AddSource2Generator(nameJpsi,genjpsi);
203 TF1 *fPtJpsi = genjpsi->GetPt();
204 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
205 fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
207 fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
212 if(fDecayMode!=kGammaEM) return;
213 if(fSelectedParticles&kGenDirectRealGamma){
214 AliGenParam *genDirectRealG=0;
215 Char_t nameDirectRealG[10];
216 snprintf(nameDirectRealG,10,"DirectRealGamma");
217 genDirectRealG = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDirectRealGamma, "DUMMY");
218 genDirectRealG->SetYRange(fYMin, fYMax);
219 AddSource2Generator(nameDirectRealG,genDirectRealG);
220 TF1 *fPtDirectRealG = genDirectRealG->GetPt();
221 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
222 fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
224 fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
228 if(fSelectedParticles&kGenDirectVirtGamma){
229 TDatabasePDG::Instance()->AddParticle("DirectVirtGamma","DirectVirtGamma",0,true,0,0,"GaugeBoson",220000);
230 AliGenParam *genDirectVirtG=0;
231 Char_t nameDirectVirtG[10];
232 snprintf(nameDirectVirtG,10,"DirectVirtGamma");
233 genDirectVirtG = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kDirectVirtGamma, "DUMMY");
234 genDirectVirtG->SetYRange(fYMin/0.775, fYMax/0.775);
235 AddSource2Generator(nameDirectVirtG,genDirectVirtG);
236 TF1 *fPtDirectVirtG = genDirectVirtG->GetPt();
237 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
238 fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,1.e-6);
240 fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
245 //-------------------------------------------------------------------
246 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
247 AliGenParam* const genSource)
249 // add sources to the cocktail
250 Double_t phiMin = fPhiMin*180./TMath::Pi();
251 Double_t phiMax = fPhiMax*180./TMath::Pi();
253 genSource->SetPtRange(fPtMin, fPtMax);
254 genSource->SetPhiRange(phiMin, phiMax);
255 genSource->SetWeighting(fWeightingMode);
256 genSource->SetForceGammaConversion(fForceConv);
257 if (!TVirtualMC::GetMC()) genSource->SetDecayer(fDecayer);
260 AddGenerator(genSource,nameSource,1.); // Adding Generator
263 //-------------------------------------------------------------------
264 void AliGenEMCocktail::Init()
267 TIter next(fEntries);
268 AliGenCocktailEntry *entry;
270 while((entry = (AliGenCocktailEntry*)next())) {
271 entry->Generator()->SetStack(fStack);
276 //_________________________________________________________________________
277 void AliGenEMCocktail::Generate()
280 TIter next(fEntries);
281 AliGenCocktailEntry *entry = 0;
282 AliGenCocktailEntry *preventry = 0;
283 AliGenerator* gen = 0;
285 if (fHeader) delete fHeader;
286 fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
288 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
290 // Generate the vertex position used by all generators
291 if(fVertexSmear == kPerEvent) Vertex();
294 AliRunLoader * runloader = AliRunLoader::Instance();
296 if (runloader->Stack())
297 runloader->Stack()->Clean();
299 // Loop over generators and generate events
303 evPlane*=TMath::Pi()*2;
304 while((entry = (AliGenCocktailEntry*)next())) {
305 gen = entry->Generator();
306 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
310 if (igen == 1) entry->SetFirst(0);
311 else entry->SetFirst((partArray->GetEntriesFast())+1);
312 gen->SetEventPlane(evPlane);
314 entry->SetLast(partArray->GetEntriesFast());
320 // Setting weights for proper absolute normalization
321 Int_t iPart, iMother;
323 Double_t weight = 0.;
325 Int_t maxPart = partArray->GetEntriesFast();
326 for(iPart=0; iPart<maxPart; iPart++){
327 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
328 iMother = part->GetFirstMother();
329 TParticle *mother = 0;
331 mother = gAlice->GetMCApp()->Particle(iMother);
332 pdgMother = mother->GetPdgCode();
335 pdgMother = part->GetPdgCode();
339 dNdy = fYieldArray[kPizero];
342 dNdy = fYieldArray[kEta];
345 dNdy = fYieldArray[kRho];
348 dNdy = fYieldArray[kOmega];
351 dNdy = fYieldArray[kEtaprime];
354 dNdy = fYieldArray[kPhi];
357 dNdy = fYieldArray[kJpsi];
360 dNdy = fYieldArray[kDirectRealGamma];
363 dNdy = fYieldArray[kDirectVirtGamma];
370 weight = dNdy*part->GetWeight();
371 part->SetWeight(weight);
374 fHeader->SetNProduced(maxPart);
379 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
381 fHeader->SetPrimaryVertex(eventVertex);
383 gAlice->SetGenEventHeader(fHeader);