Some clean-up and new method AliStack::Clean needed by AliGenMUONCocktail
[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
103ac317 27
28#include <TObjArray.h>
29#include <TParticle.h>
30#include <TF1.h>
aaa95f22 31#include <TVirtualMC.h>
103ac317 32
33#include "AliGenCocktailEventHeader.h"
34
35#include "AliGenCocktailEntry.h"
36#include "AliGenMUONCocktailpp.h"
37#include "AliGenMUONlib.h"
38#include "AliGenParam.h"
39#include "AliMC.h"
40#include "AliRun.h"
41#include "AliStack.h"
aaa95f22 42#include "AliDecayer.h"
103ac317 43#include "AliLog.h"
44#include "AliGenPythia.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 fTotalRate(0),
56 fMuonMultiplicity(0),
57 fMuonPtCut(0.),
18bc0c8e 58 fMuonThetaMinCut(0.),
1c56e311 59 fMuonThetaMaxCut(0.),
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//_________________________________________________________________________
75void AliGenMUONCocktailpp::Init()
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;
103ac317 110
aaa95f22 111 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
112
18bc0c8e 113// Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
114//----------------------------------------------------------------------
103ac317 115// Jpsi generator
18bc0c8e 116 AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp", "Jpsi");
103ac317 117// first step: generation in 4pi
118 genjpsi->SetPtRange(0.,100.);
119 genjpsi->SetYRange(-8.,8.);
120 genjpsi->SetPhiRange(0.,360.);
aaa95f22 121 genjpsi->SetForceDecay(fDecayModeResonance);
122 if (!gMC) genjpsi->SetDecayer(fDecayer);
123
103ac317 124 genjpsi->Init(); // generation in 4pi
125 ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
126 AliInfo(Form("jpsi prod. cross-section in pp(14 TeV) %5.3g b",sigmajpsi));
127 AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
128// second step: generation in selected kinematical range
129 genjpsi->SetPtRange(ptMin, ptMax);
130 genjpsi->SetYRange(yMin, yMax);
131 genjpsi->SetPhiRange(phiMin, phiMax);
132 genjpsi->Init(); // generation in selected kinematic range
133 AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
134 fTotalRate+=ratiojpsi;
135
136//------------------------------------------------------------------
137// Psi prime generator
18bc0c8e 138 AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
103ac317 139// first step: generation in 4pi
140 genpsiP->SetPtRange(0.,100.);
141 genpsiP->SetYRange(-8.,8.);
142 genpsiP->SetPhiRange(0.,360.);
aaa95f22 143 genpsiP->SetForceDecay(fDecayModeResonance);
144 if (!gMC) genpsiP->SetDecayer(fDecayer);
145
103ac317 146 genpsiP->Init(); // generation in 4pi
147 ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
148 AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP));
149 AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
150// second step: generation in selected kinematical range
151 genpsiP->SetPtRange(ptMin, ptMax);
152 genpsiP->SetYRange(yMin, yMax);
153 genpsiP->SetPhiRange(phiMin, phiMax);
154 genpsiP->Init(); // generation in selected kinematic range
155 AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
156 fTotalRate+=ratiopsiP;
157
158//------------------------------------------------------------------
159// Upsilon 1S generator
18bc0c8e 160 AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
103ac317 161// first step: generation in 4pi
162 genupsilon->SetPtRange(0.,100.);
163 genupsilon->SetYRange(-8.,8.);
164 genupsilon->SetPhiRange(0.,360.);
aaa95f22 165 genupsilon->SetForceDecay(fDecayModeResonance);
166 if (!gMC) genupsilon->SetDecayer(fDecayer);
103ac317 167 genupsilon->Init(); // generation in 4pi
168 ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
169 AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon));
170 AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
171// second step: generation in selected kinematical range
172 genupsilon->SetPtRange(ptMin, ptMax);
173 genupsilon->SetYRange(yMin, yMax);
174 genupsilon->SetPhiRange(phiMin, phiMax);
175 genupsilon->Init(); // generation in selected kinematical range
176 AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
177 fTotalRate+=ratioupsilon;
178
179//------------------------------------------------------------------
180// Upsilon 2S generator
18bc0c8e 181 AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
103ac317 182// first step: generation in 4pi
183 genupsilonP->SetPtRange(0.,100.);
184 genupsilonP->SetYRange(-8.,8.);
185 genupsilonP->SetPhiRange(0.,360.);
aaa95f22 186 genupsilonP->SetForceDecay(fDecayModeResonance);
187 if (!gMC) genupsilonP->SetDecayer(fDecayer);
103ac317 188 genupsilonP->Init(); // generation in 4pi
189 ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
190 AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP));
191 AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
192// second step: generation in the kinematical range
193 genupsilonP->SetPtRange(ptMin, ptMax);
194 genupsilonP->SetYRange(yMin, yMax);
195 genupsilonP->SetPhiRange(phiMin, phiMax);
196 genupsilonP->Init(); // generation in selected kinematical range
197 AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
198 fTotalRate+=ratioupsilonP;
199
200//------------------------------------------------------------------
201// Upsilon 3S generator
18bc0c8e 202 AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
103ac317 203// first step: generation in 4pi
204 genupsilonPP->SetPtRange(0.,100.);
205 genupsilonPP->SetYRange(-8.,8.);
206 genupsilonPP->SetPhiRange(0.,360.);
aaa95f22 207 genupsilonPP->SetForceDecay(fDecayModeResonance);
208 if (!gMC) genupsilonPP->SetDecayer(fDecayer);
103ac317 209 genupsilonPP->Init(); // generation in 4pi
210 ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
211 AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP));
212 AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
213// second step: generation in selected kinematical range
214 genupsilonPP->SetPtRange(ptMin, ptMax);
215 genupsilonPP->SetYRange(yMin, yMax);
216 genupsilonPP->SetPhiRange(phiMin, phiMax);
217 genupsilonPP->Init(); // generation in selected kinematical range
218 AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
219 fTotalRate+=ratioupsilonPP;
b44c3901 220
221//------------------------------------------------------------------
222// Generator of charm
223
224 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4);
225 gencharm->SetMomentumRange(0,9999);
226 gencharm->SetForceDecay(kAll);
227 ratioccbar = sigmaccbar/sigmaReaction;
228 gencharm->Init();
229 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
230 fTotalRate+=ratioccbar;
231
103ac317 232//------------------------------------------------------------------
b44c3901 233// Generator of beauty
234
235 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5);
236 genbeauty->SetMomentumRange(0,9999);
237 genbeauty->SetForceDecay(kAll);
238 ratiobbbar = sigmabbbar/sigmaReaction;
239 genbeauty->Init();
240 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
241 fTotalRate+=ratiobbbar;
242
243//--------------------------t----------------------------------------
103ac317 244// Pythia generator
245 AliGenPythia *pythia = new AliGenPythia(1);
246 pythia->SetProcess(kPyMbMSEL1);
247 pythia->SetStrucFunc(kCTEQ5L);
248 pythia->SetEnergyCMS(14000.);
aaa95f22 249 AliInfo(Form("\n\npythia uses the decay mode %d",fDecayModePythia));
250 pythia->SetForceDecay(fDecayModePythia);
103ac317 251 pythia->SetPtRange(0.,100.);
252 pythia->SetYRange(-8.,8.);
253 pythia->SetPhiRange(0.,360.);
f58d19bb 254 pythia->SetPtHard(2.76,-1.0);
b44c3901 255 pythia->SwitchHFOff();
103ac317 256 pythia->Init();
257 AddGenerator(pythia,"Pythia",1);
258 fTotalRate+=1.;
259
aaa95f22 260 TIter next(fEntries);
261 AliGenCocktailEntry *entry;
262 if (fStack) {
263 while((entry = (AliGenCocktailEntry*)next())) {
264 entry->Generator()->SetStack(fStack);
265 }
266 }
103ac317 267}
268
269//_________________________________________________________________________
270void AliGenMUONCocktailpp::Generate()
271{
272// Generate event
273 TIter next(fEntries);
274 AliGenCocktailEntry *entry = 0;
275 AliGenCocktailEntry *preventry = 0;
276 AliGenerator* gen = 0;
277
278 TObjArray *partArray = gAlice->GetMCApp()->Particles();
279
280// Generate the vertex position used by all generators
281 if(fVertexSmear == kPerEvent) Vertex();
282
283// Loop on primordialTrigger:
284// minimum muon multiplicity above a pt cut in a theta acceptance region
285
286 Bool_t primordialTrigger = kFALSE;
287 while(!primordialTrigger) {
288 //Reseting stack
289 AliRunLoader * runloader = gAlice->GetRunLoader();
290 if (runloader)
291 if (runloader->Stack())
292 runloader->Stack()->Reset();
293 // Loop over generators and generate events
294 Int_t igen = 0;
295 Int_t npart = 0;
296 const char* genName = 0;
297 while((entry = (AliGenCocktailEntry*)next())) {
298 gen = entry->Generator();
299 genName = entry->GetName();
300 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
301
302 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
303 gRandom->Poisson(entry->Rate());
304
305 if (npart > 0) {
306 igen++;
307 if (igen == 1) entry->SetFirst(0);
308 else entry->SetFirst((partArray->GetEntriesFast())+1);
309
310 gen->SetNumberParticles(npart);
311 gen->Generate();
312 entry->SetLast(partArray->GetEntriesFast());
313 preventry = entry;
314 }
315 }
316 next.Reset();
317
318// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
319// in the muon spectrometer acceptance
320 Int_t iPart;
321 fNGenerated++;
b44c3901 322 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22 323 for(iPart=0; iPart<maxPart; iPart++){
103ac317 324
aaa95f22 325 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
326 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
327 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
328 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
329 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
330 (part->Pt()>fMuonPtCut)) {
331 numberOfMuons++;
103ac317 332 }
aaa95f22 333 }
334 }
103ac317 335 if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
336 }
337 fNSucceded++;
338
aaa95f22 339// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 340 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
341}
342
343