add histogram to check invariant mass for pairs inside isolation cone; change printf...
[u/mrichter/AliRoot.git] / EVGEN / AliGenEMCocktail.cxx
CommitLineData
e40b9538 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
4ae1c9f0 18// Class to create the cocktail for physics with electrons, di-electrons,
e40b9538 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
4ae1c9f0 28
e40b9538 29#include <TObjArray.h>
30#include <TParticle.h>
31#include <TF1.h>
32#include <TVirtualMC.h>
33#include <TPDGCode.h>
71443190 34#include <TDatabasePDG.h>
e40b9538 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
50ClassImp(AliGenEMCocktail)
51
52//________________________________________________________________________
53AliGenEMCocktail::AliGenEMCocktail()
4ae1c9f0 54:AliGenCocktail(),
e40b9538 55 fDecayer(0),
56 fDecayMode(kAll),
57 fWeightingMode(kNonAnalog),
58 fNPart(1000),
4ae1c9f0 59 fYieldArray(),
1b589c8a 60 fPtSelect(AliGenEMlib::kPizero7TeVpp),
61 fCentrality(AliGenEMlib::kpp),
62 fV2Systematic(AliGenEMlib::kNoV2Sys),
4ae1c9f0 63 fForceConv(kFALSE),
1b589c8a 64 fSelectedParticles(kGenPizero)
e40b9538 65{
4ae1c9f0 66 // Constructor
1b589c8a 67 SetHeaviestHadron(kGenPhi);
e40b9538 68}
69
70//_________________________________________________________________________
71AliGenEMCocktail::~AliGenEMCocktail()
72{
4ae1c9f0 73 // Destructor
e40b9538 74
75}
76
77//_________________________________________________________________________
1b589c8a 78void 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//_________________________________________________________________________
e40b9538 87void AliGenEMCocktail::CreateCocktail()
88{
4ae1c9f0 89 // create and add sources to the cocktail
e40b9538 90
91 fDecayer->SetForceDecay(fDecayMode);
92 fDecayer->ForceDecay();
93
4ae1c9f0 94 // Set kinematic limits
e40b9538 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
4ae1c9f0 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
4ae1c9f0 108 // pizero
6a8b015a 109 if(fSelectedParticles&kGenPizero){
4ae1c9f0 110 AliGenParam *genpizero=0;
e40b9538 111 Char_t namePizero[10];
4ae1c9f0 112 snprintf(namePizero,10,"Pizero");
1b589c8a 113 //fNPart/0.925: increase number of particles so that we have the chosen number of particles in the chosen eta range
71443190 114 genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
1b589c8a 115 //fYMin/0.925: increase eta range, so that the electron yield is constant (<5% change) over the chosen eta range
71443190 116 genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
e40b9538 117 AddSource2Generator(namePizero,genpizero);
118 TF1 *fPtPizero = genpizero->GetPt();
faf00b0d 119#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
6a8b015a 120 fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
faf00b0d 121#else
6a8b015a 122 fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d 123#endif
6a8b015a 124 }
4ae1c9f0 125
126 // eta
6a8b015a 127 if(fSelectedParticles&kGenEta){
4ae1c9f0 128 AliGenParam *geneta=0;
e40b9538 129 Char_t nameEta[10];
4ae1c9f0 130 snprintf(nameEta,10,"Eta");
71443190 131 geneta = new AliGenParam(fNPart/0.825, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
132 geneta->SetYRange(fYMin/0.825, fYMax/0.825);
e40b9538 133 AddSource2Generator(nameEta,geneta);
134 TF1 *fPtEta = geneta->GetPt();
4ae1c9f0 135#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
6a8b015a 136 fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
4ae1c9f0 137#else
6a8b015a 138 fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d 139#endif
6a8b015a 140 }
4ae1c9f0 141
142 // rho
6a8b015a 143 if(fSelectedParticles&kGenRho){
4ae1c9f0 144 AliGenParam *genrho=0;
e40b9538 145 Char_t nameRho[10];
4ae1c9f0 146 snprintf(nameRho,10,"Rho");
71443190 147 genrho = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
148 genrho->SetYRange(fYMin/0.775, fYMax/0.775);
e40b9538 149 AddSource2Generator(nameRho,genrho);
150 TF1 *fPtRho = genrho->GetPt();
faf00b0d 151#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
6a8b015a 152 fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
faf00b0d 153#else
6a8b015a 154 fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d 155#endif
6a8b015a 156 }
4ae1c9f0 157
158 // omega
6a8b015a 159 if(fSelectedParticles&kGenOmega){
4ae1c9f0 160 AliGenParam *genomega=0;
e40b9538 161 Char_t nameOmega[10];
4ae1c9f0 162 snprintf(nameOmega,10,"Omega");
71443190 163 genomega = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
164 genomega->SetYRange(fYMin/0.775, fYMax/0.775);
e40b9538 165 AddSource2Generator(nameOmega,genomega);
166 TF1 *fPtOmega = genomega->GetPt();
faf00b0d 167#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
6a8b015a 168 fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
faf00b0d 169#else
6a8b015a 170 fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d 171#endif
6a8b015a 172 }
4ae1c9f0 173
174 // etaprime
6a8b015a 175 if(fSelectedParticles&kGenEtaprime){
4ae1c9f0 176 AliGenParam *genetaprime=0;
e40b9538 177 Char_t nameEtaprime[10];
4ae1c9f0 178 snprintf(nameEtaprime,10,"Etaprime");
71443190 179 genetaprime = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
180 genetaprime->SetYRange(fYMin/0.725, fYMax/0.725);
e40b9538 181 AddSource2Generator(nameEtaprime,genetaprime);
182 TF1 *fPtEtaprime = genetaprime->GetPt();
faf00b0d 183#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
6a8b015a 184 fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
faf00b0d 185#else
6a8b015a 186 fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d 187#endif
6a8b015a 188 }
4ae1c9f0 189
190 // phi
6a8b015a 191 if(fSelectedParticles&kGenPhi){
4ae1c9f0 192 AliGenParam *genphi=0;
e40b9538 193 Char_t namePhi[10];
4ae1c9f0 194 snprintf(namePhi,10,"Phi");
71443190 195 genphi = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
196 genphi->SetYRange(fYMin/0.725, fYMax/0.725);
e40b9538 197 AddSource2Generator(namePhi,genphi);
198 TF1 *fPtPhi = genphi->GetPt();
faf00b0d 199#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
6a8b015a 200 fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
faf00b0d 201#else
6a8b015a 202 fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
faf00b0d 203#endif
6a8b015a 204 }
4ae1c9f0 205
206 // jpsi
6a8b015a 207 if(fSelectedParticles&kGenJpsi){
4ae1c9f0 208 AliGenParam *genjpsi=0;
209 Char_t nameJpsi[10];
210 snprintf(nameJpsi,10,"Jpsi");
71443190 211 genjpsi = new AliGenParam(fNPart/0.525, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
212 genjpsi->SetYRange(fYMin/0.525, fYMax/0.525);
4ae1c9f0 213 AddSource2Generator(nameJpsi,genjpsi);
214 TF1 *fPtJpsi = genjpsi->GetPt();
215#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
6a8b015a 216 fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
4ae1c9f0 217#else
6a8b015a 218 fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
71443190 219#endif
6a8b015a 220 }
71443190 221
6a8b015a 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();
71443190 232#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
6a8b015a 233 fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
71443190 234#else
6a8b015a 235 fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
71443190 236#endif
6a8b015a 237 }
71443190 238
6a8b015a 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();
71443190 248#if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
6a8b015a 249 fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,1.e-6);
71443190 250#else
6a8b015a 251 fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
71443190 252#endif
253 }
e40b9538 254}
255
256//-------------------------------------------------------------------
257void 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);
e40b9538 265 genSource->SetPhiRange(phiMin, phiMax);
266 genSource->SetWeighting(fWeightingMode);
4ae1c9f0 267 genSource->SetForceGammaConversion(fForceConv);
2942f542 268 if (!TVirtualMC::GetMC()) genSource->SetDecayer(fDecayer);
e40b9538 269 genSource->Init();
270
271 AddGenerator(genSource,nameSource,1.); // Adding Generator
272}
273
274//-------------------------------------------------------------------
275void AliGenEMCocktail::Init()
276{
4ae1c9f0 277 // Initialisation
e40b9538 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//_________________________________________________________________________
288void AliGenEMCocktail::Generate()
289{
4ae1c9f0 290 // Generate event
e40b9538 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
4ae1c9f0 301 // Generate the vertex position used by all generators
e40b9538 302 if(fVertexSmear == kPerEvent) Vertex();
303
4ae1c9f0 304 //Reseting stack
e40b9538 305 AliRunLoader * runloader = AliRunLoader::Instance();
306 if (runloader)
307 if (runloader->Stack())
308 runloader->Stack()->Clean();
309
4ae1c9f0 310 // Loop over generators and generate events
e40b9538 311 Int_t igen = 0;
6078e216 312 Float_t evPlane;
313 Rndm(&evPlane,1);
314 evPlane*=TMath::Pi()*2;
e40b9538 315 while((entry = (AliGenCocktailEntry*)next())) {
316 gen = entry->Generator();
e40b9538 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);
6078e216 323 gen->SetEventPlane(evPlane);
e40b9538 324 gen->Generate();
325 entry->SetLast(partArray->GetEntriesFast());
326 preventry = entry;
327 }
328 }
329 next.Reset();
330
4ae1c9f0 331 // Setting weights for proper absolute normalization
e40b9538 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();
71443190 347
e40b9538 348 switch (pdgMother){
349 case 111:
6a8b015a 350 dNdy = fYieldArray[kPizero];
e40b9538 351 break;
352 case 221:
6a8b015a 353 dNdy = fYieldArray[kEta];
e40b9538 354 break;
355 case 113:
6a8b015a 356 dNdy = fYieldArray[kRho];
e40b9538 357 break;
358 case 223:
6a8b015a 359 dNdy = fYieldArray[kOmega];
e40b9538 360 break;
361 case 331:
6a8b015a 362 dNdy = fYieldArray[kEtaprime];
e40b9538 363 break;
364 case 333:
6a8b015a 365 dNdy = fYieldArray[kPhi];
e40b9538 366 break;
4ae1c9f0 367 case 443:
6a8b015a 368 dNdy = fYieldArray[kJpsi];
369 break;
370 case 22:
371 dNdy = fYieldArray[kDirectRealGamma];
4ae1c9f0 372 break;
6a8b015a 373 case 220000:
374 dNdy = fYieldArray[kDirectVirtGamma];
375 break;
376
e40b9538 377 default:
378 dNdy = 0.;
379 }
71443190 380
e40b9538 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}