]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EVGEN/AliGenMUONCocktailpp.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
reduce histogram bin size for delta eta close to 0 from 0.05 to 0.01
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
... / ...
CommitLineData
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// Class to create the coktail for physics with muons for pp collisions
20// using the The followoing sources:
21// jpsi, psiP, upsilon, upsilonP, upsilonPP, open charm and open beauty
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
28// 2009: added possibility to hide x-sections (B. Vulpescu)
29// 2009: added possibility to have the cocktail (fast generator and param.)
30// for pp @ 10 TeV or pp @ 14 TeV (N. Bastid)
31//-----------------
32// 2009: added polarization (L. Bianchi)
33//------------------
34// 11/2009: added chi_c1 & chi_c2 (P.Crochet & N.Bastid).
35// Cross-sections for charmonia are now directly taken from the Yellow Report
36// (hep-ph/0311048) Tab.9, page 19. See below for details w.r.t. beam energy.
37// usage: see example of Config in $ALICE_ROOT/prod/LHC09a10/Config.C
38//------------------------
39// 04/2010:
40// - CMS energy passed via parameter
41// i.e. gener->SetCMSEnergy(AliGenMUONCocktailpp::kCMS07TeV) in Config.C
42// - resonances now added to the cocktail via AddReso2Generator
43// - cleaning
44// B.Vulpescu & P.Crochet
45//-----------------------
46// 10/2011:
47// - added the cocktail for p-Pb & Pb-p @ 8.8 TeV with 4 centrality bins and
48// for Pb-Pb @ 2.76 TeV with 11 centrality bins. Bins should be defined also
49// in the Config.C with one AliGenMUONCocktailpp per bin. These generators
50// included in a AliGenCocktail together with an event generator (e.g. Hijing)
51// providing the underlying event and collision centrality. The bin number n
52// passed via AliGenMUONCocktailpp::SetCentralityBin(n).
53// See details in my presentation at the PWG3-Muon meeting (05.10.2011):
54// https://indico.cern.ch/conferenceDisplay.py?confId=157367
55// - simplifications and bug fix in CreateCocktail()
56// S. Grigoryan
57//-----------------------
58// 06/2012:
59// - added the cocktail for p-Pb & Pb-p @ 2.76, 4.4 & 5.03 TeV with 4 centrality
60// bins, using the EPS09-LO shadowing computed for 5.03 TeV. Energies are set by
61// AliGenMUONCocktailpp::SetCMSEnergy(int CMSEnergyCode), CMSEnergy codes are
62// defined in AliGenMUONCocktailpp.h.
63// - added functions to scale x-section of JPsi, Charmonia, Bottomonia, CCbar & BBbar
64// in Config.C, e.g. AliGenMUONCocktailpp::ScaleJPsi(2.5), to manage the statistics.
65// S. Grigoryan
66
67#include <TObjArray.h>
68#include <TParticle.h>
69#include <TF1.h>
70#include <TVirtualMC.h>
71#include "AliGenCocktailEventHeader.h"
72
73#include "AliGenCocktailEntry.h"
74#include "AliGenMUONCocktailpp.h"
75#include "AliGenMUONlib.h"
76#include "AliGenParam.h"
77#include "AliMC.h"
78#include "AliRun.h"
79#include "AliStack.h"
80#include "AliDecayer.h"
81#include "AliLog.h"
82#include "AliGenCorrHF.h"
83#include "AliDecayerPolarized.h"
84
85ClassImp(AliGenMUONCocktailpp)
86
87//________________________________________________________________________
88AliGenMUONCocktailpp::AliGenMUONCocktailpp()
89 :AliGenCocktail(),
90 fDecayer(0),
91 fDecayModeResonance(kAll),
92 fDecayModePythia(kAll),
93 fMuonMultiplicity(0),
94 fMuonPtCut(0.),
95 fMuonPCut(0.),
96 fMuonThetaMinCut(0.),
97 fMuonThetaMaxCut(180.),
98 fMuonOriginCut(-999.),
99 fNSucceded(0),
100 fNGenerated(0),
101 fCentralityBin(0),
102 fScaleJPsi(1),
103 fScaleCharmonia(1),
104 fScaleBottomonia(1),
105 fScaleCCbar(1),
106 fScaleBBbar(1),
107
108 fJpsiPol(0),
109 fChic1Pol(0),
110 fChic2Pol(0),
111 fPsiPPol(0),
112 fUpsPol(0),
113 fUpsPPol(0),
114 fUpsPPPol(0),
115 fPolFrame(0),
116
117 fCMSEnergyTeV(0),
118 fCMSEnergyTeVArray(),
119 fSigmaReaction(0),
120 fSigmaReactionArray(),
121 fSigmaJPsi(0),
122 fSigmaJPsiArray(),
123 fSigmaChic1(0),
124 fSigmaChic1Array(),
125 fSigmaChic2(0),
126 fSigmaChic2Array(),
127 fSigmaPsiP(0),
128 fSigmaPsiPArray(),
129 fSigmaUpsilon(0),
130 fSigmaUpsilonArray(),
131 fSigmaUpsilonP(0),
132 fSigmaUpsilonPArray(),
133 fSigmaUpsilonPP(0),
134 fSigmaUpsilonPPArray(),
135 fSigmaCCbar(0),
136 fSigmaCCbarArray(),
137 fSigmaBBbar(0),
138 fSigmaBBbarArray(),
139
140 fSigmaSilent(kFALSE)
141{
142// Constructor
143
144// x-sections for pp @ 7 TeV:
145// -charmonia: 4pi integral of fit function for inclusive J/psi dsigma/dy LHC data
146// gives 60 mub; so sigma_prompt = 54 mub, while Ref = R.Vogt_arXiv:1003.3497 (Table 2)
147// gives 35 mub. Below we use sigma_direct from the Ref scaled by the factor 54/35.
148// -bottomonia: 4pi integral of fit function for inclusive Upsilon1S dsigma/dy LHC data
149// gives 0.56 mub, sigmas for 2S & 3S obtained using LHCb data for ratios 2S/1S & 3S/1S
150// -ccbar & bbbar: NLO pQCD computations - http://www-alice.gsi.de/ana/MNR/results.html
151 fCMSEnergyTeVArray[0] = 7.00;
152 fSigmaReactionArray[0] = 0.070;
153 fSigmaJPsiArray[0] = 33.6e-6;
154 fSigmaChic1Array[0] = 32.6e-6;
155 fSigmaChic2Array[0] = 53.8e-6;
156 fSigmaPsiPArray[0] = 7.6e-6;
157 fSigmaUpsilonArray[0] = 0.56e-6;
158 fSigmaUpsilonPArray[0] = 0.18e-6;
159 fSigmaUpsilonPPArray[0] = 0.08e-6;
160 fSigmaCCbarArray[0] = 6.91e-3;
161 fSigmaBBbarArray[0] = 0.232e-3;
162
163//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
164// scaled down according to ccbar and bbbar cross-sections
165 fCMSEnergyTeVArray[1] = 10.00;
166 fSigmaReactionArray[1] = 0.070;
167 fSigmaJPsiArray[1] = 26.06e-6;
168 fSigmaChic1Array[1] = 25.18e-6;
169 fSigmaChic2Array[1] = 41.58e-6;
170 fSigmaPsiPArray[1] = 5.88e-6;
171 fSigmaUpsilonArray[1] = 0.658e-6;
172 fSigmaUpsilonPArray[1] = 0.218e-6;
173 fSigmaUpsilonPPArray[1] = 0.122e-6;
174 fSigmaCCbarArray[1] = 8.9e-3;
175 fSigmaBBbarArray[1] = 0.33e-3;
176
177//x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
178// bottomonium from hep-ph/0311048 Tab.9, page 19 taken into account that
179// feed-down from chib is included
180 fCMSEnergyTeVArray[2] = 14.00;
181 fSigmaReactionArray[2] = 0.070;
182 fSigmaJPsiArray[2] = 32.9e-6;
183 fSigmaChic1Array[2] = 31.8e-6;
184 fSigmaChic2Array[2] = 52.5e-6;
185 fSigmaPsiPArray[2] = 7.43e-6;
186 fSigmaUpsilonArray[2] = 0.989e-6;
187 fSigmaUpsilonPArray[2] = 0.502e-6;
188 fSigmaUpsilonPPArray[2] = 0.228e-6;
189 fSigmaCCbarArray[2] = 11.2e-3;
190 fSigmaBBbarArray[2] = 0.445e-3;
191
192// x-sections for Min. Bias p-Pb & Pb-p @ 2.76 TeV: charmonia and bottomonia
193// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
194// and with Glauber scaling
195 fCMSEnergyTeVArray[3] = 2.00; // for 2.76 TeV
196 fSigmaReactionArray[3] = 2.10;
197 fSigmaJPsiArray[3] = 3.54e-3; // 208*0.507*33.6e-6
198 fSigmaChic1Array[3] = 3.44e-3;
199 fSigmaChic2Array[3] = 5.66e-3;
200 fSigmaPsiPArray[3] = 0.80e-3;
201 fSigmaUpsilonArray[3] = 0.045e-3; // 208*0.384*0.56e-6
202 fSigmaUpsilonPArray[3] = 0.014e-3;
203 fSigmaUpsilonPPArray[3] = 0.006e-3;
204 fSigmaCCbarArray[3] = 0.73; // 208*3.50e-3
205 fSigmaBBbarArray[3] = 0.019; // 208*0.089e-3
206
207 fCMSEnergyTeVArray[4] = -fCMSEnergyTeVArray[3];
208 fSigmaReactionArray[4] = fSigmaReactionArray[3];
209 fSigmaJPsiArray[4] = fSigmaJPsiArray[3];
210 fSigmaChic1Array[4] = fSigmaChic1Array[3];
211 fSigmaChic2Array[4] = fSigmaChic2Array[3];
212 fSigmaPsiPArray[4] = fSigmaPsiPArray[3];
213 fSigmaUpsilonArray[4] = fSigmaUpsilonArray[3];
214 fSigmaUpsilonPArray[4] = fSigmaUpsilonPArray[3];
215 fSigmaUpsilonPPArray[4] = fSigmaUpsilonPPArray[3];
216 fSigmaCCbarArray[4] = fSigmaCCbarArray[3];
217 fSigmaBBbarArray[4] = fSigmaBBbarArray[3];
218
219// x-sections for Min. Bias p-Pb & Pb-p @ 4.4 TeV: charmonia and bottomonia
220// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
221// and with Glauber scaling
222 fCMSEnergyTeVArray[5] = 4.00; // for 4.4 TeV
223 fSigmaReactionArray[5] = 2.10;
224 fSigmaJPsiArray[5] = 5.00e-3; // 208*0.715*33.6e-6
225 fSigmaChic1Array[5] = 4.86e-3;
226 fSigmaChic2Array[5] = 7.99e-3;
227 fSigmaPsiPArray[5] = 1.12e-3;
228 fSigmaUpsilonArray[5] = 0.074e-3; // 208*0.629*0.56e-6
229 fSigmaUpsilonPArray[5] = 0.023e-3;
230 fSigmaUpsilonPPArray[5] = 0.010e-3;
231 fSigmaCCbarArray[5] = 1.03; // 208*4.94e-3
232 fSigmaBBbarArray[5] = 0.030; // 208*0.146e-3
233
234 fCMSEnergyTeVArray[6] = -fCMSEnergyTeVArray[5];
235 fSigmaReactionArray[6] = fSigmaReactionArray[5];
236 fSigmaJPsiArray[6] = fSigmaJPsiArray[5];
237 fSigmaChic1Array[6] = fSigmaChic1Array[5];
238 fSigmaChic2Array[6] = fSigmaChic2Array[5];
239 fSigmaPsiPArray[6] = fSigmaPsiPArray[5];
240 fSigmaUpsilonArray[6] = fSigmaUpsilonArray[5];
241 fSigmaUpsilonPArray[6] = fSigmaUpsilonPArray[5];
242 fSigmaUpsilonPPArray[6] = fSigmaUpsilonPPArray[5];
243 fSigmaCCbarArray[6] = fSigmaCCbarArray[5];
244 fSigmaBBbarArray[6] = fSigmaBBbarArray[5];
245
246// x-sections for Min. Bias p-Pb & Pb-p @ 5.03 TeV: charmonia and bottomonia
247// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
248// and with Glauber scaling
249 fCMSEnergyTeVArray[7] = 5.00; // for 5.03 TeV
250 fSigmaReactionArray[7] = 2.10;
251 fSigmaJPsiArray[7] = 5.50e-3; // 208*0.787*33.6e-6
252 fSigmaChic1Array[7] = 5.35e-3;
253 fSigmaChic2Array[7] = 8.79e-3;
254 fSigmaPsiPArray[7] = 1.23e-3;
255 fSigmaUpsilonArray[7] = 0.083e-3; // 208*0.716*0.56e-6
256 fSigmaUpsilonPArray[7] = 0.026e-3;
257 fSigmaUpsilonPPArray[7] = 0.011e-3;
258 fSigmaCCbarArray[7] = 1.13; // 208*5.44e-3
259 fSigmaBBbarArray[7] = 0.035; // 208*0.166e-3
260
261 fCMSEnergyTeVArray[8] = -fCMSEnergyTeVArray[7];
262 fSigmaReactionArray[8] = fSigmaReactionArray[7];
263 fSigmaJPsiArray[8] = fSigmaJPsiArray[7];
264 fSigmaChic1Array[8] = fSigmaChic1Array[7];
265 fSigmaChic2Array[8] = fSigmaChic2Array[7];
266 fSigmaPsiPArray[8] = fSigmaPsiPArray[7];
267 fSigmaUpsilonArray[8] = fSigmaUpsilonArray[7];
268 fSigmaUpsilonPArray[8] = fSigmaUpsilonPArray[7];
269 fSigmaUpsilonPPArray[8] = fSigmaUpsilonPPArray[7];
270 fSigmaCCbarArray[8] = fSigmaCCbarArray[7];
271 fSigmaBBbarArray[8] = fSigmaBBbarArray[7];
272
273// x-sections for Min. Bias p-Pb & Pb-p @ 8.8 TeV: charmonia and bottomonia
274// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
275// and with Glauber scaling
276 fCMSEnergyTeVArray[9] = 9.00; // for 8.8 TeV
277 fSigmaReactionArray[9] = 2.10;
278 fSigmaJPsiArray[9] = 8.19e-3; // 208*1.172*33.6e-6
279 fSigmaChic1Array[9] = 7.95e-3;
280 fSigmaChic2Array[9] = 13.1e-3;
281 fSigmaPsiPArray[9] = 1.85e-3;
282 fSigmaUpsilonArray[9] = 0.146e-3; // 208*1.25*0.56e-6
283 fSigmaUpsilonPArray[9] = 0.047e-3;
284 fSigmaUpsilonPPArray[9] = 0.021e-3;
285 fSigmaCCbarArray[9] = 1.68; // 208*8.1e-3
286 fSigmaBBbarArray[9] = 0.061; // 208*0.29e-3
287
288 fCMSEnergyTeVArray[10] = -fCMSEnergyTeVArray[9];
289 fSigmaReactionArray[10] = fSigmaReactionArray[9];
290 fSigmaJPsiArray[10] = fSigmaJPsiArray[9];
291 fSigmaChic1Array[10] = fSigmaChic1Array[9];
292 fSigmaChic2Array[10] = fSigmaChic2Array[9];
293 fSigmaPsiPArray[10] = fSigmaPsiPArray[9];
294 fSigmaUpsilonArray[10] = fSigmaUpsilonArray[9];
295 fSigmaUpsilonPArray[10] = fSigmaUpsilonPArray[9];
296 fSigmaUpsilonPPArray[10] = fSigmaUpsilonPPArray[9];
297 fSigmaCCbarArray[10] = fSigmaCCbarArray[9];
298 fSigmaBBbarArray[10] = fSigmaBBbarArray[9];
299
300// x-sections for Min. Bias Pb-Pb @ 2.76 TeV: charmonia and bottomonia
301// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
302// and with Glauber scaling
303 fCMSEnergyTeVArray[11] = 3.00; // for 2.76 TeV
304 fSigmaReactionArray[11] = 7.65;
305 fSigmaJPsiArray[11] = 0.737; // 208*208*0.507*33.6e-6
306 fSigmaChic1Array[11] = 0.715;
307 fSigmaChic2Array[11] = 1.179;
308 fSigmaPsiPArray[11] = 0.166;
309 fSigmaUpsilonArray[11] = 0.0093; // 208*208*0.384*0.56e-6
310 fSigmaUpsilonPArray[11] = 0.0030;
311 fSigmaUpsilonPPArray[11] = 0.0013;
312 fSigmaCCbarArray[11] = 151.; // 208*208*3.50e-3
313 fSigmaBBbarArray[11] = 3.8; // 208*208*0.089e-3
314
315}
316
317//_________________________________________________________________________
318AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
319{
320// Destructor
321
322}
323
324//_________________________________________________________________________
325void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy)
326{
327// setter for CMSEnergy and corresponding cross-sections
328 fCMSEnergyTeV = fCMSEnergyTeVArray[cmsEnergy];
329 fSigmaReaction = fSigmaReactionArray[cmsEnergy];
330 fSigmaJPsi = fSigmaJPsiArray[cmsEnergy];
331 fSigmaChic1 = fSigmaChic1Array[cmsEnergy];
332 fSigmaChic2 = fSigmaChic2Array[cmsEnergy];
333 fSigmaPsiP = fSigmaPsiPArray[cmsEnergy];
334 fSigmaUpsilon = fSigmaUpsilonArray[cmsEnergy];
335 fSigmaUpsilonP = fSigmaUpsilonPArray[cmsEnergy];
336 fSigmaUpsilonPP = fSigmaUpsilonPPArray[cmsEnergy];
337 fSigmaCCbar = fSigmaCCbarArray[cmsEnergy];
338 fSigmaBBbar = fSigmaBBbarArray[cmsEnergy];
339}
340
341//_________________________________________________________________________
342void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
343 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
344// setter for resonances polarization
345 if (strcmp(PolFrame,"kColSop")==0){
346 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
347 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
348 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
349 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
350 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
351 fPolFrame = 0;
352 } else if (strcmp(PolFrame,"kHelicity")==0){
353 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
354 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
355 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
356 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
357 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
358 fPolFrame = 1;
359
360 } else {
361 AliInfo(Form("The polarization frame is not valid"));
362 AliInfo(Form("No polarization will be set"));
363 fJpsiPol=0.;
364 fPsiPPol=0.;
365 fUpsPol=0.;
366 fUpsPPol=0.;
367 fUpsPPPol=0.;
368 }
369}
370
371//_________________________________________________________________________
372void AliGenMUONCocktailpp::CreateCocktail()
373{
374// create and add resonances and open HF to the coctail
375 Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
376
377 // For temporary use of p-Pb & Pb-p shadowing at 5.03 TeV for lower energies
378 if (cmsEnergy == 2 || cmsEnergy == 4) cmsEnergy = 5;
379 if (cmsEnergy == -2 || cmsEnergy == -4) cmsEnergy = -5;
380
381// These limits are only used for renormalization of quarkonia cross section
382// Pythia events are generated in 4pi
383 Double_t ptMin = fPtMin;
384 Double_t ptMax = fPtMax;
385 Double_t yMin = fYMin;;
386 Double_t yMax = fYMax;;
387 Double_t phiMin = fPhiMin*180./TMath::Pi();
388 Double_t phiMax = fPhiMax*180./TMath::Pi();
389 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));
390
391// Cross sections in barns
392
393 Double_t sigmajpsi = fSigmaJPsi;
394 Double_t sigmachic1 = fSigmaChic1;
395 Double_t sigmachic2 = fSigmaChic2;
396 Double_t sigmapsiP = fSigmaPsiP;
397 Double_t sigmaupsilon = fSigmaUpsilon;
398 Double_t sigmaupsilonP = fSigmaUpsilonP;
399 Double_t sigmaupsilonPP = fSigmaUpsilonPP;
400 Double_t sigmaccbar = fSigmaCCbar;
401 Double_t sigmabbbar = fSigmaBBbar;
402
403// Cross sections corrected with the BR in mu+mu-
404// (only in case of use of AliDecayerPolarized)
405
406 if(TMath::Abs(fJpsiPol) > 1.e-30) {sigmajpsi = fSigmaJPsi*0.0593;}
407 if(TMath::Abs(fChic1Pol) > 1.e-30) {sigmachic1 = fSigmaChic1*0.;} // tb consistent
408 if(TMath::Abs(fChic2Pol) > 1.e-30) {sigmachic2 = fSigmaChic2*0.;} // tb consistent
409 if(TMath::Abs(fPsiPPol) > 1.e-30) {sigmapsiP = fSigmaPsiP*0.0075;}
410 if(TMath::Abs(fUpsPol) > 1.e-30) {sigmaupsilon = fSigmaUpsilon*0.0248;}
411 if(TMath::Abs(fUpsPPol) > 1.e-30) {sigmaupsilonP = fSigmaUpsilonP*0.0193;}
412 if(TMath::Abs(fUpsPPPol) > 1.e-30) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
413
414// Cross sections scaled to manage the statistics
415
416 sigmajpsi *= fScaleJPsi*fScaleCharmonia;
417 sigmachic1 *= fScaleCharmonia;
418 sigmachic2 *= fScaleCharmonia;
419 sigmapsiP *= fScaleCharmonia;
420 sigmaupsilon *= fScaleBottomonia;
421 sigmaupsilonP *= fScaleBottomonia;
422 sigmaupsilonPP *= fScaleBottomonia;
423 sigmaccbar *= fScaleCCbar;
424 sigmabbbar *= fScaleBBbar;
425
426 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
427
428// Create and add resonances to the generator
429 AliGenParam * genjpsi=0;
430 AliGenParam * genchic1=0;
431 AliGenParam * genchic2=0;
432 AliGenParam * genpsiP=0;
433 AliGenParam * genupsilon=0;
434 AliGenParam * genupsilonP=0;
435 AliGenParam * genupsilonPP=0;
436
437 Char_t nameJpsi[10];
438 Char_t nameChic1[10];
439 Char_t nameChic2[10];
440 Char_t namePsiP[10];
441 Char_t nameUps[10];
442 Char_t nameUpsP[10];
443 Char_t nameUpsPP[10];
444
445 snprintf(nameJpsi,10, "Jpsi");
446 snprintf(nameChic1,10, "Chic1");
447 snprintf(nameChic2,10, "Chic2");
448 snprintf(namePsiP,10, "PsiP");
449 snprintf(nameUps,10, "Ups");
450 snprintf(nameUpsP,10, "UpsP");
451 snprintf(nameUpsPP,10, "UpsPP");
452
453 Char_t tname[40] = "";
454 if(cmsEnergy == 10) {snprintf(tname, 40, "CDF pp 10");
455 } else if (cmsEnergy == 14){snprintf(tname, 40, "CDF pp");
456 } else if (cmsEnergy == 7) {snprintf(tname, 40, "pp 7");
457 // } else if (cmsEnergy == 2) {snprintf(tname, 40, "pp 2.76");
458 } else if (cmsEnergy == 5) {snprintf(tname, 40, "pPb 5.03");
459 if (fCentralityBin > 0) snprintf(tname, 40, "pPb 5.03c%d",fCentralityBin);
460 } else if (cmsEnergy == -5){snprintf(tname, 40, "Pbp 5.03");
461 if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 5.03c%d",fCentralityBin);
462 } else if (cmsEnergy == 9) {snprintf(tname, 40, "pPb 8.8");
463 if (fCentralityBin > 0) snprintf(tname, 40, "pPb 8.8c%d",fCentralityBin);
464 } else if (cmsEnergy == -9){snprintf(tname, 40, "Pbp 8.8");
465 if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 8.8c%d",fCentralityBin);
466 } else if (cmsEnergy == 3) {snprintf(tname, 40, "PbPb 2.76");
467 if (fCentralityBin > 0) snprintf(tname, 40, "PbPb 2.76c%d",fCentralityBin);
468 } else {
469 AliError("Initialisation failed, wrong cmsEnergy");
470 return;
471 }
472 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, tname, "Jpsi");
473 genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, tname, "Chic1");
474 genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, tname, "Chic2");
475 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, tname, "PsiP");
476 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, tname, "Upsilon");
477 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, tname, "UpsilonP");
478 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, tname, "UpsilonPP");
479
480// Hard process yield per pA or AA collision for i-th centrality bin is R*r[i]*shad[i]
481// where R is the ratio of hard and geometrical x-sections, r[i] is the ratio of these
482// x-section fractions for given centrality and shad[i] is the shadowing factor (in 4pi).
483// The latter is assumed to be the same for HF-hadrons & quarkonia of the same flavour.
484 Int_t i = 0;
485 Double_t chard[20] = {0}; // charm & beauty shadowing factors are different
486 Double_t bhard[20] = {0};
487 chard[0] = 1; // 1st element for pp and min. bias (MB) collisions
488 bhard[0] = 1;
489
490// 4 centrality bins for p-Pb & Pb-p at 5.03 TeV: 0-20-40-60-100 %
491 if (cmsEnergy == 5 || cmsEnergy == -5) {
492 const Int_t n5 = 5; // 1st element for MB collisions
493 Double_t r5[n5] = {1, 1.936, 1.473, 0.914, 0.333}; // ratio of hard-over-geo fractions
494 Double_t cshad5[n5] = {0.806, 0.742, 0.796, 0.870, 0.955};// EPS09-LO shadowing factors
495 Double_t bshad5[n5] = {0.917, 0.889, 0.913, 0.944, 0.981};
496 for(i=0; i<n5; i++) {
497 chard[i] = cshad5[i]*r5[i];
498 bhard[i] = bshad5[i]*r5[i];
499 }
500 }
501
502// 4 centrality bins for p-Pb & Pb-p at 8.8 TeV: 0-20-40-60-100 %
503 if (cmsEnergy == 9 || cmsEnergy == -9) {
504 const Int_t n9 = 5; // 1st element for MB collisions
505 Double_t r9[n9] = {1, 1.936, 1.473, 0.914, 0.333}; // ratio of hard-over-geo fractions
506 Double_t cshad9[n9] = {0.785, 0.715, 0.775, 0.856, 0.951};// EKS98 shadowing factors
507 Double_t bshad9[n9] = {0.889, 0.853, 0.884, 0.926, 0.975};
508 for(i=0; i<n9; i++) {
509 chard[i] = cshad9[i]*r9[i];
510 bhard[i] = bshad9[i]*r9[i];
511 }
512 }
513
514// 11 centrality bins for Pb-Pb at 2.76 TeV: 0-5-10-20-30-40-50-60-70-80-90-100 %
515 if (cmsEnergy == 3) {
516 const Int_t n3 = 12; // 1st element for MB collisions
517 Double_t r3[n3] = {1, 4.661, 3.647, 2.551, 1.544, 0.887, 0.474,
518 0.235, 0.106, 0.044, 0.017, 0.007}; // ratio of hard-over-geo fractions
519 Double_t cshad3[n3] = {0.662, 0.622, 0.631, 0.650, 0.681, 0.718,
520 0.760, 0.805, 0.849, 0.888, 0.918, 0.944};// EKS98 shadowing factors
521 Double_t bshad3[n3] = {0.874, 0.856, 0.861, 0.869, 0.883, 0.898,
522 0.915, 0.932, 0.948, 0.962, 0.972, 0.981};
523 for(i=0; i<n3; i++) {
524 chard[i] = cshad3[i]*r3[i];
525 bhard[i] = bshad3[i]*r3[i];
526 }
527 }
528
529 AddReso2Generator(nameJpsi,genjpsi,chard[fCentralityBin]*sigmajpsi,fJpsiPol);
530 AddReso2Generator(nameChic1,genchic1,chard[fCentralityBin]*sigmachic1,fChic1Pol);
531 AddReso2Generator(nameChic2,genchic2,chard[fCentralityBin]*sigmachic2,fChic2Pol);
532 AddReso2Generator(namePsiP,genpsiP,chard[fCentralityBin]*sigmapsiP,fPsiPPol);
533
534 AddReso2Generator(nameUps,genupsilon,bhard[fCentralityBin]*sigmaupsilon,fUpsPol);
535 AddReso2Generator(nameUpsP,genupsilonP,bhard[fCentralityBin]*sigmaupsilonP,fUpsPPol);
536 AddReso2Generator(nameUpsPP,genupsilonPP,bhard[fCentralityBin]*sigmaupsilonPP,fUpsPPPol);
537
538//------------------------------------------------------------------
539// Generator of charm
540 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
541 gencharm->SetMomentumRange(0,9999);
542 gencharm->SetForceDecay(kAll);
543 Double_t ratioccbar = chard[fCentralityBin]*sigmaccbar/fSigmaReaction;
544 if (!gMC) gencharm->SetDecayer(fDecayer);
545 gencharm->Init();
546 if (!fSigmaSilent) {
547 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
548 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
549 }
550 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
551//------------------------------------------------------------------
552// Generator of beauty
553 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
554 genbeauty->SetMomentumRange(0,9999);
555 genbeauty->SetForceDecay(kAll);
556 Double_t ratiobbbar = bhard[fCentralityBin]*sigmabbbar/fSigmaReaction;
557 if (!gMC) genbeauty->SetDecayer(fDecayer);
558 genbeauty->Init();
559 if (!fSigmaSilent) {
560 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
561 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
562 }
563 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
564
565//-------------------------------------------------------------------
566// Pythia generator
567//
568// This has to go into the Config.C
569//
570// AliGenPythia *pythia = new AliGenPythia(1);
571// pythia->SetProcess(kPyMbMSEL1);
572// pythia->SetStrucFunc(kCTEQ5L);
573// pythia->SetEnergyCMS(14000.);
574// AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
575// Decay_t dt = gener->GetDecayModePythia();
576// pythia->SetForceDecay(dt);
577// pythia->SetPtRange(0.,100.);
578// pythia->SetYRange(-8.,8.);
579// pythia->SetPhiRange(0.,360.);
580// pythia->SetPtHard(2.76,-1.0);
581// pythia->SwitchHFOff();
582// pythia->Init();
583// AddGenerator(pythia,"Pythia",1);
584
585}
586
587//-------------------------------------------------------------------
588void AliGenMUONCocktailpp::AddReso2Generator(Char_t* nameReso,
589 AliGenParam* const genReso,
590 Double_t sigmaReso,
591 Double_t polReso)
592{
593// add resonances to the cocktail
594 Double_t phiMin = fPhiMin*180./TMath::Pi();
595 Double_t phiMax = fPhiMax*180./TMath::Pi();
596
597 // first step: generation in 4pi
598 genReso->SetPtRange(0.,100.);
599 genReso->SetYRange(-8.,8.);
600 genReso->SetPhiRange(0.,360.);
601 genReso->SetForceDecay(fDecayModeResonance);
602 if (!gMC) genReso->SetDecayer(fDecayer);
603 genReso->Init(); // generation in 4pi
604// Ratios with respect to the reaction cross-section in the
605// kinematics limit of the MUONCocktail
606 Double_t ratioReso = sigmaReso / fSigmaReaction * genReso->GetRelativeArea(fPtMin,fPtMax,fYMin,fYMax,phiMin,phiMax);
607 if (!fSigmaSilent) {
608 AliInfo(Form("%s prod. cross-section in pp %5.3g b",nameReso,sigmaReso));
609 AliInfo(Form("%s prod. probability per collision in acceptance %5.3g",nameReso,ratioReso));
610 }
611// second step: generation in selected kinematical range
612 genReso->SetPtRange(fPtMin, fPtMax);
613 genReso->SetYRange(fYMin, fYMax);
614 genReso->SetPhiRange(phiMin, phiMax);
615 genReso->Init(); // generation in selected kinematical range
616
617 TVirtualMCDecayer *decReso = 0;
618 if(TMath::Abs(polReso) > 1.e-30){
619 AliInfo(Form("******Setting polarized decayer for %s''",nameReso));
620 if(fPolFrame==0) {
621 decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
622 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",polReso));
623 }
624 if(fPolFrame==1) {
625 decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
626 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",polReso));
627 }
628 if (decReso) {
629 decReso->SetForceDecay(kAll);
630 decReso->Init();
631 genReso->SetDecayer(decReso);
632 }
633 }
634
635 AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
636}
637
638//-------------------------------------------------------------------
639void AliGenMUONCocktailpp::Init()
640{
641 // Initialisation
642 TIter next(fEntries);
643 AliGenCocktailEntry *entry;
644 if (fStack) {
645 while((entry = (AliGenCocktailEntry*)next())) {
646 entry->Generator()->SetStack(fStack);
647 }
648 }
649}
650
651//_________________________________________________________________________
652void AliGenMUONCocktailpp::Generate()
653{
654// Generate event
655 TIter next(fEntries);
656 AliGenCocktailEntry *entry = 0;
657 AliGenCocktailEntry *preventry = 0;
658 AliGenerator* gen = 0;
659
660 if (fHeader) delete fHeader;
661 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
662
663 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
664
665// Generate the vertex position used by all generators
666 if(fVertexSmear == kPerEvent) Vertex();
667
668// Loop on primordialTrigger:
669// minimum muon multiplicity above a pt cut in a theta acceptance region
670
671 Bool_t primordialTrigger = kFALSE;
672 while(!primordialTrigger) {
673 //Reseting stack
674 AliRunLoader * runloader = AliRunLoader::Instance();
675 if (runloader)
676 if (runloader->Stack())
677 runloader->Stack()->Clean();
678 // Loop over generators and generate events
679 Int_t igen = 0;
680 Int_t npart = 0;
681 const char* genName = 0;
682 while((entry = (AliGenCocktailEntry*)next())) {
683 gen = entry->Generator();
684 genName = entry->GetName();
685 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
686 gen->SetTime(fTime);
687
688 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
689 gRandom->Poisson(entry->Rate());
690
691 if (npart > 0) {
692 igen++;
693 if (igen == 1) entry->SetFirst(0);
694 else entry->SetFirst((partArray->GetEntriesFast())+1);
695
696 gen->SetNumberParticles(npart);
697 gen->Generate();
698 entry->SetLast(partArray->GetEntriesFast());
699 preventry = entry;
700 }
701 }
702 next.Reset();
703
704
705// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
706// in the muon spectrometer acceptance
707 Int_t iPart;
708 fNGenerated++;
709 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
710 for(iPart=0; iPart<maxPart; iPart++){
711
712 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
713 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
714 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
715 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
716 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
717 (part->Pt()>fMuonPtCut) &&
718 (part->P()>fMuonPCut)) {
719 numberOfMuons++;
720 }
721 }
722 }
723 if (numberOfMuons >= fMuonMultiplicity) {
724 primordialTrigger = kTRUE;
725 fHeader->SetNProduced(maxPart);
726 }
727
728 }
729 fNSucceded++;
730
731 TArrayF eventVertex;
732 eventVertex.Set(3);
733 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
734
735 fHeader->SetPrimaryVertex(eventVertex);
736 fHeader->SetInteractionTime(fTime);
737
738 gAlice->SetGenEventHeader(fHeader);
739
740// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
741 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
742}
743
744