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 "AliGenCocktailEventHeader.h"
36 #include "AliGenCocktailEntry.h"
37 #include "AliGenEMCocktail.h"
38 #include "AliGenEMlib.h"
39 #include "AliGenBox.h"
40 #include "AliGenParam.h"
44 #include "AliDecayer.h"
45 #include "AliDecayerPythia.h"
47 #include "AliGenCorrHF.h"
49 ClassImp(AliGenEMCocktail)
51 //________________________________________________________________________
52 AliGenEMCocktail::AliGenEMCocktail()
56 fWeightingMode(kNonAnalog),
63 fHeaviestParticle(kGENs)
69 //_________________________________________________________________________
70 AliGenEMCocktail::~AliGenEMCocktail()
76 //_________________________________________________________________________
77 void AliGenEMCocktail::CreateCocktail()
79 // create and add sources to the cocktail
81 fDecayer->SetForceDecay(fDecayMode);
82 fDecayer->ForceDecay();
84 // Set kinematic limits
85 Double_t ptMin = fPtMin;
86 Double_t ptMax = fPtMax;
87 Double_t yMin = fYMin;;
88 Double_t yMax = fYMax;;
89 Double_t phiMin = fPhiMin*180./TMath::Pi();
90 Double_t phiMax = fPhiMax*180./TMath::Pi();
91 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));
92 AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
94 //Initialize user selection for Pt Parameterization and centrality:
95 AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
97 // Create and add electron sources to the generator
100 AliGenParam *genpizero=0;
101 Char_t namePizero[10];
102 snprintf(namePizero,10,"Pizero");
103 genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
104 AddSource2Generator(namePizero,genpizero);
105 TF1 *fPtPizero = genpizero->GetPt();
106 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
107 fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
109 fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
113 if(fHeaviestParticle<kGenEta)return;
114 AliGenParam *geneta=0;
116 snprintf(nameEta,10,"Eta");
117 geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
118 AddSource2Generator(nameEta,geneta);
119 TF1 *fPtEta = geneta->GetPt();
120 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
121 fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
123 fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
127 if(fHeaviestParticle<kGenRho)return;
128 AliGenParam *genrho=0;
130 snprintf(nameRho,10,"Rho");
131 genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
132 AddSource2Generator(nameRho,genrho);
133 TF1 *fPtRho = genrho->GetPt();
134 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
135 fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
137 fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
141 if(fHeaviestParticle<kGenOmega)return;
142 AliGenParam *genomega=0;
143 Char_t nameOmega[10];
144 snprintf(nameOmega,10,"Omega");
145 genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
146 AddSource2Generator(nameOmega,genomega);
147 TF1 *fPtOmega = genomega->GetPt();
148 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
149 fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
151 fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
155 if(fHeaviestParticle<kGenEtaprime)return;
156 AliGenParam *genetaprime=0;
157 Char_t nameEtaprime[10];
158 snprintf(nameEtaprime,10,"Etaprime");
159 genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
160 AddSource2Generator(nameEtaprime,genetaprime);
161 TF1 *fPtEtaprime = genetaprime->GetPt();
162 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
163 fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
165 fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
169 if(fHeaviestParticle<kGenPhi)return;
170 AliGenParam *genphi=0;
172 snprintf(namePhi,10,"Phi");
173 genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
174 AddSource2Generator(namePhi,genphi);
175 TF1 *fPtPhi = genphi->GetPt();
176 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
177 fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
179 fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
183 if(fHeaviestParticle<kGenJpsi)return;
184 AliGenParam *genjpsi=0;
186 snprintf(nameJpsi,10,"Jpsi");
187 genjpsi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
188 AddSource2Generator(nameJpsi,genjpsi);
189 TF1 *fPtJpsi = genjpsi->GetPt();
190 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
191 fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
193 fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
198 //-------------------------------------------------------------------
199 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource,
200 AliGenParam* const genSource)
202 // add sources to the cocktail
203 Double_t phiMin = fPhiMin*180./TMath::Pi();
204 Double_t phiMax = fPhiMax*180./TMath::Pi();
206 genSource->SetPtRange(fPtMin, fPtMax);
207 genSource->SetYRange(fYMin, fYMax);
208 genSource->SetPhiRange(phiMin, phiMax);
209 genSource->SetWeighting(fWeightingMode);
210 genSource->SetForceGammaConversion(fForceConv);
211 if (!gMC) genSource->SetDecayer(fDecayer);
214 AddGenerator(genSource,nameSource,1.); // Adding Generator
217 //-------------------------------------------------------------------
218 void AliGenEMCocktail::Init()
221 TIter next(fEntries);
222 AliGenCocktailEntry *entry;
224 while((entry = (AliGenCocktailEntry*)next())) {
225 entry->Generator()->SetStack(fStack);
230 //_________________________________________________________________________
231 void AliGenEMCocktail::Generate()
234 TIter next(fEntries);
235 AliGenCocktailEntry *entry = 0;
236 AliGenCocktailEntry *preventry = 0;
237 AliGenerator* gen = 0;
239 if (fHeader) delete fHeader;
240 fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
242 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
244 // Generate the vertex position used by all generators
245 if(fVertexSmear == kPerEvent) Vertex();
248 AliRunLoader * runloader = AliRunLoader::Instance();
250 if (runloader->Stack())
251 runloader->Stack()->Clean();
253 // Loop over generators and generate events
257 evPlane*=TMath::Pi()*2;
258 while((entry = (AliGenCocktailEntry*)next())) {
259 gen = entry->Generator();
260 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
264 if (igen == 1) entry->SetFirst(0);
265 else entry->SetFirst((partArray->GetEntriesFast())+1);
266 gen->SetEventPlane(evPlane);
267 gen->SetNumberParticles(fNPart);
269 entry->SetLast(partArray->GetEntriesFast());
275 // Setting weights for proper absolute normalization
276 Int_t iPart, iMother;
278 Double_t weight = 0.;
280 Int_t maxPart = partArray->GetEntriesFast();
281 for(iPart=0; iPart<maxPart; iPart++){
282 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
283 iMother = part->GetFirstMother();
284 TParticle *mother = 0;
286 mother = gAlice->GetMCApp()->Particle(iMother);
287 pdgMother = mother->GetPdgCode();
290 pdgMother = part->GetPdgCode();
293 dNdy = fYieldArray[kGenPizero];
296 dNdy = fYieldArray[kGenEta];
299 dNdy = fYieldArray[kGenRho];
302 dNdy = fYieldArray[kGenOmega];
305 dNdy = fYieldArray[kGenEtaprime];
308 dNdy = fYieldArray[kGenPhi];
311 dNdy = fYieldArray[kGenJpsi];
317 weight = dNdy*part->GetWeight();
318 part->SetWeight(weight);
321 fHeader->SetNProduced(maxPart);
326 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
328 fHeader->SetPrimaryVertex(eventVertex);
330 gAlice->SetGenEventHeader(fHeader);