Added AliGenHepMCEventHeader class and put it in STEERBase library. Updated AliGenRea...
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
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   }
79 AliGenMUONCocktailpp* 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
106 ClassImp(AliGenMUONCocktailpp)  
107   
108 //________________________________________________________________________
109 AliGenMUONCocktailpp::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 //_________________________________________________________________________
339 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
340 {
341 // Destructor
342
343 }
344
345 //_________________________________________________________________________
346 void 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 //_________________________________________________________________________
363 void 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 //_________________________________________________________________________
393 void 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 (!TVirtualMC::GetMC()) 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 (!TVirtualMC::GetMC()) 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 //-------------------------------------------------------------------
609 void 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 (!TVirtualMC::GetMC()) 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 //-------------------------------------------------------------------
660 void 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 //_________________________________________________________________________
673 void 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