From Ivana: generalize geom_gentle_muon.C macro thus removing the need for additional...
[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),
9ff768ee 62 fNGenerated(0),
63 fSigmaJPsi(49.44e-6),
64 fSigmaPsiP(7.67e-6),
65 fSigmaUpsilon(0.989e-6),
66 fSigmaUpsilonP(0.502e-6),
67 fSigmaUpsilonPP(0.228e-6),
68 fSigmaCCbar(11.2e-3),
69 fSigmaBBbar(0.51e-3),
70 fSigmaSilent(kFALSE)
103ac317 71{
72// Constructor
103ac317 73}
93a2041b 74
103ac317 75//_________________________________________________________________________
76AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
77// Destructor
78{
79
80}
81
82//_________________________________________________________________________
8058ead1 83void AliGenMUONCocktailpp::CreateCocktail()
103ac317 84{
85// MinBias NN cross section @ pp 14 TeV -PR Vol II p:473
86 Double_t sigmaReaction = 0.070;
87
88// These limits are only used for renormalization of quarkonia cross section
89// Pythia events are generated in 4pi
90 Double_t ptMin = fPtMin;
91 Double_t ptMax = fPtMax;
92 Double_t yMin = fYMin;;
93 Double_t yMax = fYMax;;
94 Double_t phiMin = fPhiMin*180./TMath::Pi();
95 Double_t phiMax = fPhiMax*180./TMath::Pi();
96 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));
97
98// Ratios with respect to the reaction cross-section in the
99// kinematics limit of the MUONCocktail
100 Double_t ratiojpsi;
101 Double_t ratiopsiP;
102 Double_t ratioupsilon;
103 Double_t ratioupsilonP;
104 Double_t ratioupsilonPP;
b44c3901 105 Double_t ratioccbar;
106 Double_t ratiobbbar;
107
18bc0c8e 108// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
109// corrected from feed down of higher resonances
103ac317 110
9ff768ee 111 Double_t sigmajpsi = fSigmaJPsi;
112 Double_t sigmapsiP = fSigmaPsiP;
113 Double_t sigmaupsilon = fSigmaUpsilon;
114 Double_t sigmaupsilonP = fSigmaUpsilonP;
115 Double_t sigmaupsilonPP = fSigmaUpsilonPP;
116 Double_t sigmaccbar = fSigmaCCbar;
117 Double_t sigmabbbar = fSigmaBBbar;
3285b419 118
aaa95f22 119 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
120
18bc0c8e 121// Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
122//----------------------------------------------------------------------
103ac317 123// Jpsi generator
18bc0c8e 124 AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp", "Jpsi");
103ac317 125// first step: generation in 4pi
126 genjpsi->SetPtRange(0.,100.);
127 genjpsi->SetYRange(-8.,8.);
128 genjpsi->SetPhiRange(0.,360.);
aaa95f22 129 genjpsi->SetForceDecay(fDecayModeResonance);
130 if (!gMC) genjpsi->SetDecayer(fDecayer);
131
103ac317 132 genjpsi->Init(); // generation in 4pi
133 ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
9ff768ee 134 if (!fSigmaSilent) {
3285b419 135 AliInfo(Form("jpsi prod. cross-section in pp(14 TeV) %5.3g b",sigmajpsi));
136 AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
137 }
103ac317 138// second step: generation in selected kinematical range
139 genjpsi->SetPtRange(ptMin, ptMax);
140 genjpsi->SetYRange(yMin, yMax);
141 genjpsi->SetPhiRange(phiMin, phiMax);
142 genjpsi->Init(); // generation in selected kinematic range
143 AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
103ac317 144//------------------------------------------------------------------
145// Psi prime generator
18bc0c8e 146 AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
103ac317 147// first step: generation in 4pi
148 genpsiP->SetPtRange(0.,100.);
149 genpsiP->SetYRange(-8.,8.);
150 genpsiP->SetPhiRange(0.,360.);
aaa95f22 151 genpsiP->SetForceDecay(fDecayModeResonance);
152 if (!gMC) genpsiP->SetDecayer(fDecayer);
153
103ac317 154 genpsiP->Init(); // generation in 4pi
155 ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 156 if (!fSigmaSilent) {
3285b419 157 AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP));
158 AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
159 }
103ac317 160// second step: generation in selected kinematical range
161 genpsiP->SetPtRange(ptMin, ptMax);
162 genpsiP->SetYRange(yMin, yMax);
163 genpsiP->SetPhiRange(phiMin, phiMax);
164 genpsiP->Init(); // generation in selected kinematic range
165 AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
103ac317 166//------------------------------------------------------------------
167// Upsilon 1S generator
18bc0c8e 168 AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
103ac317 169// first step: generation in 4pi
170 genupsilon->SetPtRange(0.,100.);
171 genupsilon->SetYRange(-8.,8.);
172 genupsilon->SetPhiRange(0.,360.);
aaa95f22 173 genupsilon->SetForceDecay(fDecayModeResonance);
174 if (!gMC) genupsilon->SetDecayer(fDecayer);
103ac317 175 genupsilon->Init(); // generation in 4pi
176 ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 177 if (!fSigmaSilent) {
3285b419 178 AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon));
179 AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
180 }
103ac317 181// second step: generation in selected kinematical range
182 genupsilon->SetPtRange(ptMin, ptMax);
183 genupsilon->SetYRange(yMin, yMax);
184 genupsilon->SetPhiRange(phiMin, phiMax);
185 genupsilon->Init(); // generation in selected kinematical range
186 AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
103ac317 187//------------------------------------------------------------------
188// Upsilon 2S generator
18bc0c8e 189 AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
103ac317 190// first step: generation in 4pi
191 genupsilonP->SetPtRange(0.,100.);
192 genupsilonP->SetYRange(-8.,8.);
193 genupsilonP->SetPhiRange(0.,360.);
aaa95f22 194 genupsilonP->SetForceDecay(fDecayModeResonance);
195 if (!gMC) genupsilonP->SetDecayer(fDecayer);
103ac317 196 genupsilonP->Init(); // generation in 4pi
197 ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 198 if (!fSigmaSilent) {
3285b419 199 AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP));
200 AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
201 }
103ac317 202// second step: generation in the kinematical range
203 genupsilonP->SetPtRange(ptMin, ptMax);
204 genupsilonP->SetYRange(yMin, yMax);
205 genupsilonP->SetPhiRange(phiMin, phiMax);
206 genupsilonP->Init(); // generation in selected kinematical range
207 AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
103ac317 208
209//------------------------------------------------------------------
210// Upsilon 3S generator
18bc0c8e 211 AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
103ac317 212// first step: generation in 4pi
213 genupsilonPP->SetPtRange(0.,100.);
214 genupsilonPP->SetYRange(-8.,8.);
215 genupsilonPP->SetPhiRange(0.,360.);
aaa95f22 216 genupsilonPP->SetForceDecay(fDecayModeResonance);
217 if (!gMC) genupsilonPP->SetDecayer(fDecayer);
103ac317 218 genupsilonPP->Init(); // generation in 4pi
219 ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 220 if (!fSigmaSilent) {
3285b419 221 AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP));
222 AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
223 }
103ac317 224// second step: generation in selected kinematical range
225 genupsilonPP->SetPtRange(ptMin, ptMax);
226 genupsilonPP->SetYRange(yMin, yMax);
227 genupsilonPP->SetPhiRange(phiMin, phiMax);
228 genupsilonPP->Init(); // generation in selected kinematical range
229 AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
b44c3901 230
231//------------------------------------------------------------------
232// Generator of charm
233
234 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4);
8058ead1 235 gencharm->SetMomentumRange(0,9999);
236 gencharm->SetForceDecay(kAll);
237 ratioccbar = sigmaccbar/sigmaReaction;
238 if (!gMC) gencharm->SetDecayer(fDecayer);
239 gencharm->Init();
9ff768ee 240 if (!fSigmaSilent) {
3285b419 241 AliInfo(Form("c-cbar prod. cross-section in pp(14 TeV) %5.3g b",sigmaccbar));
242 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
243 }
8058ead1 244 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
b44c3901 245
103ac317 246//------------------------------------------------------------------
b44c3901 247// Generator of beauty
248
8058ead1 249 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5);
250 genbeauty->SetMomentumRange(0,9999);
251 genbeauty->SetForceDecay(kAll);
252 ratiobbbar = sigmabbbar/sigmaReaction;
253 if (!gMC) genbeauty->SetDecayer(fDecayer);
254 genbeauty->Init();
9ff768ee 255 if (!fSigmaSilent) {
3285b419 256 AliInfo(Form("b-bbar prod. cross-section in pp(14 TeV) %5.3g b",sigmabbbar));
257 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
258 }
8058ead1 259 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
b44c3901 260
8058ead1 261//-------------------------------------------------------------------
103ac317 262// Pythia generator
8058ead1 263//
264// This has to go into the Config.C
265//
266// AliGenPythia *pythia = new AliGenPythia(1);
267// pythia->SetProcess(kPyMbMSEL1);
268// pythia->SetStrucFunc(kCTEQ5L);
269// pythia->SetEnergyCMS(14000.);
270// AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
271// Decay_t dt = gener->GetDecayModePythia();
272// pythia->SetForceDecay(dt);
273// pythia->SetPtRange(0.,100.);
274// pythia->SetYRange(-8.,8.);
275// pythia->SetPhiRange(0.,360.);
276// pythia->SetPtHard(2.76,-1.0);
277// pythia->SwitchHFOff();
278// pythia->Init();
279// AddGenerator(pythia,"Pythia",1);
280}
103ac317 281
8058ead1 282void AliGenMUONCocktailpp::Init()
283{
284 //
285 // Initialisation
aaa95f22 286 TIter next(fEntries);
287 AliGenCocktailEntry *entry;
288 if (fStack) {
289 while((entry = (AliGenCocktailEntry*)next())) {
290 entry->Generator()->SetStack(fStack);
291 }
292 }
103ac317 293}
294
295//_________________________________________________________________________
296void AliGenMUONCocktailpp::Generate()
297{
298// Generate event
299 TIter next(fEntries);
300 AliGenCocktailEntry *entry = 0;
301 AliGenCocktailEntry *preventry = 0;
302 AliGenerator* gen = 0;
303
d94af0c1 304 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
103ac317 305
306// Generate the vertex position used by all generators
307 if(fVertexSmear == kPerEvent) Vertex();
308
309// Loop on primordialTrigger:
310// minimum muon multiplicity above a pt cut in a theta acceptance region
311
312 Bool_t primordialTrigger = kFALSE;
313 while(!primordialTrigger) {
314 //Reseting stack
33c3c91a 315 AliRunLoader * runloader = AliRunLoader::Instance();
103ac317 316 if (runloader)
317 if (runloader->Stack())
34da8494 318 runloader->Stack()->Clean();
103ac317 319 // Loop over generators and generate events
320 Int_t igen = 0;
321 Int_t npart = 0;
322 const char* genName = 0;
323 while((entry = (AliGenCocktailEntry*)next())) {
324 gen = entry->Generator();
325 genName = entry->GetName();
326 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
327
328 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
329 gRandom->Poisson(entry->Rate());
330
331 if (npart > 0) {
332 igen++;
333 if (igen == 1) entry->SetFirst(0);
334 else entry->SetFirst((partArray->GetEntriesFast())+1);
335
336 gen->SetNumberParticles(npart);
337 gen->Generate();
338 entry->SetLast(partArray->GetEntriesFast());
339 preventry = entry;
340 }
341 }
342 next.Reset();
343
344// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
345// in the muon spectrometer acceptance
346 Int_t iPart;
347 fNGenerated++;
b44c3901 348 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22 349 for(iPart=0; iPart<maxPart; iPart++){
103ac317 350
aaa95f22 351 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
352 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
353 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
354 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
355 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
f81a4aca 356 (part->Pt()>fMuonPtCut) &&
357 (part->P()>fMuonPCut)) {
aaa95f22 358 numberOfMuons++;
103ac317 359 }
aaa95f22 360 }
361 }
103ac317 362 if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
363 }
364 fNSucceded++;
365
aaa95f22 366// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 367 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
368}
369
370