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),
60 fPtSelect(AliGenEMlib::kPizero7TeVpp),
61 fCentrality(AliGenEMlib::kpp),
62 fV2Systematic(AliGenEMlib::kNoV2Sys),
64 fSelectedParticles(kGenPizero)
67 SetHeaviestHadron(kGenPhi);
70 //_________________________________________________________________________
71 AliGenEMCocktail::~AliGenEMCocktail()
77 //_________________________________________________________________________
78 void AliGenEMCocktail::SetHeaviestHadron(ParticeGenerator_t part)
81 while(val<part) val|=val<<1;
83 fSelectedParticles=val;
86 //_________________________________________________________________________
87 void AliGenEMCocktail::CreateCocktail()
89 // create and add sources to the cocktail
91 fDecayer->SetForceDecay(fDecayMode);
92 fDecayer->ForceDecay();
94 // Set kinematic limits
95 Double_t ptMin = fPtMin;
96 Double_t ptMax = fPtMax;
97 Double_t yMin = fYMin;;
98 Double_t yMax = fYMax;;
99 Double_t phiMin = fPhiMin*180./TMath::Pi();
100 Double_t phiMax = fPhiMax*180./TMath::Pi();
101 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));
102 AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
104 //Initialize user selection for Pt Parameterization and centrality:
105 AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
107 // Create and add electron sources to the generator
109 if(fSelectedParticles&kGenPizero){
110 AliGenParam *genpizero=0;
111 Char_t namePizero[10];
112 snprintf(namePizero,10,"Pizero");
113 //fNPart/0.925: increase number of particles so that we have the chosen number of particles in the chosen eta range
114 genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
115 //fYMin/0.925: increase eta range, so that the electron yield is constant (<5% change) over the chosen eta range
116 genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
117 AddSource2Generator(namePizero,genpizero);
118 TF1 *fPtPizero = genpizero->GetPt();
119 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
120 fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
122 fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
127 if(fSelectedParticles&kGenEta){
128 AliGenParam *geneta=0;
130 snprintf(nameEta,10,"Eta");
131 geneta = new AliGenParam(fNPart/0.825, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
132 geneta->SetYRange(fYMin/0.825, fYMax/0.825);
133 AddSource2Generator(nameEta,geneta);
134 TF1 *fPtEta = geneta->GetPt();
135 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
136 fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
138 fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
143 if(fSelectedParticles&kGenRho){
144 AliGenParam *genrho=0;
146 snprintf(nameRho,10,"Rho");
147 genrho = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
148 genrho->SetYRange(fYMin/0.775, fYMax/0.775);
149 AddSource2Generator(nameRho,genrho);
150 TF1 *fPtRho = genrho->GetPt();
151 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
152 fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
154 fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
159 if(fSelectedParticles&kGenOmega){
160 AliGenParam *genomega=0;
161 Char_t nameOmega[10];
162 snprintf(nameOmega,10,"Omega");
163 genomega = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
164 genomega->SetYRange(fYMin/0.775, fYMax/0.775);
165 AddSource2Generator(nameOmega,genomega);
166 TF1 *fPtOmega = genomega->GetPt();
167 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
168 fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
170 fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
175 if(fSelectedParticles&kGenEtaprime){
176 AliGenParam *genetaprime=0;
177 Char_t nameEtaprime[10];
178 snprintf(nameEtaprime,10,"Etaprime");
179 genetaprime = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
180 genetaprime->SetYRange(fYMin/0.725, fYMax/0.725);
181 AddSource2Generator(nameEtaprime,genetaprime);
182 TF1 *fPtEtaprime = genetaprime->GetPt();
183 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
184 fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
186 fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
191 if(fSelectedParticles&kGenPhi){
192 AliGenParam *genphi=0;
194 snprintf(namePhi,10,"Phi");
195 genphi = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
196 genphi->SetYRange(fYMin/0.725, fYMax/0.725);
197 AddSource2Generator(namePhi,genphi);
198 TF1 *fPtPhi = genphi->GetPt();
199 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
200 fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
202 fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
207 if(fSelectedParticles&kGenJpsi){
208 AliGenParam *genjpsi=0;
210 snprintf(nameJpsi,10,"Jpsi");
211 genjpsi = new AliGenParam(fNPart/0.525, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
212 genjpsi->SetYRange(fYMin/0.525, fYMax/0.525);
213 AddSource2Generator(nameJpsi,genjpsi);
214 TF1 *fPtJpsi = genjpsi->GetPt();
215 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
216 fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
218 fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
223 if(fDecayMode!=kGammaEM) return;
224 if(fSelectedParticles&kGenDirectRealGamma){
225 AliGenParam *genDirectRealG=0;
226 Char_t nameDirectRealG[10];
227 snprintf(nameDirectRealG,10,"DirectRealGamma");
228 genDirectRealG = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDirectRealGamma, "DUMMY");
229 genDirectRealG->SetYRange(fYMin, fYMax);
230 AddSource2Generator(nameDirectRealG,genDirectRealG);
231 TF1 *fPtDirectRealG = genDirectRealG->GetPt();
232 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
233 fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
235 fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
239 if(fSelectedParticles&kGenDirectVirtGamma){
240 TDatabasePDG::Instance()->AddParticle("DirectVirtGamma","DirectVirtGamma",0,true,0,0,"GaugeBoson",220000);
241 AliGenParam *genDirectVirtG=0;
242 Char_t nameDirectVirtG[10];
243 snprintf(nameDirectVirtG,10,"DirectVirtGamma");
244 genDirectVirtG = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kDirectVirtGamma, "DUMMY");
245 genDirectVirtG->SetYRange(fYMin/0.775, fYMax/0.775);
246 AddSource2Generator(nameDirectVirtG,genDirectVirtG);
247 TF1 *fPtDirectVirtG = genDirectVirtG->GetPt();
248 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
249 fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,1.e-6);
251 fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
256 //-------------------------------------------------------------------
257 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
258 AliGenParam* const genSource)
260 // add sources to the cocktail
261 Double_t phiMin = fPhiMin*180./TMath::Pi();
262 Double_t phiMax = fPhiMax*180./TMath::Pi();
264 genSource->SetPtRange(fPtMin, fPtMax);
265 genSource->SetPhiRange(phiMin, phiMax);
266 genSource->SetWeighting(fWeightingMode);
267 genSource->SetForceGammaConversion(fForceConv);
268 if (!TVirtualMC::GetMC()) genSource->SetDecayer(fDecayer);
271 AddGenerator(genSource,nameSource,1.); // Adding Generator
274 //-------------------------------------------------------------------
275 void AliGenEMCocktail::Init()
278 TIter next(fEntries);
279 AliGenCocktailEntry *entry;
281 while((entry = (AliGenCocktailEntry*)next())) {
282 entry->Generator()->SetStack(fStack);
287 //_________________________________________________________________________
288 void AliGenEMCocktail::Generate()
291 TIter next(fEntries);
292 AliGenCocktailEntry *entry = 0;
293 AliGenCocktailEntry *preventry = 0;
294 AliGenerator* gen = 0;
296 if (fHeader) delete fHeader;
297 fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
299 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
301 // Generate the vertex position used by all generators
302 if(fVertexSmear == kPerEvent) Vertex();
305 AliRunLoader * runloader = AliRunLoader::Instance();
307 if (runloader->Stack())
308 runloader->Stack()->Clean();
310 // Loop over generators and generate events
314 evPlane*=TMath::Pi()*2;
315 while((entry = (AliGenCocktailEntry*)next())) {
316 gen = entry->Generator();
317 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
321 if (igen == 1) entry->SetFirst(0);
322 else entry->SetFirst((partArray->GetEntriesFast())+1);
323 gen->SetEventPlane(evPlane);
325 entry->SetLast(partArray->GetEntriesFast());
331 // Setting weights for proper absolute normalization
332 Int_t iPart, iMother;
334 Double_t weight = 0.;
336 Int_t maxPart = partArray->GetEntriesFast();
337 for(iPart=0; iPart<maxPart; iPart++){
338 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
339 iMother = part->GetFirstMother();
340 TParticle *mother = 0;
342 mother = gAlice->GetMCApp()->Particle(iMother);
343 pdgMother = mother->GetPdgCode();
346 pdgMother = part->GetPdgCode();
350 dNdy = fYieldArray[kPizero];
353 dNdy = fYieldArray[kEta];
356 dNdy = fYieldArray[kRho];
359 dNdy = fYieldArray[kOmega];
362 dNdy = fYieldArray[kEtaprime];
365 dNdy = fYieldArray[kPhi];
368 dNdy = fYieldArray[kJpsi];
371 dNdy = fYieldArray[kDirectRealGamma];
374 dNdy = fYieldArray[kDirectVirtGamma];
381 weight = dNdy*part->GetWeight();
382 part->SetWeight(weight);
385 fHeader->SetNProduced(maxPart);
390 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
392 fHeader->SetPrimaryVertex(eventVertex);
394 gAlice->SetGenEventHeader(fHeader);