- correct minor issues with forcred gamma decays, and speed up
[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(0),
61   fCentrality(0),
62   fV2Systematic(0),
63   fForceConv(kFALSE),
64   fHeaviestParticle(kGenJpsi)
65 {
66   // Constructor
67
68 }
69
70 //_________________________________________________________________________
71 AliGenEMCocktail::~AliGenEMCocktail()
72 {
73   // Destructor
74
75 }
76
77 //_________________________________________________________________________
78 void AliGenEMCocktail::CreateCocktail()
79 {
80   // create and add sources to the cocktail
81
82   fDecayer->SetForceDecay(fDecayMode);
83   fDecayer->ForceDecay();
84
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));
94
95   //Initialize user selection for Pt Parameterization and centrality:
96   AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
97
98   // Create and add electron sources to the generator
99   // pizero
100   if(fHeaviestParticle<kGenPizero)return;
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[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
110 #else
111   fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
112 #endif
113
114   // eta  
115   if(fHeaviestParticle<kGenEta)return;
116   AliGenParam *geneta=0;
117   Char_t nameEta[10];    
118   snprintf(nameEta,10,"Eta");    
119   geneta = new AliGenParam(fNPart/0.825, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
120   geneta->SetYRange(fYMin/0.825, fYMax/0.825);
121   AddSource2Generator(nameEta,geneta);
122   TF1 *fPtEta = geneta->GetPt();
123 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
124   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
125 #else
126   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
127 #endif
128
129   // rho  
130   if(fHeaviestParticle<kGenRho)return;
131   AliGenParam *genrho=0;
132   Char_t nameRho[10];    
133   snprintf(nameRho,10,"Rho");    
134   genrho = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
135   genrho->SetYRange(fYMin/0.775, fYMax/0.775);
136   AddSource2Generator(nameRho,genrho);
137   TF1 *fPtRho = genrho->GetPt();
138 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
139   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
140 #else
141   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
142 #endif
143   
144   // omega
145   if(fHeaviestParticle<kGenOmega)return;
146   AliGenParam *genomega=0;
147   Char_t nameOmega[10];    
148   snprintf(nameOmega,10,"Omega");    
149   genomega = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
150   genomega->SetYRange(fYMin/0.775, fYMax/0.775);
151   AddSource2Generator(nameOmega,genomega);
152   TF1 *fPtOmega = genomega->GetPt();
153 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
154   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
155 #else
156   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
157 #endif
158
159   // etaprime
160   if(fHeaviestParticle<kGenEtaprime)return;
161   AliGenParam *genetaprime=0;
162   Char_t nameEtaprime[10];    
163   snprintf(nameEtaprime,10,"Etaprime");    
164   genetaprime = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
165   genetaprime->SetYRange(fYMin/0.725, fYMax/0.725);
166   AddSource2Generator(nameEtaprime,genetaprime);
167   TF1 *fPtEtaprime = genetaprime->GetPt();
168 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
169   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
170 #else
171   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
172 #endif
173
174   // phi  
175   if(fHeaviestParticle<kGenPhi)return;
176   AliGenParam *genphi=0;
177   Char_t namePhi[10];    
178   snprintf(namePhi,10,"Phi");    
179   genphi = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
180   genphi->SetYRange(fYMin/0.725, fYMax/0.725);
181   AddSource2Generator(namePhi,genphi);
182   TF1 *fPtPhi = genphi->GetPt();
183 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
184   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
185 #else
186   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
187 #endif
188
189   // jpsi  
190   if(fHeaviestParticle<kGenJpsi)return;
191   AliGenParam *genjpsi=0;
192   Char_t nameJpsi[10];    
193   snprintf(nameJpsi,10,"Jpsi");    
194   genjpsi = new AliGenParam(fNPart/0.525, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
195   genjpsi->SetYRange(fYMin/0.525, fYMax/0.525);
196   AddSource2Generator(nameJpsi,genjpsi);
197   TF1 *fPtJpsi = genjpsi->GetPt();
198 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
199   fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
200 #else
201   fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
202 #endif
203
204   // prompt gamma
205   if(fDecayMode==kGammaEM){
206     if(fHeaviestParticle<kGenPromptRealGamma)return;
207     TDatabasePDG::Instance()->AddParticle("PromptRealGamma","PromptRealGamma",0,true,0,0,"GaugeBoson",221000);
208     //gMC->DefineParticle(221000, "PromptGamma", kPTGamma, 0, 0, 0,"Gamma", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kFALSE);
209     AliGenParam *genPromptRealG=0;
210     Char_t namePromptRealG[10];    
211     snprintf(namePromptRealG,10,"PromptRealGamma");    
212     genPromptRealG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kPromptRealGamma, "DUMMY");
213     genPromptRealG->SetYRange(fYMin, fYMax);
214     AddSource2Generator(namePromptRealG,genPromptRealG);
215     TF1 *fPtPromptRealG = genPromptRealG->GetPt();
216 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
217     fYieldArray[kGenPromptRealGamma] = fPtPromptRealG->Integral(fPtMin,fPtMax,1.e-6);
218 #else
219     fYieldArray[kGenPromptRealGamma] = fPtPromptRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
220 #endif
221
222     if(fHeaviestParticle<kGenThermRealGamma)return;
223     TDatabasePDG::Instance()->AddParticle("ThermRealGamma","ThermRealGamma",0,true,0,0,"GaugeBoson",222000);
224     //gMC->DefineParticle(221000, "ThermGamma", kPTGamma, 0, 0, 0,"Gamma", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kFALSE);
225     AliGenParam *genThermRealG=0;
226     Char_t nameThermRealG[10];    
227     snprintf(nameThermRealG,10,"ThermRealGamma");    
228     genThermRealG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kThermRealGamma, "DUMMY");
229     genThermRealG->SetYRange(fYMin, fYMax);
230     AddSource2Generator(nameThermRealG,genThermRealG);
231     TF1 *fPtThermRealG = genThermRealG->GetPt();
232 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
233     fYieldArray[kGenThermRealGamma] = fPtThermRealG->Integral(fPtMin,fPtMax,1.e-6);
234 #else
235     fYieldArray[kGenThermRealGamma] = fPtThermRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
236 #endif
237
238     if(fHeaviestParticle<kGenPromptVirtGamma)return;
239     TDatabasePDG::Instance()->AddParticle("PromptVirtGamma","PromptVirtGamma",0,true,0,0,"GaugeBoson",223000);
240     AliGenParam *genPromptVirtG=0;
241     Char_t namePromptVirtG[10];    
242     snprintf(namePromptVirtG,10,"PromptVirtGamma");    
243     genPromptVirtG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kPromptVirtGamma, "DUMMY");
244     genPromptVirtG->SetYRange(fYMin, fYMax);
245     AddSource2Generator(namePromptVirtG,genPromptVirtG);
246     TF1 *fPtPromptVirtG = genPromptVirtG->GetPt();
247 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
248     fYieldArray[kGenPromptVirtGamma] = fPtPromptVirtG->Integral(fPtMin,fPtMax,1.e-6);
249 #else
250     fYieldArray[kGenPromptVirtGamma] = fPtPromptVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
251 #endif
252
253     if(fHeaviestParticle<kGenThermVirtGamma)return;
254     TDatabasePDG::Instance()->AddParticle("ThermVirtGamma","ThermVirtGamma",0,true,0,0,"GaugeBoson",224000);
255     AliGenParam *genThermVirtG=0;
256     Char_t nameThermVirtG[10];    
257     snprintf(nameThermVirtG,10,"ThermVirtGamma");    
258     genThermVirtG = new AliGenParam(fNPart*0.5, new AliGenEMlib(), AliGenEMlib::kThermVirtGamma, "DUMMY");
259     genThermVirtG->SetYRange(fYMin, fYMax);
260     AddSource2Generator(nameThermVirtG,genThermVirtG);
261     TF1 *fPtThermVirtG = genThermVirtG->GetPt();
262 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
263     fYieldArray[kGenThermVirtGamma] = fPtThermVirtG->Integral(fPtMin,fPtMax,1.e-6);
264 #else
265     fYieldArray[kGenThermVirtGamma] = fPtThermVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
266 #endif
267   }
268   
269 }
270
271 //-------------------------------------------------------------------
272 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource, 
273                                          AliGenParam* const genSource)
274 {
275 // add sources to the cocktail
276   Double_t phiMin = fPhiMin*180./TMath::Pi();
277   Double_t phiMax = fPhiMax*180./TMath::Pi();
278
279   genSource->SetPtRange(fPtMin, fPtMax);  
280   genSource->SetPhiRange(phiMin, phiMax);
281   genSource->SetWeighting(fWeightingMode);
282   genSource->SetForceGammaConversion(fForceConv);
283   if (!gMC) genSource->SetDecayer(fDecayer);  
284   genSource->Init();
285     
286   AddGenerator(genSource,nameSource,1.); // Adding Generator    
287 }
288
289 //-------------------------------------------------------------------
290 void AliGenEMCocktail::Init()
291 {
292   // Initialisation
293   TIter next(fEntries);
294   AliGenCocktailEntry *entry;
295   if (fStack) {
296     while((entry = (AliGenCocktailEntry*)next())) {
297       entry->Generator()->SetStack(fStack);
298     }
299   }
300 }
301
302 //_________________________________________________________________________
303 void AliGenEMCocktail::Generate()
304 {
305   // Generate event 
306   TIter next(fEntries);
307   AliGenCocktailEntry *entry = 0;
308   AliGenCocktailEntry *preventry = 0;
309   AliGenerator* gen = 0;
310
311   if (fHeader) delete fHeader;
312   fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
313
314   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
315     
316   // Generate the vertex position used by all generators    
317   if(fVertexSmear == kPerEvent) Vertex();
318
319   //Reseting stack
320   AliRunLoader * runloader = AliRunLoader::Instance();
321   if (runloader)
322     if (runloader->Stack())
323       runloader->Stack()->Clean();
324   
325   // Loop over generators and generate events
326   Int_t igen = 0;
327   Float_t evPlane;
328   Rndm(&evPlane,1);
329   evPlane*=TMath::Pi()*2;
330   while((entry = (AliGenCocktailEntry*)next())) {
331     gen = entry->Generator();
332     gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
333     
334     if (fNPart > 0) {
335       igen++;   
336       if (igen == 1) entry->SetFirst(0);                
337       else  entry->SetFirst((partArray->GetEntriesFast())+1);
338       gen->SetEventPlane(evPlane);
339       gen->Generate();
340       entry->SetLast(partArray->GetEntriesFast());
341       preventry = entry;
342     }
343   }  
344   next.Reset();
345
346   // Setting weights for proper absolute normalization
347   Int_t iPart, iMother;
348   Int_t pdgMother = 0;
349   Double_t weight = 0.;
350   Double_t dNdy = 0.;
351   Int_t maxPart = partArray->GetEntriesFast();
352   for(iPart=0; iPart<maxPart; iPart++){      
353     TParticle *part = gAlice->GetMCApp()->Particle(iPart);
354     iMother = part->GetFirstMother();
355     TParticle *mother = 0;
356     if (iMother>=0){
357       mother = gAlice->GetMCApp()->Particle(iMother);
358       pdgMother = mother->GetPdgCode();
359     }
360     else
361       pdgMother = part->GetPdgCode();
362
363     switch (pdgMother){
364     case 111:
365       dNdy = fYieldArray[kGenPizero];
366       break;
367     case 221:
368       dNdy = fYieldArray[kGenEta];
369       break;
370     case 113:
371       dNdy = fYieldArray[kGenRho];
372       break;
373     case 223:
374       dNdy = fYieldArray[kGenOmega];
375       break;
376     case 331:
377       dNdy = fYieldArray[kGenEtaprime];
378       break;
379     case 333:
380       dNdy = fYieldArray[kGenPhi];
381       break;
382     case 443:
383       dNdy = fYieldArray[kGenJpsi];
384       break;
385     default:
386       dNdy = 0.;
387     }
388
389     switch (pdgMother){
390     case 221000:
391       dNdy = fYieldArray[kGenPromptRealGamma];
392       break;
393     case 223000:
394       dNdy = fYieldArray[kGenPromptVirtGamma];
395       break;
396     case 222000:
397       dNdy = fYieldArray[kGenThermRealGamma];
398       break;      
399     case 224000:
400       dNdy = fYieldArray[kGenThermVirtGamma];
401       break;   
402     }
403     weight = dNdy*part->GetWeight();
404     part->SetWeight(weight);
405   }     
406   
407   fHeader->SetNProduced(maxPart);
408
409
410   TArrayF eventVertex;
411   eventVertex.Set(3);
412   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
413   
414   fHeader->SetPrimaryVertex(eventVertex);
415
416   gAlice->SetGenEventHeader(fHeader);
417 }