Removed memory leaks
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
CommitLineData
103ac317 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$ */
17
18//
19// Classe to create the coktail for physics with muons for pp at 14 TeV
20// The followoing sources:
21// jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia
22// The free parameeters are :
23// pp reaction cross-section
24// production cross-sections in pp collisions
b44c3901 25// July 07:added heavy quark production from AliGenCorrHF and heavy quark
26// production switched off in Pythia
f81a4aca 27// Aug. 07: added trigger cut on total momentum
103ac317 28
29#include <TObjArray.h>
30#include <TParticle.h>
31#include <TF1.h>
aaa95f22 32#include <TVirtualMC.h>
103ac317 33
34#include "AliGenCocktailEventHeader.h"
35
36#include "AliGenCocktailEntry.h"
37#include "AliGenMUONCocktailpp.h"
38#include "AliGenMUONlib.h"
39#include "AliGenParam.h"
40#include "AliMC.h"
41#include "AliRun.h"
42#include "AliStack.h"
aaa95f22 43#include "AliDecayer.h"
103ac317 44#include "AliLog.h"
b44c3901 45#include "AliGenCorrHF.h"
103ac317 46
47ClassImp(AliGenMUONCocktailpp)
48
49//________________________________________________________________________
50AliGenMUONCocktailpp::AliGenMUONCocktailpp()
1c56e311 51 :AliGenCocktail(),
aaa95f22 52 fDecayer(0),
53 fDecayModeResonance(kAll),
54 fDecayModePythia(kAll),
1c56e311 55 fMuonMultiplicity(0),
56 fMuonPtCut(0.),
f81a4aca 57 fMuonPCut(0.),
18bc0c8e 58 fMuonThetaMinCut(0.),
f81a4aca 59 fMuonThetaMaxCut(180.),
aaa95f22 60 fMuonOriginCut(-999.),
18bc0c8e 61 fNSucceded(0),
1c56e311 62 fNGenerated(0)
103ac317 63{
64// Constructor
103ac317 65}
93a2041b 66
103ac317 67//_________________________________________________________________________
68AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
69// Destructor
70{
71
72}
73
74//_________________________________________________________________________
8058ead1 75void AliGenMUONCocktailpp::CreateCocktail()
103ac317 76{
77// MinBias NN cross section @ pp 14 TeV -PR Vol II p:473
78 Double_t sigmaReaction = 0.070;
79
80// These limits are only used for renormalization of quarkonia cross section
81// Pythia events are generated in 4pi
82 Double_t ptMin = fPtMin;
83 Double_t ptMax = fPtMax;
84 Double_t yMin = fYMin;;
85 Double_t yMax = fYMax;;
86 Double_t phiMin = fPhiMin*180./TMath::Pi();
87 Double_t phiMax = fPhiMax*180./TMath::Pi();
88 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));
89
90// Ratios with respect to the reaction cross-section in the
91// kinematics limit of the MUONCocktail
92 Double_t ratiojpsi;
93 Double_t ratiopsiP;
94 Double_t ratioupsilon;
95 Double_t ratioupsilonP;
96 Double_t ratioupsilonPP;
b44c3901 97 Double_t ratioccbar;
98 Double_t ratiobbbar;
99
18bc0c8e 100// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
101// corrected from feed down of higher resonances
103ac317 102
18bc0c8e 103 Double_t sigmajpsi = 49.44e-6;
103ac317 104 Double_t sigmapsiP = 7.67e-6;
18bc0c8e 105 Double_t sigmaupsilon = 0.989e-6;
106 Double_t sigmaupsilonP = 0.502e-6;
103ac317 107 Double_t sigmaupsilonPP = 0.228e-6;
b44c3901 108 Double_t sigmaccbar = 11.2e-3;
109 Double_t sigmabbbar = 0.51e-3;
3285b419 110
111 Bool_t showSIGJPSI = kTRUE;
112 if (gSystem->Getenv("COCKTAIL_SIGJPSI")) {
113 sigmajpsi = atof(gSystem->Getenv("COCKTAIL_SIGJPSI"));
114 AliInfo("Cross-section for JPsi set from external value");
115 if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGJPSI = kFALSE;
116 }
117 Bool_t showSIGPSIP = kTRUE;
118 if (gSystem->Getenv("COCKTAIL_SIGPSIP")) {
119 sigmapsiP = atof(gSystem->Getenv("COCKTAIL_SIGPSIP"));
120 AliInfo("Cross-section for Psi-prime set from external value");
121 if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGPSIP = kFALSE;
122 }
123 Bool_t showSIGUPSI = kTRUE;
124 if (gSystem->Getenv("COCKTAIL_SIGUPSI")) {
125 sigmaupsilon = atof(gSystem->Getenv("COCKTAIL_SIGUPSI"));
126 AliInfo("Cross-section for Upsilon 1S set from external value");
127 if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSI = kFALSE;
128 }
129 Bool_t showSIGUPSIP = kTRUE;
130 if (gSystem->Getenv("COCKTAIL_SIGUPSIP")) {
131 sigmaupsilonP = atof(gSystem->Getenv("COCKTAIL_SIGUPSIP"));
132 AliInfo("Cross-section for Upsilon 2S set from external value");
133 if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSIP = kFALSE;
134 }
135 Bool_t showSIGUPSIPP = kTRUE;
136 if (gSystem->Getenv("COCKTAIL_SIGUPSIPP")) {
137 sigmaupsilonPP = atof(gSystem->Getenv("COCKTAIL_SIGUPSIPP"));
138 AliInfo("Cross-section for Upsilon 3S set from external value");
139 if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSIPP = kFALSE;
140 }
141 Bool_t showSIGCCBAR = kTRUE;
142 if (gSystem->Getenv("COCKTAIL_SIGCCBAR")) {
143 sigmaccbar = atof(gSystem->Getenv("COCKTAIL_SIGCCBAR"));
144 AliInfo("Cross-section for c-cbar set from external value");
145 if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGCCBAR = kFALSE;
146 }
147 Bool_t showSIGBBBAR = kTRUE;
148 if (gSystem->Getenv("COCKTAIL_SIGBBBAR")) {
149 sigmabbbar = atof(gSystem->Getenv("COCKTAIL_SIGBBBAR"));
150 AliInfo("Cross-section for b-bbar set from external value");
151 if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGBBBAR = kFALSE;
152 }
153
aaa95f22 154 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
155
18bc0c8e 156// Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
157//----------------------------------------------------------------------
103ac317 158// Jpsi generator
18bc0c8e 159 AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp", "Jpsi");
103ac317 160// first step: generation in 4pi
161 genjpsi->SetPtRange(0.,100.);
162 genjpsi->SetYRange(-8.,8.);
163 genjpsi->SetPhiRange(0.,360.);
aaa95f22 164 genjpsi->SetForceDecay(fDecayModeResonance);
165 if (!gMC) genjpsi->SetDecayer(fDecayer);
166
103ac317 167 genjpsi->Init(); // generation in 4pi
168 ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
3285b419 169 if (showSIGJPSI) {
170 AliInfo(Form("jpsi prod. cross-section in pp(14 TeV) %5.3g b",sigmajpsi));
171 AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
172 }
103ac317 173// second step: generation in selected kinematical range
174 genjpsi->SetPtRange(ptMin, ptMax);
175 genjpsi->SetYRange(yMin, yMax);
176 genjpsi->SetPhiRange(phiMin, phiMax);
177 genjpsi->Init(); // generation in selected kinematic range
178 AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
103ac317 179//------------------------------------------------------------------
180// Psi prime generator
18bc0c8e 181 AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
103ac317 182// first step: generation in 4pi
183 genpsiP->SetPtRange(0.,100.);
184 genpsiP->SetYRange(-8.,8.);
185 genpsiP->SetPhiRange(0.,360.);
aaa95f22 186 genpsiP->SetForceDecay(fDecayModeResonance);
187 if (!gMC) genpsiP->SetDecayer(fDecayer);
188
103ac317 189 genpsiP->Init(); // generation in 4pi
190 ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
3285b419 191 if (showSIGPSIP) {
192 AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP));
193 AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
194 }
103ac317 195// second step: generation in selected kinematical range
196 genpsiP->SetPtRange(ptMin, ptMax);
197 genpsiP->SetYRange(yMin, yMax);
198 genpsiP->SetPhiRange(phiMin, phiMax);
199 genpsiP->Init(); // generation in selected kinematic range
200 AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
103ac317 201//------------------------------------------------------------------
202// Upsilon 1S generator
18bc0c8e 203 AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
103ac317 204// first step: generation in 4pi
205 genupsilon->SetPtRange(0.,100.);
206 genupsilon->SetYRange(-8.,8.);
207 genupsilon->SetPhiRange(0.,360.);
aaa95f22 208 genupsilon->SetForceDecay(fDecayModeResonance);
209 if (!gMC) genupsilon->SetDecayer(fDecayer);
103ac317 210 genupsilon->Init(); // generation in 4pi
211 ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
3285b419 212 if (showSIGUPSI) {
213 AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon));
214 AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
215 }
103ac317 216// second step: generation in selected kinematical range
217 genupsilon->SetPtRange(ptMin, ptMax);
218 genupsilon->SetYRange(yMin, yMax);
219 genupsilon->SetPhiRange(phiMin, phiMax);
220 genupsilon->Init(); // generation in selected kinematical range
221 AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
103ac317 222//------------------------------------------------------------------
223// Upsilon 2S generator
18bc0c8e 224 AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
103ac317 225// first step: generation in 4pi
226 genupsilonP->SetPtRange(0.,100.);
227 genupsilonP->SetYRange(-8.,8.);
228 genupsilonP->SetPhiRange(0.,360.);
aaa95f22 229 genupsilonP->SetForceDecay(fDecayModeResonance);
230 if (!gMC) genupsilonP->SetDecayer(fDecayer);
103ac317 231 genupsilonP->Init(); // generation in 4pi
232 ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
3285b419 233 if (showSIGUPSIP) {
234 AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP));
235 AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
236 }
103ac317 237// second step: generation in the kinematical range
238 genupsilonP->SetPtRange(ptMin, ptMax);
239 genupsilonP->SetYRange(yMin, yMax);
240 genupsilonP->SetPhiRange(phiMin, phiMax);
241 genupsilonP->Init(); // generation in selected kinematical range
242 AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
103ac317 243
244//------------------------------------------------------------------
245// Upsilon 3S generator
18bc0c8e 246 AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
103ac317 247// first step: generation in 4pi
248 genupsilonPP->SetPtRange(0.,100.);
249 genupsilonPP->SetYRange(-8.,8.);
250 genupsilonPP->SetPhiRange(0.,360.);
aaa95f22 251 genupsilonPP->SetForceDecay(fDecayModeResonance);
252 if (!gMC) genupsilonPP->SetDecayer(fDecayer);
103ac317 253 genupsilonPP->Init(); // generation in 4pi
254 ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
3285b419 255 if (showSIGUPSIPP) {
256 AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP));
257 AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
258 }
103ac317 259// second step: generation in selected kinematical range
260 genupsilonPP->SetPtRange(ptMin, ptMax);
261 genupsilonPP->SetYRange(yMin, yMax);
262 genupsilonPP->SetPhiRange(phiMin, phiMax);
263 genupsilonPP->Init(); // generation in selected kinematical range
264 AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
b44c3901 265
266//------------------------------------------------------------------
267// Generator of charm
268
269 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4);
8058ead1 270 gencharm->SetMomentumRange(0,9999);
271 gencharm->SetForceDecay(kAll);
272 ratioccbar = sigmaccbar/sigmaReaction;
273 if (!gMC) gencharm->SetDecayer(fDecayer);
274 gencharm->Init();
3285b419 275 if (showSIGCCBAR) {
276 AliInfo(Form("c-cbar prod. cross-section in pp(14 TeV) %5.3g b",sigmaccbar));
277 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
278 }
8058ead1 279 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
b44c3901 280
103ac317 281//------------------------------------------------------------------
b44c3901 282// Generator of beauty
283
8058ead1 284 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5);
285 genbeauty->SetMomentumRange(0,9999);
286 genbeauty->SetForceDecay(kAll);
287 ratiobbbar = sigmabbbar/sigmaReaction;
288 if (!gMC) genbeauty->SetDecayer(fDecayer);
289 genbeauty->Init();
3285b419 290 if (showSIGBBBAR) {
291 AliInfo(Form("b-bbar prod. cross-section in pp(14 TeV) %5.3g b",sigmabbbar));
292 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
293 }
8058ead1 294 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
b44c3901 295
8058ead1 296//-------------------------------------------------------------------
103ac317 297// Pythia generator
8058ead1 298//
299// This has to go into the Config.C
300//
301// AliGenPythia *pythia = new AliGenPythia(1);
302// pythia->SetProcess(kPyMbMSEL1);
303// pythia->SetStrucFunc(kCTEQ5L);
304// pythia->SetEnergyCMS(14000.);
305// AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
306// Decay_t dt = gener->GetDecayModePythia();
307// pythia->SetForceDecay(dt);
308// pythia->SetPtRange(0.,100.);
309// pythia->SetYRange(-8.,8.);
310// pythia->SetPhiRange(0.,360.);
311// pythia->SetPtHard(2.76,-1.0);
312// pythia->SwitchHFOff();
313// pythia->Init();
314// AddGenerator(pythia,"Pythia",1);
315}
103ac317 316
8058ead1 317void AliGenMUONCocktailpp::Init()
318{
319 //
320 // Initialisation
aaa95f22 321 TIter next(fEntries);
322 AliGenCocktailEntry *entry;
323 if (fStack) {
324 while((entry = (AliGenCocktailEntry*)next())) {
325 entry->Generator()->SetStack(fStack);
326 }
327 }
103ac317 328}
329
330//_________________________________________________________________________
331void AliGenMUONCocktailpp::Generate()
332{
333// Generate event
334 TIter next(fEntries);
335 AliGenCocktailEntry *entry = 0;
336 AliGenCocktailEntry *preventry = 0;
337 AliGenerator* gen = 0;
338
d94af0c1 339 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
103ac317 340
341// Generate the vertex position used by all generators
342 if(fVertexSmear == kPerEvent) Vertex();
343
344// Loop on primordialTrigger:
345// minimum muon multiplicity above a pt cut in a theta acceptance region
346
347 Bool_t primordialTrigger = kFALSE;
348 while(!primordialTrigger) {
349 //Reseting stack
33c3c91a 350 AliRunLoader * runloader = AliRunLoader::Instance();
103ac317 351 if (runloader)
352 if (runloader->Stack())
34da8494 353 runloader->Stack()->Clean();
103ac317 354 // Loop over generators and generate events
355 Int_t igen = 0;
356 Int_t npart = 0;
357 const char* genName = 0;
358 while((entry = (AliGenCocktailEntry*)next())) {
359 gen = entry->Generator();
360 genName = entry->GetName();
361 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
362
363 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
364 gRandom->Poisson(entry->Rate());
365
366 if (npart > 0) {
367 igen++;
368 if (igen == 1) entry->SetFirst(0);
369 else entry->SetFirst((partArray->GetEntriesFast())+1);
370
371 gen->SetNumberParticles(npart);
372 gen->Generate();
373 entry->SetLast(partArray->GetEntriesFast());
374 preventry = entry;
375 }
376 }
377 next.Reset();
378
379// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
380// in the muon spectrometer acceptance
381 Int_t iPart;
382 fNGenerated++;
b44c3901 383 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22 384 for(iPart=0; iPart<maxPart; iPart++){
103ac317 385
aaa95f22 386 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
387 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
388 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
389 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
390 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
f81a4aca 391 (part->Pt()>fMuonPtCut) &&
392 (part->P()>fMuonPCut)) {
aaa95f22 393 numberOfMuons++;
103ac317 394 }
aaa95f22 395 }
396 }
103ac317 397 if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
398 }
399 fNSucceded++;
400
aaa95f22 401// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 402 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
403}
404
405