Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / EVGEN / AliGenEMCocktail.cxx
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 "AliDecayer.h"
46 #include "AliDecayerPythia.h"
47 #include "AliLog.h"
48 #include "AliGenCorrHF.h"
49
50 ClassImp(AliGenEMCocktail)  
51   
52 //________________________________________________________________________
53 AliGenEMCocktail::AliGenEMCocktail()
54 :AliGenCocktail(),
55    fDecayer(0),
56    fDecayMode(kAll),
57    fWeightingMode(kNonAnalog),
58    fNPart(1000),
59   fYieldArray(),
60   fPtSelect(AliGenEMlib::kPizero7TeVpp),
61   fCentrality(AliGenEMlib::kpp),
62   fV2Systematic(AliGenEMlib::kNoV2Sys),
63   fForceConv(kFALSE),
64   fSelectedParticles(kGenPizero)
65 {
66   // Constructor
67   SetHeaviestHadron(kGenPhi);
68 }
69
70 //_________________________________________________________________________
71 AliGenEMCocktail::~AliGenEMCocktail()
72 {
73   // Destructor
74
75 }
76
77 //_________________________________________________________________________
78 void AliGenEMCocktail::SetHeaviestHadron(ParticeGenerator_t part)
79 {
80   Int_t val=kGenPizero;
81   while(val<part) val|=val<<1;
82
83   fSelectedParticles=val;
84 }
85
86 //_________________________________________________________________________
87 void AliGenEMCocktail::CreateCocktail()
88 {
89   // create and add sources to the cocktail
90
91   fDecayer->SetForceDecay(fDecayMode);
92   fDecayer->ForceDecay();
93
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));
103
104   //Initialize user selection for Pt Parameterization and centrality:
105   AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
106
107   // Create and add electron sources to the generator
108   // pizero
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);
121 #else
122     fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
123 #endif
124   }
125
126   // eta  
127   if(fSelectedParticles&kGenEta){
128   AliGenParam *geneta=0;
129   Char_t nameEta[10];    
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);
137 #else
138     fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
139 #endif
140   }
141
142   // rho  
143   if(fSelectedParticles&kGenRho){
144   AliGenParam *genrho=0;
145   Char_t nameRho[10];    
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);
153 #else
154     fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
155 #endif
156   }
157   
158   // omega
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);
169 #else
170     fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
171 #endif
172   }
173
174   // etaprime
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);
185 #else
186     fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
187 #endif
188   }
189
190   // phi  
191   if(fSelectedParticles&kGenPhi){
192   AliGenParam *genphi=0;
193   Char_t namePhi[10];    
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);
201 #else
202     fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
203 #endif
204   }
205
206   // jpsi  
207   if(fSelectedParticles&kGenJpsi){
208   AliGenParam *genjpsi=0;
209   Char_t nameJpsi[10];    
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);
217 #else
218     fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
219 #endif
220   }
221
222   // direct gamma
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);
234 #else
235     fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
236 #endif
237   }
238
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);
250 #else
251     fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
252 #endif
253   }
254 }
255
256 //-------------------------------------------------------------------
257 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource, 
258                                          AliGenParam* const genSource)
259 {
260 // add sources to the cocktail
261   Double_t phiMin = fPhiMin*180./TMath::Pi();
262   Double_t phiMax = fPhiMax*180./TMath::Pi();
263
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);  
269   genSource->Init();
270     
271   AddGenerator(genSource,nameSource,1.); // Adding Generator    
272 }
273
274 //-------------------------------------------------------------------
275 void AliGenEMCocktail::Init()
276 {
277   // Initialisation
278   TIter next(fEntries);
279   AliGenCocktailEntry *entry;
280   if (fStack) {
281     while((entry = (AliGenCocktailEntry*)next())) {
282       entry->Generator()->SetStack(fStack);
283     }
284   }
285 }
286
287 //_________________________________________________________________________
288 void AliGenEMCocktail::Generate()
289 {
290   // Generate event 
291   TIter next(fEntries);
292   AliGenCocktailEntry *entry = 0;
293   AliGenCocktailEntry *preventry = 0;
294   AliGenerator* gen = 0;
295
296   if (fHeader) delete fHeader;
297   fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
298
299   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
300     
301   // Generate the vertex position used by all generators    
302   if(fVertexSmear == kPerEvent) Vertex();
303
304   //Reseting stack
305   AliRunLoader * runloader = AliRunLoader::Instance();
306   if (runloader)
307     if (runloader->Stack())
308       runloader->Stack()->Clean();
309   
310   // Loop over generators and generate events
311   Int_t igen = 0;
312   Float_t evPlane;
313   Rndm(&evPlane,1);
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));
318     
319     if (fNPart > 0) {
320       igen++;   
321       if (igen == 1) entry->SetFirst(0);                
322       else  entry->SetFirst((partArray->GetEntriesFast())+1);
323       gen->SetEventPlane(evPlane);
324       gen->Generate();
325       entry->SetLast(partArray->GetEntriesFast());
326       preventry = entry;
327     }
328   }  
329   next.Reset();
330
331   // Setting weights for proper absolute normalization
332   Int_t iPart, iMother;
333   Int_t pdgMother = 0;
334   Double_t weight = 0.;
335   Double_t dNdy = 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;
341     if (iMother>=0){
342       mother = gAlice->GetMCApp()->Particle(iMother);
343       pdgMother = mother->GetPdgCode();
344     }
345     else
346       pdgMother = part->GetPdgCode();
347
348     switch (pdgMother){
349     case 111:
350       dNdy = fYieldArray[kPizero];
351       break;
352     case 221:
353       dNdy = fYieldArray[kEta];
354       break;
355     case 113:
356       dNdy = fYieldArray[kRho];
357       break;
358     case 223:
359       dNdy = fYieldArray[kOmega];
360       break;
361     case 331:
362       dNdy = fYieldArray[kEtaprime];
363       break;
364     case 333:
365       dNdy = fYieldArray[kPhi];
366       break;
367     case 443:
368       dNdy = fYieldArray[kJpsi];
369       break;
370     case 22:
371       dNdy = fYieldArray[kDirectRealGamma];
372       break;
373     case 220000:
374       dNdy = fYieldArray[kDirectVirtGamma];
375       break;
376       
377     default:
378       dNdy = 0.;
379     }
380
381     weight = dNdy*part->GetWeight();
382     part->SetWeight(weight);
383   }     
384   
385   fHeader->SetNProduced(maxPart);
386
387
388   TArrayF eventVertex;
389   eventVertex.Set(3);
390   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
391   
392   fHeader->SetPrimaryVertex(eventVertex);
393
394   gAlice->SetGenEventHeader(fHeader);
395 }