reduce histogram bin size for delta eta close to 0 from 0.05 to 0.01
[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//
ba762ab4 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
103ac317 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
00e6c5ee 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)
ba8bf3f5 31//-----------------
00e6c5ee 32// 2009: added polarization (L. Bianchi)
ba8bf3f5 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
ba762ab4 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
7455632e 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
4bac9bd3 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
ba762ab4 66
103ac317 67#include <TObjArray.h>
68#include <TParticle.h>
69#include <TF1.h>
aaa95f22 70#include <TVirtualMC.h>
103ac317 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"
aaa95f22 80#include "AliDecayer.h"
103ac317 81#include "AliLog.h"
b44c3901 82#include "AliGenCorrHF.h"
00e6c5ee 83#include "AliDecayerPolarized.h"
103ac317 84
85ClassImp(AliGenMUONCocktailpp)
86
87//________________________________________________________________________
88AliGenMUONCocktailpp::AliGenMUONCocktailpp()
1c56e311 89 :AliGenCocktail(),
aaa95f22 90 fDecayer(0),
91 fDecayModeResonance(kAll),
92 fDecayModePythia(kAll),
1c56e311 93 fMuonMultiplicity(0),
94 fMuonPtCut(0.),
f81a4aca 95 fMuonPCut(0.),
18bc0c8e 96 fMuonThetaMinCut(0.),
f81a4aca 97 fMuonThetaMaxCut(180.),
aaa95f22 98 fMuonOriginCut(-999.),
18bc0c8e 99 fNSucceded(0),
9ff768ee 100 fNGenerated(0),
7455632e 101 fCentralityBin(0),
4bac9bd3 102 fScaleJPsi(1),
103 fScaleCharmonia(1),
104 fScaleBottomonia(1),
105 fScaleCCbar(1),
106 fScaleBBbar(1),
00e6c5ee 107
108 fJpsiPol(0),
ba8bf3f5 109 fChic1Pol(0),
110 fChic2Pol(0),
00e6c5ee 111 fPsiPPol(0),
112 fUpsPol(0),
113 fUpsPPol(0),
114 fUpsPPPol(0),
115 fPolFrame(0),
116
ba762ab4 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
7455632e 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
4bac9bd3 149// gives 0.56 mub, sigmas for 2S & 3S obtained using LHCb data for ratios 2S/1S & 3S/1S
7455632e 150// -ccbar & bbbar: NLO pQCD computations - http://www-alice.gsi.de/ana/MNR/results.html
ba762ab4 151 fCMSEnergyTeVArray[0] = 7.00;
7455632e 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;
4bac9bd3 158 fSigmaUpsilonPArray[0] = 0.18e-6;
159 fSigmaUpsilonPPArray[0] = 0.08e-6;
ba762ab4 160 fSigmaCCbarArray[0] = 6.91e-3;
161 fSigmaBBbarArray[0] = 0.232e-3;
162
ba8bf3f5 163//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
7455632e 164// scaled down according to ccbar and bbbar cross-sections
ba762ab4 165 fCMSEnergyTeVArray[1] = 10.00;
7455632e 166 fSigmaReactionArray[1] = 0.070;
ba762ab4 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;
7455632e 176
ba8bf3f5 177//x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
7455632e 178// bottomonium from hep-ph/0311048 Tab.9, page 19 taken into account that
ba8bf3f5 179// feed-down from chib is included
ba762ab4 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;
7455632e 190 fSigmaBBbarArray[2] = 0.445e-3;
191
4bac9bd3 192// x-sections for Min. Bias p-Pb & Pb-p @ 2.76 TeV: charmonia and bottomonia
7455632e 193// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
194// and with Glauber scaling
4bac9bd3 195 fCMSEnergyTeVArray[3] = 2.00; // for 2.76 TeV
7455632e 196 fSigmaReactionArray[3] = 2.10;
4bac9bd3 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
7455632e 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
4bac9bd3 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
7455632e 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
4bac9bd3 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
ba762ab4 314
103ac317 315}
93a2041b 316
103ac317 317//_________________________________________________________________________
318AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
ba762ab4 319{
103ac317 320// Destructor
ba762ab4 321
322}
323
324//_________________________________________________________________________
325void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy)
103ac317 326{
ba762ab4 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];
103ac317 339}
340
341//_________________________________________________________________________
00e6c5ee 342void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
343 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
ba762ab4 344// setter for resonances polarization
00e6c5ee 345 if (strcmp(PolFrame,"kColSop")==0){
ba762ab4 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;
00e6c5ee 352 } else if (strcmp(PolFrame,"kHelicity")==0){
ba762ab4 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
00e6c5ee 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//_________________________________________________________________________
8058ead1 372void AliGenMUONCocktailpp::CreateCocktail()
103ac317 373{
ba762ab4 374// create and add resonances and open HF to the coctail
375 Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
376
4bac9bd3 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
103ac317 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
4bac9bd3 391// Cross sections in barns
103ac317 392
9ff768ee 393 Double_t sigmajpsi = fSigmaJPsi;
ba8bf3f5 394 Double_t sigmachic1 = fSigmaChic1;
395 Double_t sigmachic2 = fSigmaChic2;
9ff768ee 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;
3285b419 402
00e6c5ee 403// Cross sections corrected with the BR in mu+mu-
404// (only in case of use of AliDecayerPolarized)
405
ba762ab4 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;}
00e6c5ee 413
4bac9bd3 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
aaa95f22 426 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
427
ba762ab4 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
c478a271 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");
ba762ab4 452
7455632e 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");
4bac9bd3 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);
7455632e 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);
418e091f 468 } else {
7455632e 469 AliError("Initialisation failed, wrong cmsEnergy");
418e091f 470 return;
00e6c5ee 471 }
7455632e 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
4bac9bd3 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 %
7455632e 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
4bac9bd3 514// 11 centrality bins for Pb-Pb at 2.76 TeV: 0-5-10-20-30-40-50-60-70-80-90-100 %
7455632e 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);
00e6c5ee 533
7455632e 534 AddReso2Generator(nameUps,genupsilon,bhard[fCentralityBin]*sigmaupsilon,fUpsPol);
535 AddReso2Generator(nameUpsP,genupsilonP,bhard[fCentralityBin]*sigmaupsilonP,fUpsPPol);
536 AddReso2Generator(nameUpsPP,genupsilonPP,bhard[fCentralityBin]*sigmaupsilonPP,fUpsPPPol);
b44c3901 537
538//------------------------------------------------------------------
539// Generator of charm
ba762ab4 540 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
8058ead1 541 gencharm->SetMomentumRange(0,9999);
542 gencharm->SetForceDecay(kAll);
7455632e 543 Double_t ratioccbar = chard[fCentralityBin]*sigmaccbar/fSigmaReaction;
ba762ab4 544 if (!gMC) gencharm->SetDecayer(fDecayer);
8058ead1 545 gencharm->Init();
9ff768ee 546 if (!fSigmaSilent) {
00e6c5ee 547 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
3285b419 548 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
549 }
8058ead1 550 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
103ac317 551//------------------------------------------------------------------
b44c3901 552// Generator of beauty
ba762ab4 553 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
8058ead1 554 genbeauty->SetMomentumRange(0,9999);
555 genbeauty->SetForceDecay(kAll);
7455632e 556 Double_t ratiobbbar = bhard[fCentralityBin]*sigmabbbar/fSigmaReaction;
ba762ab4 557 if (!gMC) genbeauty->SetDecayer(fDecayer);
8058ead1 558 genbeauty->Init();
9ff768ee 559 if (!fSigmaSilent) {
00e6c5ee 560 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
3285b419 561 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
562 }
ba762ab4 563 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
b44c3901 564
8058ead1 565//-------------------------------------------------------------------
103ac317 566// Pythia generator
8058ead1 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);
ba8bf3f5 584
8058ead1 585}
103ac317 586
ba762ab4 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 }
418e091f 628 if (decReso) {
629 decReso->SetForceDecay(kAll);
630 decReso->Init();
631 genReso->SetDecayer(decReso);
632 }
ba762ab4 633 }
634
635 AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
636}
637
638//-------------------------------------------------------------------
8058ead1 639void AliGenMUONCocktailpp::Init()
640{
8058ead1 641 // Initialisation
aaa95f22 642 TIter next(fEntries);
643 AliGenCocktailEntry *entry;
644 if (fStack) {
645 while((entry = (AliGenCocktailEntry*)next())) {
646 entry->Generator()->SetStack(fStack);
647 }
648 }
103ac317 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
01c89ce1 660 if (fHeader) delete fHeader;
661 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
662
d94af0c1 663 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
103ac317 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
33c3c91a 674 AliRunLoader * runloader = AliRunLoader::Instance();
103ac317 675 if (runloader)
676 if (runloader->Stack())
34da8494 677 runloader->Stack()->Clean();
103ac317 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));
21391258 686 gen->SetTime(fTime);
103ac317 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
ba762ab4 704
103ac317 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++;
b44c3901 709 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22 710 for(iPart=0; iPart<maxPart; iPart++){
103ac317 711
aaa95f22 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) &&
f81a4aca 717 (part->Pt()>fMuonPtCut) &&
718 (part->P()>fMuonPCut)) {
aaa95f22 719 numberOfMuons++;
103ac317 720 }
aaa95f22 721 }
722 }
01c89ce1 723 if (numberOfMuons >= fMuonMultiplicity) {
724 primordialTrigger = kTRUE;
725 fHeader->SetNProduced(maxPart);
726 }
727
103ac317 728 }
729 fNSucceded++;
730
01c89ce1 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);
21391258 736 fHeader->SetInteractionTime(fTime);
01c89ce1 737
738 gAlice->SetGenEventHeader(fHeader);
739
aaa95f22 740// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 741 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
742}
743
744