1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
25 // July 07:added heavy quark production from AliGenCorrHF and heavy quark
26 // production switched off in Pythia
27 // Aug. 07: added trigger cut on total momentum
29 #include <TObjArray.h>
30 #include <TParticle.h>
32 #include <TVirtualMC.h>
34 #include "AliGenCocktailEventHeader.h"
36 #include "AliGenCocktailEntry.h"
37 #include "AliGenMUONCocktailpp.h"
38 #include "AliGenMUONlib.h"
39 #include "AliGenParam.h"
43 #include "AliDecayer.h"
45 #include "AliGenCorrHF.h"
47 ClassImp(AliGenMUONCocktailpp)
49 //________________________________________________________________________
50 AliGenMUONCocktailpp::AliGenMUONCocktailpp()
53 fDecayModeResonance(kAll),
54 fDecayModePythia(kAll),
59 fMuonThetaMaxCut(180.),
60 fMuonOriginCut(-999.),
67 //_________________________________________________________________________
68 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
74 //_________________________________________________________________________
75 void AliGenMUONCocktailpp::CreateCocktail()
77 // MinBias NN cross section @ pp 14 TeV -PR Vol II p:473
78 Double_t sigmaReaction = 0.070;
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));
90 // Ratios with respect to the reaction cross-section in the
91 // kinematics limit of the MUONCocktail
94 Double_t ratioupsilon;
95 Double_t ratioupsilonP;
96 Double_t ratioupsilonPP;
100 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
101 // corrected from feed down of higher resonances
103 Double_t sigmajpsi = 49.44e-6;
104 Double_t sigmapsiP = 7.67e-6;
105 Double_t sigmaupsilon = 0.989e-6;
106 Double_t sigmaupsilonP = 0.502e-6;
107 Double_t sigmaupsilonPP = 0.228e-6;
108 Double_t sigmaccbar = 11.2e-3;
109 Double_t sigmabbbar = 0.51e-3;
111 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
113 // Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
114 //----------------------------------------------------------------------
116 AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp", "Jpsi");
117 // first step: generation in 4pi
118 genjpsi->SetPtRange(0.,100.);
119 genjpsi->SetYRange(-8.,8.);
120 genjpsi->SetPhiRange(0.,360.);
121 genjpsi->SetForceDecay(fDecayModeResonance);
122 if (!gMC) genjpsi->SetDecayer(fDecayer);
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 //------------------------------------------------------------------
135 // Psi prime generator
136 AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
137 // first step: generation in 4pi
138 genpsiP->SetPtRange(0.,100.);
139 genpsiP->SetYRange(-8.,8.);
140 genpsiP->SetPhiRange(0.,360.);
141 genpsiP->SetForceDecay(fDecayModeResonance);
142 if (!gMC) genpsiP->SetDecayer(fDecayer);
144 genpsiP->Init(); // generation in 4pi
145 ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
146 AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP));
147 AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
148 // second step: generation in selected kinematical range
149 genpsiP->SetPtRange(ptMin, ptMax);
150 genpsiP->SetYRange(yMin, yMax);
151 genpsiP->SetPhiRange(phiMin, phiMax);
152 genpsiP->Init(); // generation in selected kinematic range
153 AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
154 //------------------------------------------------------------------
155 // Upsilon 1S generator
156 AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
157 // first step: generation in 4pi
158 genupsilon->SetPtRange(0.,100.);
159 genupsilon->SetYRange(-8.,8.);
160 genupsilon->SetPhiRange(0.,360.);
161 genupsilon->SetForceDecay(fDecayModeResonance);
162 if (!gMC) genupsilon->SetDecayer(fDecayer);
163 genupsilon->Init(); // generation in 4pi
164 ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
165 AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon));
166 AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
167 // second step: generation in selected kinematical range
168 genupsilon->SetPtRange(ptMin, ptMax);
169 genupsilon->SetYRange(yMin, yMax);
170 genupsilon->SetPhiRange(phiMin, phiMax);
171 genupsilon->Init(); // generation in selected kinematical range
172 AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
173 //------------------------------------------------------------------
174 // Upsilon 2S generator
175 AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
176 // first step: generation in 4pi
177 genupsilonP->SetPtRange(0.,100.);
178 genupsilonP->SetYRange(-8.,8.);
179 genupsilonP->SetPhiRange(0.,360.);
180 genupsilonP->SetForceDecay(fDecayModeResonance);
181 if (!gMC) genupsilonP->SetDecayer(fDecayer);
182 genupsilonP->Init(); // generation in 4pi
183 ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
184 AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP));
185 AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
186 // second step: generation in the kinematical range
187 genupsilonP->SetPtRange(ptMin, ptMax);
188 genupsilonP->SetYRange(yMin, yMax);
189 genupsilonP->SetPhiRange(phiMin, phiMax);
190 genupsilonP->Init(); // generation in selected kinematical range
191 AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
193 //------------------------------------------------------------------
194 // Upsilon 3S generator
195 AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
196 // first step: generation in 4pi
197 genupsilonPP->SetPtRange(0.,100.);
198 genupsilonPP->SetYRange(-8.,8.);
199 genupsilonPP->SetPhiRange(0.,360.);
200 genupsilonPP->SetForceDecay(fDecayModeResonance);
201 if (!gMC) genupsilonPP->SetDecayer(fDecayer);
202 genupsilonPP->Init(); // generation in 4pi
203 ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
204 AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP));
205 AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
206 // second step: generation in selected kinematical range
207 genupsilonPP->SetPtRange(ptMin, ptMax);
208 genupsilonPP->SetYRange(yMin, yMax);
209 genupsilonPP->SetPhiRange(phiMin, phiMax);
210 genupsilonPP->Init(); // generation in selected kinematical range
211 AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
213 //------------------------------------------------------------------
214 // Generator of charm
216 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4);
217 gencharm->SetMomentumRange(0,9999);
218 gencharm->SetForceDecay(kAll);
219 ratioccbar = sigmaccbar/sigmaReaction;
220 if (!gMC) gencharm->SetDecayer(fDecayer);
222 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
224 //------------------------------------------------------------------
225 // Generator of beauty
227 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5);
228 genbeauty->SetMomentumRange(0,9999);
229 genbeauty->SetForceDecay(kAll);
230 ratiobbbar = sigmabbbar/sigmaReaction;
231 if (!gMC) genbeauty->SetDecayer(fDecayer);
233 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
235 //-------------------------------------------------------------------
238 // This has to go into the Config.C
240 // AliGenPythia *pythia = new AliGenPythia(1);
241 // pythia->SetProcess(kPyMbMSEL1);
242 // pythia->SetStrucFunc(kCTEQ5L);
243 // pythia->SetEnergyCMS(14000.);
244 // AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
245 // Decay_t dt = gener->GetDecayModePythia();
246 // pythia->SetForceDecay(dt);
247 // pythia->SetPtRange(0.,100.);
248 // pythia->SetYRange(-8.,8.);
249 // pythia->SetPhiRange(0.,360.);
250 // pythia->SetPtHard(2.76,-1.0);
251 // pythia->SwitchHFOff();
253 // AddGenerator(pythia,"Pythia",1);
256 void AliGenMUONCocktailpp::Init()
260 TIter next(fEntries);
261 AliGenCocktailEntry *entry;
263 while((entry = (AliGenCocktailEntry*)next())) {
264 entry->Generator()->SetStack(fStack);
269 //_________________________________________________________________________
270 void AliGenMUONCocktailpp::Generate()
273 TIter next(fEntries);
274 AliGenCocktailEntry *entry = 0;
275 AliGenCocktailEntry *preventry = 0;
276 AliGenerator* gen = 0;
278 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
280 // Generate the vertex position used by all generators
281 if(fVertexSmear == kPerEvent) Vertex();
283 // Loop on primordialTrigger:
284 // minimum muon multiplicity above a pt cut in a theta acceptance region
286 Bool_t primordialTrigger = kFALSE;
287 while(!primordialTrigger) {
289 AliRunLoader * runloader = gAlice->GetRunLoader();
291 if (runloader->Stack())
292 runloader->Stack()->Clean();
293 // Loop over generators and generate events
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));
302 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
303 gRandom->Poisson(entry->Rate());
307 if (igen == 1) entry->SetFirst(0);
308 else entry->SetFirst((partArray->GetEntriesFast())+1);
310 gen->SetNumberParticles(npart);
312 entry->SetLast(partArray->GetEntriesFast());
318 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
319 // in the muon spectrometer acceptance
322 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
323 for(iPart=0; iPart<maxPart; iPart++){
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 (part->P()>fMuonPCut)) {
336 if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
340 // AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
341 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));