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