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