I have added Chic_1 and Chic_2 production in the muon cocktail.
[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 // Classe to create the coktail for physics with muons for pp at 14 TeV
20 // The followoing sources: 
21 // jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia
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 #include <TObjArray.h>
40 #include <TParticle.h>
41 #include <TF1.h>
42 #include <TVirtualMC.h>
43
44 #include "AliGenCocktailEventHeader.h"
45
46 #include "AliGenCocktailEntry.h"
47 #include "AliGenMUONCocktailpp.h"
48 #include "AliGenMUONlib.h"
49 #include "AliGenParam.h"
50 #include "AliMC.h"
51 #include "AliRun.h"
52 #include "AliStack.h"
53 #include "AliDecayer.h"
54 #include "AliLog.h"
55 #include "AliGenCorrHF.h"
56 #include "AliDecayerPolarized.h"
57
58 ClassImp(AliGenMUONCocktailpp)  
59   
60 //________________________________________________________________________
61 AliGenMUONCocktailpp::AliGenMUONCocktailpp()
62     :AliGenCocktail(),
63      fDecayer(0),
64      fDecayModeResonance(kAll),
65      fDecayModePythia(kAll),
66      fMuonMultiplicity(0),
67      fMuonPtCut(0.),
68      fMuonPCut(0.),
69      fMuonThetaMinCut(0.),
70      fMuonThetaMaxCut(180.),
71      fMuonOriginCut(-999.),
72      fNSucceded(0),
73      fNGenerated(0),
74
75      fJpsiPol(0), 
76      fChic1Pol(0), 
77      fChic2Pol(0), 
78      fPsiPPol(0), 
79      fUpsPol(0), 
80      fUpsPPol(0), 
81      fUpsPPPol(0),
82      fPolFrame(0),
83
84 // x-sections for pp @ 7 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
85 // bottomnium as for 10 TeV
86 /*     fCMSEnergy(7),
87      fSigmaReaction(0.0695),
88      fSigmaJPsi(21.8e-6),
89      fSigmaChic1(21.1e-6),
90      fSigmaChic2(34.9e-6),
91      fSigmaPsiP(4.93e-6),
92      fSigmaUpsilon(0.463e-6),
93      fSigmaUpsilonP(0.154e-6),
94      fSigmaUpsilonPP(0.0886e-6),
95      fSigmaCCbar(6.91e-3),
96      fSigmaBBbar(0.232e-3),
97 */
98
99 //x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
100 // scaled down according to ccbar and bbar cross-sections
101      fCMSEnergy(10),
102      fSigmaReaction(0.0695),
103      fSigmaJPsi(26.06e-6),
104      fSigmaChic1(25.18e-6),
105      fSigmaChic2(41.58e-6),
106      fSigmaPsiP(5.88e-6),
107      fSigmaUpsilon(0.658e-6),
108      fSigmaUpsilonP(0.218e-6),
109      fSigmaUpsilonPP(0.122e-6),
110      fSigmaCCbar(8.9e-3),
111      fSigmaBBbar(0.33e-3),
112
113
114 //x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
115 // bottomonium from hep-ph/0311048 Tab.9, page 19 taken inton account that 
116 // feed-down from chib is included
117 /*     fCMSEnergy(14),
118      fSigmaReaction(0.070),
119      fSigmaJPsi(32.9e-6),
120      fSigmaChic1(31.8e-6),
121      fSigmaChic2(52.5e-6),
122      fSigmaPsiP(7.43e-6),
123      fSigmaUpsilon(0.989e-6),
124      fSigmaUpsilonP(0.502e-6),
125      fSigmaUpsilonPP(0.228e-6),
126      fSigmaCCbar(11.2e-3),
127      fSigmaBBbar(0.51e-3),
128 */
129      fSigmaSilent(kFALSE)
130 {
131 // Constructor
132 }
133
134 //_________________________________________________________________________
135 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
136 // Destructor
137 {
138     
139 }
140
141 //_________________________________________________________________________
142 void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
143                                 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
144    
145    if (strcmp(PolFrame,"kColSop")==0){
146       fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
147       fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
148       fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
149       fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
150       fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
151       fPolFrame = 0;
152    } else if (strcmp(PolFrame,"kHelicity")==0){
153         fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
154         fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
155         fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
156         fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
157         fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
158         fPolFrame = 1;
159    
160    } else {
161        AliInfo(Form("The polarization frame is not valid"));
162        AliInfo(Form("No polarization will be set"));
163        fJpsiPol=0.;
164        fPsiPPol=0.;
165        fUpsPol=0.;
166        fUpsPPol=0.;
167        fUpsPPPol=0.;
168    }
169 }
170
171 //_________________________________________________________________________
172 void AliGenMUONCocktailpp::CreateCocktail()
173 {
174     
175 // These limits are only used for renormalization of quarkonia cross section
176 // Pythia events are generated in 4pi  
177     Double_t ptMin  = fPtMin;
178     Double_t ptMax  = fPtMax;
179     Double_t yMin   = fYMin;;
180     Double_t yMax   = fYMax;;
181     Double_t phiMin = fPhiMin*180./TMath::Pi();
182     Double_t phiMax = fPhiMax*180./TMath::Pi();
183     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));
184     
185 // Ratios with respect to the reaction cross-section in the 
186 // kinematics limit of the MUONCocktail
187     Double_t ratiojpsi;
188     Double_t ratiochic1;
189     Double_t ratiochic2;
190     Double_t ratiopsiP;
191     Double_t ratioupsilon;
192     Double_t ratioupsilonP;
193     Double_t ratioupsilonPP;
194     Double_t ratioccbar;
195     Double_t ratiobbbar;
196
197 // Beam energy
198     Int_t cmsEnergy = fCMSEnergy; 
199
200 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and 
201 // corrected from feed down of higher resonances 
202
203     Double_t sigmajpsi      = fSigmaJPsi;  
204     Double_t sigmachic1     = fSigmaChic1;  
205     Double_t sigmachic2     = fSigmaChic2;  
206     Double_t sigmapsiP      = fSigmaPsiP;  
207     Double_t sigmaupsilon   = fSigmaUpsilon;  
208     Double_t sigmaupsilonP  = fSigmaUpsilonP;  
209     Double_t sigmaupsilonPP = fSigmaUpsilonPP;
210     Double_t sigmaccbar     = fSigmaCCbar;
211     Double_t sigmabbbar     = fSigmaBBbar;
212
213 // Cross sections corrected with the BR in mu+mu-
214 // (only in case of use of AliDecayerPolarized)
215
216     if(fJpsiPol  != 0)  {sigmajpsi      = fSigmaJPsi*0.0593;}
217     if(fChic1Pol != 0)  {sigmachic1     = fSigmaChic1*0.;} // tb consistent
218     if(fChic2Pol != 0)  {sigmachic2     = fSigmaChic2*0.;} // tb consistent 
219     if(fPsiPPol  != 0)  {sigmapsiP      = fSigmaPsiP*0.0075;}
220     if(fUpsPol   != 0)  {sigmaupsilon   = fSigmaUpsilon*0.0248;}
221     if(fUpsPPol  != 0)  {sigmaupsilonP  = fSigmaUpsilonP*0.0193;}
222     if(fUpsPPPol != 0)  {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
223
224 // MinBias NN cross section @ pp 14 TeV -PR  Vol II p:473
225     Double_t sigmaReaction = fSigmaReaction;
226
227     Int_t eincStart = 10;
228
229     AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
230
231 // Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
232 //----------------------------------------------------------------------
233 // Jpsi generator
234     AliGenParam * genjpsi;
235     if(cmsEnergy == eincStart){
236         genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
237     } else {
238         genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi");
239     }
240 // first step: generation in 4pi
241     genjpsi->SetPtRange(0.,100.);
242     genjpsi->SetYRange(-8.,8.);
243     genjpsi->SetPhiRange(0.,360.);
244     genjpsi->SetForceDecay(fDecayModeResonance);
245     //if (!gMC) genjpsi->SetDecayer(fDecayer);
246
247     genjpsi->Init();  // generation in 4pi
248     ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
249     if (!fSigmaSilent) {
250       AliInfo(Form("jpsi prod. cross-section in pp %5.3g b",sigmajpsi));
251       AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
252     }
253 // second step: generation in selected kinematical range
254     genjpsi->SetPtRange(ptMin, ptMax);  
255     genjpsi->SetYRange(yMin, yMax);
256     genjpsi->SetPhiRange(phiMin, phiMax);
257     genjpsi->Init(); // generation in selected kinematic range
258
259     TVirtualMCDecayer *Jpsidec = 0;
260     if(fJpsiPol != 0){
261     AliInfo(Form("******Setting polarized decayer for J/psi"));
262     if(fPolFrame==0) {
263          Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
264          AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fJpsiPol));
265            }
266     if(fPolFrame==1) {
267          Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
268          AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fJpsiPol));
269            }
270     Jpsidec->SetForceDecay(kAll);
271     Jpsidec->Init();
272     genjpsi->SetDecayer(Jpsidec);
273     }
274
275     AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
276
277
278 //----------------------------------------------------------------------
279 // Chic1 generator
280     AliGenParam * genchic1;
281     if(cmsEnergy == eincStart){
282         genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 10", "Chic1");
283     } else {
284         genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp ", "Chic1");
285     }
286 // first step: generation in 4pi
287     genchic1->SetPtRange(0.,100.);
288     genchic1->SetYRange(-8.,8.);
289     genchic1->SetPhiRange(0.,360.);
290     genchic1->SetForceDecay(fDecayModeResonance);
291     //if (!gMC) genchic1->SetDecayer(fDecayer);
292
293     genchic1->Init();  // generation in 4pi
294     ratiochic1 = sigmachic1 / sigmaReaction * genchic1->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
295     if (!fSigmaSilent) {
296       AliInfo(Form("chic1 prod. cross-section in pp %5.3g b",sigmachic1));
297       AliInfo(Form("chic1 prod. probability per collision in acceptance %5.3g",ratiochic1));
298     }
299 // second step: generation in selected kinematical range
300     genchic1->SetPtRange(ptMin, ptMax);  
301     genchic1->SetYRange(yMin, yMax);
302     genchic1->SetPhiRange(phiMin, phiMax);
303     genchic1->Init(); // generation in selected kinematic range
304
305     TVirtualMCDecayer *Chic1dec = 0;
306     if(fChic1Pol != 0){
307     AliInfo(Form("******Setting polarized decayer for Chic1"));
308     if(fPolFrame==0) {
309          Chic1dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
310          AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic1Pol));
311            }
312     if(fPolFrame==1) {
313          Chic1dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
314          AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic1Pol));
315            }
316     Chic1dec->SetForceDecay(kAll);
317     Chic1dec->Init();
318     genchic1->SetDecayer(Chic1dec);
319     }
320
321     AddGenerator(genchic1, "Chic1", ratiochic1); // Adding Generator
322
323 //----------------------------------------------------------------------
324 // Chic2 generator
325     AliGenParam * genchic2;
326     if(cmsEnergy == eincStart){
327         genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 10", "Chic2");
328     } else {
329         genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp ", "Chic2");
330     }
331 // first step: generation in 4pi
332     genchic2->SetPtRange(0.,100.);
333     genchic2->SetYRange(-8.,8.);
334     genchic2->SetPhiRange(0.,360.);
335     genchic2->SetForceDecay(fDecayModeResonance);
336     //if (!gMC) genchic1->SetDecayer(fDecayer);
337
338     genchic2->Init();  // generation in 4pi
339     ratiochic2 = sigmachic2 / sigmaReaction * genchic2->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
340     if (!fSigmaSilent) {
341       AliInfo(Form("chic2 prod. cross-section in pp %5.3g b",sigmachic2));
342       AliInfo(Form("chic2 prod. probability per collision in acceptance %5.3g",ratiochic2));
343     }
344 // second step: generation in selected kinematical range
345     genchic2->SetPtRange(ptMin, ptMax);  
346     genchic2->SetYRange(yMin, yMax);
347     genchic2->SetPhiRange(phiMin, phiMax);
348     genchic2->Init(); // generation in selected kinematic range
349
350     TVirtualMCDecayer *Chic2dec = 0;
351     if(fChic2Pol != 0){
352     AliInfo(Form("******Setting polarized decayer for Chic2"));
353     if(fPolFrame==0) {
354          Chic2dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
355          AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic2Pol));
356            }
357     if(fPolFrame==1) {
358          Chic2dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
359          AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic2Pol));
360            }
361     Chic2dec->SetForceDecay(kAll);
362     Chic2dec->Init();
363     genchic2->SetDecayer(Chic2dec);
364     }
365
366     AddGenerator(genchic2, "Chic2", ratiochic2); // Adding Generator
367
368 //------------------------------------------------------------------
369 // Psi prime generator
370     AliGenParam * genpsiP;
371     if(cmsEnergy == eincStart){    
372         genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
373     } else {
374         genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
375     }
376 // first step: generation in 4pi
377     genpsiP->SetPtRange(0.,100.);
378     genpsiP->SetYRange(-8.,8.);
379     genpsiP->SetPhiRange(0.,360.);
380     genpsiP->SetForceDecay(fDecayModeResonance);
381     //if (!gMC) genpsiP->SetDecayer(fDecayer);
382
383     genpsiP->Init();  // generation in 4pi
384     ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
385     if (!fSigmaSilent) {
386       AliInfo(Form("psiP prod. cross-section in pp %5.3g b",sigmapsiP));
387       AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
388     }
389 // second step: generation in selected kinematical range
390     genpsiP->SetPtRange(ptMin, ptMax);  
391     genpsiP->SetYRange(yMin, yMax);
392     genpsiP->SetPhiRange(phiMin, phiMax);
393     genpsiP->Init(); // generation in selected kinematic range
394
395      TVirtualMCDecayer *PsiPdec = 0;
396      if(fPsiPPol != 0){
397       AliInfo(Form("******Setting polarized decayer for psi'"));
398       if(fPolFrame==0) {
399         PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
400         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fPsiPPol));
401           }
402       if(fPolFrame==1) {
403         PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
404         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fPsiPPol));
405           }
406       PsiPdec->SetForceDecay(kAll);
407       PsiPdec->Init();
408       genpsiP->SetDecayer(PsiPdec);
409      }
410
411     AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
412
413 //------------------------------------------------------------------
414 // Upsilon 1S generator
415     AliGenParam * genupsilon;
416     if(cmsEnergy == eincStart) {
417         genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
418     } else {
419         genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
420     }
421 // first step: generation in 4pi
422     genupsilon->SetPtRange(0.,100.); 
423     genupsilon->SetYRange(-8.,8.);
424     genupsilon->SetPhiRange(0.,360.);
425     genupsilon->SetForceDecay(fDecayModeResonance);
426     //if (!gMC) genupsilon->SetDecayer(fDecayer);
427     genupsilon->Init();  // generation in 4pi
428     ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
429     if (!fSigmaSilent) {
430       AliInfo(Form("Upsilon 1S prod. cross-section in pp %5.3g b",sigmaupsilon));
431       AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
432     }
433 // second step: generation in selected kinematical range
434     genupsilon->SetPtRange(ptMin, ptMax);  
435     genupsilon->SetYRange(yMin, yMax);
436     genupsilon->SetPhiRange(phiMin, phiMax);
437     genupsilon->Init(); // generation in selected kinematical range
438      
439      TVirtualMCDecayer *Upsdec = 0;
440      if(fUpsPol != 0){
441       AliInfo(Form("******Setting polarized decayer for Upsilon"));
442       if(fPolFrame==0) {
443         Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
444         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPol));
445           }
446       if(fPolFrame==1) {
447         Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
448         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPol));
449           }
450       Upsdec->SetForceDecay(kAll);
451       Upsdec->Init();
452       genupsilon->SetDecayer(Upsdec);
453      }
454
455     AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
456 //------------------------------------------------------------------
457 // Upsilon 2S generator
458     AliGenParam * genupsilonP;
459     if(cmsEnergy == eincStart){
460         genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
461     } else {
462         genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
463     }
464         
465 // first step: generation in 4pi
466     genupsilonP->SetPtRange(0.,100.);
467     genupsilonP->SetYRange(-8.,8.);
468     genupsilonP->SetPhiRange(0.,360.);
469     genupsilonP->SetForceDecay(fDecayModeResonance);
470     //if (!gMC) genupsilonP->SetDecayer(fDecayer);  
471     genupsilonP->Init();  // generation in 4pi
472     ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
473     if (!fSigmaSilent) {
474       AliInfo(Form("Upsilon 2S prod. cross-section in pp %5.3g b",sigmaupsilonP));
475       AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
476     }
477 // second step: generation in the kinematical range
478     genupsilonP->SetPtRange(ptMin, ptMax);  
479     genupsilonP->SetYRange(yMin, yMax);
480     genupsilonP->SetPhiRange(phiMin, phiMax);
481     genupsilonP->Init(); // generation in selected kinematical range
482
483      TVirtualMCDecayer *UpsPdec = 0;
484      if(fUpsPPol != 0){
485       AliInfo(Form("******Setting polarized decayer for Upsilon'"));
486       if(fPolFrame==0) {
487         UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
488         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPol));
489           }
490       if(fPolFrame==1) {
491         UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
492         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPol));
493           }
494       UpsPdec->SetForceDecay(kAll);
495       UpsPdec->Init();
496       genupsilonP->SetDecayer(UpsPdec);
497      }
498     
499     AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
500     
501 //------------------------------------------------------------------
502 // Upsilon 3S generator
503     AliGenParam * genupsilonPP;
504     if(cmsEnergy == eincStart){
505         genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
506     }
507     else {
508         genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");    
509     }
510
511 // first step: generation in 4pi
512     genupsilonPP->SetPtRange(0.,100.); 
513     genupsilonPP->SetYRange(-8.,8.);
514     genupsilonPP->SetPhiRange(0.,360.);
515     genupsilonPP->SetForceDecay(fDecayModeResonance);
516     //if (!gMC) genupsilonPP->SetDecayer(fDecayer);  
517     genupsilonPP->Init();  // generation in 4pi
518     ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
519     if (!fSigmaSilent) {
520       AliInfo(Form("Upsilon 3S prod. cross-section in pp %5.3g b",sigmaupsilonPP));
521       AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
522     }
523 // second step: generation in selected kinematical range
524     genupsilonPP->SetPtRange(ptMin, ptMax);  
525     genupsilonPP->SetYRange(yMin, yMax);
526     genupsilonPP->SetPhiRange(phiMin, phiMax);
527     genupsilonPP->Init(); // generation in selected kinematical range
528
529      TVirtualMCDecayer *UpsPPdec = 0;
530      if(fUpsPPPol != 0){
531       AliInfo(Form("******Setting polarized decayer for Upsilon''"));
532       if(fPolFrame==0) {
533         UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
534         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPPol));
535           }
536       if(fPolFrame==1) {
537         UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
538         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPPol));
539           }
540       UpsPPdec->SetForceDecay(kAll);
541       UpsPPdec->Init();
542       genupsilonPP->SetDecayer(UpsPPdec);
543      }
544     
545     AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
546
547 //------------------------------------------------------------------
548 // Generator of charm
549     
550     AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);  
551     gencharm->SetMomentumRange(0,9999);
552     gencharm->SetForceDecay(kAll);
553     ratioccbar = sigmaccbar/sigmaReaction;
554     //if (!gMC) gencharm->SetDecayer(fDecayer);  
555     gencharm->Init();
556     if (!fSigmaSilent) {
557       AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
558       AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
559     }
560     AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
561
562 //------------------------------------------------------------------
563 // Generator of beauty
564
565     AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);  
566     genbeauty->SetMomentumRange(0,9999);
567     genbeauty->SetForceDecay(kAll);
568     ratiobbbar = sigmabbbar/sigmaReaction;
569     //if (!gMC) genbeauty->SetDecayer(fDecayer);  
570     genbeauty->Init();
571     if (!fSigmaSilent) {
572       AliInfo(Form("b-bbar prod. cross-section in pp  %5.3g b",sigmabbbar));
573       AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
574     }
575     AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar); 
576
577 //-------------------------------------------------------------------
578 // Pythia generator
579 //
580 // This has to go into the Config.C
581 //
582 //    AliGenPythia *pythia = new AliGenPythia(1);
583 //    pythia->SetProcess(kPyMbMSEL1);
584 //    pythia->SetStrucFunc(kCTEQ5L);
585 //    pythia->SetEnergyCMS(14000.);
586 //    AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
587 //    Decay_t dt = gener->GetDecayModePythia();
588 //    pythia->SetForceDecay(dt);
589 //    pythia->SetPtRange(0.,100.);
590 //    pythia->SetYRange(-8.,8.);
591 //    pythia->SetPhiRange(0.,360.);
592 //    pythia->SetPtHard(2.76,-1.0);
593 //    pythia->SwitchHFOff();
594 //    pythia->Init(); 
595 //    AddGenerator(pythia,"Pythia",1);
596
597 }
598
599 void AliGenMUONCocktailpp::Init()
600 {
601     //
602     // Initialisation
603     TIter next(fEntries);
604     AliGenCocktailEntry *entry;
605     if (fStack) {
606         while((entry = (AliGenCocktailEntry*)next())) {
607             entry->Generator()->SetStack(fStack);
608         }
609     }
610 }
611
612 //_________________________________________________________________________
613 void AliGenMUONCocktailpp::Generate()
614 {
615 // Generate event 
616     TIter next(fEntries);
617     AliGenCocktailEntry *entry = 0;
618     AliGenCocktailEntry *preventry = 0;
619     AliGenerator* gen = 0;
620
621     if (fHeader) delete fHeader;
622     fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
623
624     const TObjArray *partArray = gAlice->GetMCApp()->Particles();
625     
626 // Generate the vertex position used by all generators    
627     if(fVertexSmear == kPerEvent) Vertex();
628
629 // Loop on primordialTrigger: 
630 // minimum muon multiplicity above a pt cut in a theta acceptance region
631
632     Bool_t primordialTrigger = kFALSE;
633     while(!primordialTrigger) {
634         //Reseting stack
635         AliRunLoader * runloader = AliRunLoader::Instance();
636         if (runloader)
637             if (runloader->Stack())
638                 runloader->Stack()->Clean();
639         // Loop over generators and generate events
640         Int_t igen = 0;
641         Int_t npart = 0;        
642         const char* genName = 0;
643         while((entry = (AliGenCocktailEntry*)next())) {
644             gen = entry->Generator();
645             genName = entry->GetName();
646             gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
647
648             npart = (strcmp(genName,"Pythia") == 0) ? 1 :
649                 gRandom->Poisson(entry->Rate());
650
651             if (npart > 0) {
652                 igen++; 
653                 if (igen == 1) entry->SetFirst(0);              
654                 else  entry->SetFirst((partArray->GetEntriesFast())+1);
655                 
656                 gen->SetNumberParticles(npart);         
657                 gen->Generate();
658                 entry->SetLast(partArray->GetEntriesFast());
659                 preventry = entry;
660             }
661         }  
662         next.Reset();
663
664 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut 
665 // in the muon spectrometer acceptance
666         Int_t iPart;
667         fNGenerated++;
668         Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
669         for(iPart=0; iPart<maxPart; iPart++){      
670             
671           TParticle *part = gAlice->GetMCApp()->Particle(iPart);
672           if ( TMath::Abs(part->GetPdgCode()) == 13 ){  
673             if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
674                (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
675                (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
676                (part->Pt()>fMuonPtCut) &&
677                (part->P()>fMuonPCut)) {
678               numberOfMuons++;
679             }
680           }
681         }       
682         if (numberOfMuons >= fMuonMultiplicity) {
683           primordialTrigger = kTRUE;
684           fHeader->SetNProduced(maxPart);
685         }
686
687     }
688     fNSucceded++;
689
690     TArrayF eventVertex;
691     eventVertex.Set(3);
692     for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
693
694     fHeader->SetPrimaryVertex(eventVertex);
695
696     gAlice->SetGenEventHeader(fHeader);
697
698 //     AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
699     AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
700 }
701
702