]> 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
- Added threshold for spline systematics and different sys. scale factors below and...
[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 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()
86// S. Grigoryan
87
88#include <TObjArray.h>
89#include <TParticle.h>
90#include <TF1.h>
91#include <TVirtualMC.h>
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"
101#include "AliDecayer.h"
102#include "AliLog.h"
103#include "AliGenCorrHF.h"
104#include "AliDecayerPolarized.h"
105
106ClassImp(AliGenMUONCocktailpp)
107
108//________________________________________________________________________
109AliGenMUONCocktailpp::AliGenMUONCocktailpp()
110 :AliGenCocktail(),
111 fDecayer(0),
112 fDecayModeResonance(kAll),
113 fDecayModePythia(kAll),
114 fMuonMultiplicity(0),
115 fMuonPtCut(0.),
116 fMuonPCut(0.),
117 fMuonThetaMinCut(0.),
118 fMuonThetaMaxCut(180.),
119 fMuonOriginCut(-999.),
120 fNSucceded(0),
121 fNGenerated(0),
122 fCentralityBin(0),
123 fScaleJPsi(1),
124 fScaleCharmonia(1),
125 fScaleBottomonia(1),
126 fScaleCCbar(1),
127 fScaleBBbar(1),
128
129 fJpsiPol(0),
130 fChic1Pol(0),
131 fChic2Pol(0),
132 fPsiPPol(0),
133 fUpsPol(0),
134 fUpsPPol(0),
135 fUpsPPPol(0),
136 fPolFrame(0),
137
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
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
170// gives 0.56 mub, sigmas for 2S & 3S obtained using LHCb data for ratios 2S/1S & 3S/1S
171// -ccbar & bbbar: NLO pQCD computations - http://www-alice.gsi.de/ana/MNR/results.html
172 fCMSEnergyTeVArray[0] = 7.00;
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;
179 fSigmaUpsilonPArray[0] = 0.18e-6;
180 fSigmaUpsilonPPArray[0] = 0.08e-6;
181 fSigmaCCbarArray[0] = 6.91e-3;
182 fSigmaBBbarArray[0] = 0.232e-3;
183
184//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
185// scaled down according to ccbar and bbbar cross-sections
186 fCMSEnergyTeVArray[1] = 10.00;
187 fSigmaReactionArray[1] = 0.070;
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;
197
198//x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
199// bottomonium from hep-ph/0311048 Tab.9, page 19 taken into account that
200// feed-down from chib is included
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;
211 fSigmaBBbarArray[2] = 0.445e-3;
212
213// x-sections for Min. Bias p-Pb & Pb-p @ 2.76 TeV: charmonia and bottomonia
214// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
215// and with Glauber scaling
216 fCMSEnergyTeVArray[3] = 2.00; // for 2.76 TeV
217 fSigmaReactionArray[3] = 2.10;
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
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
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
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
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
335
336}
337
338//_________________________________________________________________________
339AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
340{
341// Destructor
342
343}
344
345//_________________________________________________________________________
346void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy)
347{
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];
360}
361
362//_________________________________________________________________________
363void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
364 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
365// setter for resonances polarization
366 if (strcmp(PolFrame,"kColSop")==0){
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;
373 } else if (strcmp(PolFrame,"kHelicity")==0){
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
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//_________________________________________________________________________
393void AliGenMUONCocktailpp::CreateCocktail()
394{
395// create and add resonances and open HF to the coctail
396 Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
397
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
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
412// Cross sections in barns
413
414 Double_t sigmajpsi = fSigmaJPsi;
415 Double_t sigmachic1 = fSigmaChic1;
416 Double_t sigmachic2 = fSigmaChic2;
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;
423
424// Cross sections corrected with the BR in mu+mu-
425// (only in case of use of AliDecayerPolarized)
426
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;}
434
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
447 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
448
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
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");
473
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");
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);
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);
489 } else {
490 AliError("Initialisation failed, wrong cmsEnergy");
491 return;
492 }
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
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 %
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
535// 11 centrality bins for Pb-Pb at 2.76 TeV: 0-5-10-20-30-40-50-60-70-80-90-100 %
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);
554
555 AddReso2Generator(nameUps,genupsilon,bhard[fCentralityBin]*sigmaupsilon,fUpsPol);
556 AddReso2Generator(nameUpsP,genupsilonP,bhard[fCentralityBin]*sigmaupsilonP,fUpsPPol);
557 AddReso2Generator(nameUpsPP,genupsilonPP,bhard[fCentralityBin]*sigmaupsilonPP,fUpsPPPol);
558
559//------------------------------------------------------------------
560// Generator of charm
561 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
562 gencharm->SetMomentumRange(0,9999);
563 gencharm->SetForceDecay(kAll);
564 Double_t ratioccbar = chard[fCentralityBin]*sigmaccbar/fSigmaReaction;
565 if (!gMC) gencharm->SetDecayer(fDecayer);
566 gencharm->Init();
567 if (!fSigmaSilent) {
568 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
569 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
570 }
571 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
572//------------------------------------------------------------------
573// Generator of beauty
574 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
575 genbeauty->SetMomentumRange(0,9999);
576 genbeauty->SetForceDecay(kAll);
577 Double_t ratiobbbar = bhard[fCentralityBin]*sigmabbbar/fSigmaReaction;
578 if (!gMC) genbeauty->SetDecayer(fDecayer);
579 genbeauty->Init();
580 if (!fSigmaSilent) {
581 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
582 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
583 }
584 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
585
586//-------------------------------------------------------------------
587// Pythia generator
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);
605
606}
607
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 }
649 if (decReso) {
650 decReso->SetForceDecay(kAll);
651 decReso->Init();
652 genReso->SetDecayer(decReso);
653 }
654 }
655
656 AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
657}
658
659//-------------------------------------------------------------------
660void AliGenMUONCocktailpp::Init()
661{
662 // Initialisation
663 TIter next(fEntries);
664 AliGenCocktailEntry *entry;
665 if (fStack) {
666 while((entry = (AliGenCocktailEntry*)next())) {
667 entry->Generator()->SetStack(fStack);
668 }
669 }
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
681 if (fHeader) delete fHeader;
682 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
683
684 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
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) {
694 //Reseting stack (if fMuonMultiplicity > 0 : S. Grigoryan, June 2012)
695 AliRunLoader * runloader = AliRunLoader::Instance();
696 // if (runloader)
697 if (runloader && fMuonMultiplicity > 0)
698 if (runloader->Stack())
699 runloader->Stack()->Clean();
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));
708 gen->SetTime(fTime);
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
726
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++;
731 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
732 for(iPart=0; iPart<maxPart; iPart++){
733
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) &&
739 (part->Pt()>fMuonPtCut) &&
740 (part->P()>fMuonPCut)) {
741 numberOfMuons++;
742 }
743 }
744 }
745 if (numberOfMuons >= fMuonMultiplicity) {
746 primordialTrigger = kTRUE;
747 fHeader->SetNProduced(maxPart);
748 }
749
750 }
751 fNSucceded++;
752
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);
758 fHeader->SetInteractionTime(fTime);
759
760 gAlice->SetGenEventHeader(fHeader);
761
762// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
763 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
764}
765
766