Adding cut on invariant mass values in the pre-triiger decision of the AliMUONCocktai...
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktail.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 MUON coktail for physics in the Alice muon spectrometer
20 // The followoing muons sources are included in this cocktail: 
21 //     jpsi, upsilon, non-correlated open and beauty.
22 // The free parameeters are :
23 //      pp reaction cross-section
24 //      production cross-sections in pp collisions and 
25 //      branching ratios in the muon channel and shadowing
26 // Hard probes are supposed to scale with Ncoll and hadronic muon production with (0.8Ncoll+0.2*Npart)
27 // There is a primordial trigger which requires :
28 //      a minimum muon multiplicity above a pT cut in a theta acceptance cone
29 //
30
31 #include <TObjArray.h>
32 #include <TParticle.h>
33 #include <TF1.h>
34
35 #include "AliGenCocktailEntry.h"
36 #include "AliGenMUONCocktail.h"
37 #include "AliGenMUONlib.h"
38 #include "AliFastGlauber.h"
39 #include "AliGenParam.h"
40 #include "AliMC.h"
41 #include "AliRun.h"
42 #include "AliStack.h"
43 #include "AliLog.h"
44
45 ClassImp(AliGenMUONCocktail)
46   
47   
48 //________________________________________________________________________
49 AliGenMUONCocktail::AliGenMUONCocktail()
50   :AliGenCocktail()
51 {
52 // Constructor
53   fFastGlauber =0x0;
54   fTotalRate =0;  
55   fNSucceded=0; 
56   fNGenerated=0; 
57   fMuonMultiplicity=1;
58   fMuonPtCut= 1.;
59   fMuonThetaMinCut=171.; 
60   fMuonThetaMaxCut=178.;
61   fLowImpactParameter = 0.;
62   fHighImpactParameter = 5.;
63   fAverageImpactParameter =0.;
64   fNumberOfCollisions   = 0.; 
65   fNumberOfParticipants = 0.;
66   fHadronicMuons = kTRUE;
67   fInvMassCut = kFALSE;
68   fInvMassMinCut = 0.;
69   fInvMassMaxCut = 100.;
70 }
71 //_________________________________________________________________________
72 AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
73     AliGenCocktail(cocktail)
74 {
75 // Copy constructor
76   fFastGlauber =0x0;
77   fTotalRate =0; 
78   fNSucceded=0; 
79   fNGenerated=0; 
80   fMuonMultiplicity=1;
81   fMuonPtCut= 1.;
82   fMuonThetaMinCut=171.; 
83   fMuonThetaMaxCut=178.;  
84   fLowImpactParameter = 0.;
85   fHighImpactParameter = 5.;
86   fAverageImpactParameter =0.;
87   fNumberOfCollisions = 0.; 
88   fNumberOfParticipants = 0.;
89   fHadronicMuons = kTRUE;
90   fInvMassCut = kFALSE;
91   fInvMassMinCut = 0.;
92   fInvMassMaxCut = 100.;
93 }
94 //_________________________________________________________________________
95 AliGenMUONCocktail::~AliGenMUONCocktail()
96 {
97 // Destructor
98   if (fFastGlauber) delete fFastGlauber;
99 }
100
101 //_________________________________________________________________________
102 void AliGenMUONCocktail::Init()
103 {
104   // NN cross section
105   Double_t sigmaReaction =   0.072;   // MinBias NN cross section for PbPb LHC energies  http://arxiv.org/pdf/nucl-ex/0302016
106
107   // Initialising Fast Glauber object
108   fFastGlauber = new AliFastGlauber();
109   fFastGlauber->SetPbPbLHC();
110   fFastGlauber->SetNNCrossSection(sigmaReaction*1000.); //Expected NN cross-section in mb at LHC with diffractive part http://arxiv.org/pdf/nucl-ex/0302016 )
111   fFastGlauber->Init(1);
112  
113   // Calculating average number of collisions
114   Int_t ib=0;
115   Int_t ibmax=10000;
116   Double_t b       = 0.;
117   fAverageImpactParameter=0.;
118   fNumberOfCollisions   = 0.;
119   fNumberOfParticipants = 0.;
120   for(ib=0; ib<ibmax; ib++) {
121     b = fFastGlauber->GetRandomImpactParameter(fLowImpactParameter,fHighImpactParameter);
122     fAverageImpactParameter+=b;
123     fNumberOfCollisions    += fFastGlauber->GetNumberOfCollisions( b )/(1.-TMath::Exp(-fFastGlauber->GetNumberOfCollisions(b)));
124     fNumberOfParticipants  += fFastGlauber->GetNumberOfParticipants( b );
125   } 
126   fAverageImpactParameter/= ((Double_t) ibmax);
127   fNumberOfCollisions    /= ((Double_t) ibmax);
128   fNumberOfParticipants  /= ((Double_t) ibmax);;
129   AliInfo(Form("<b>=%4.2f, <Ncoll>=%5.1f and and <Npart>=%5.1f",b, fNumberOfCollisions, fNumberOfParticipants));
130
131   // Estimating shadowing on charm a beaty production
132   // -----------------------------------------------------
133   // Extrapolation of the cross sections from $p-p$ to \mbox{Pb--Pb}
134   // interactions
135   // is done by means of the Glauber model. For the impact parameter dependence
136   // of the shadowing factor we use a simple formula:
137   // $C_{sh}(b) = C_{sh}(0) + (1 - C_{sh}(0))(b/16~fm)4$,
138   // motivated by the theoretical predictions (see e.g.
139   // V. Emelyanov et al., Phys. Rev. C61, 044904 (2000)) and HIJING
140   // simulations showing an almost flat behaviour
141   // up to 10~$fm$ and a rapid increase to 1 for larger impact parameters.
142   // C_{sh}(0)  = 0.60 for Psi and 0.76 for Upsilon (Smba communication). 
143   // for open charm and beauty is 0.65 and 0.84
144   // ----------------------------------------------------- 
145   Double_t charmshadowing       = 0.65 + (1.0-0.65)*TMath::Power(fAverageImpactParameter/16.,4);
146   Double_t beautyshadowing      = 0.84 + (1.0-0.84)*TMath::Power(fAverageImpactParameter/16.,4);
147   Double_t charmoniumshadowing  = 0.60 + (1.0-0.60)*TMath::Power(fAverageImpactParameter/16.,4);
148   Double_t beautoniumshadowing  = 0.76 + (1.0-0.76)*TMath::Power(fAverageImpactParameter/16.,4);
149   if (fAverageImpactParameter>16.) {
150     charmoniumshadowing = 1.0;
151     beautoniumshadowing = 1.0;
152     charmshadowing = 1.0;
153     beautyshadowing= 1.0;
154   }
155   AliInfo(Form("Shadowing for charmonium and beautonium production are %4.2f and %4.2f respectively",charmoniumshadowing,beautoniumshadowing));
156   AliInfo(Form("Shadowing for charm and beauty production are %4.2f and %4.2f respectively",charmshadowing,beautyshadowing));
157
158   // Defining MUON physics cocktail
159   // Kinematical limits for particle generation
160   Double_t ptMin  = fPtMin;
161   Double_t ptMax  = fPtMax;
162   Double_t yMin   = fYMin;;
163   Double_t yMax   = fYMax;;
164   Double_t phiMin = fPhiMin*180./TMath::Pi();
165   Double_t phiMax = fPhiMax*180./TMath::Pi();
166   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));
167
168   // Generating J/Psi Physics 
169   // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
170   AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi");
171   genjpsi->SetPtRange(0,100.); // 4pi generation
172   genjpsi->SetYRange(-8.,8);
173   genjpsi->SetPhiRange(0.,360.);
174   genjpsi->SetForceDecay(kDiMuon);
175   genjpsi->SetTrackingFlag(1);
176   // Calculation of the particle multiplicity per event in the muonic channel
177   Double_t ratiojpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
178   Double_t sigmajpsi  = 31.0e-6 * charmoniumshadowing;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
179   Double_t brjpsi     = 0.0588;              // Branching Ratio for JPsi PDG PRC15 (200)
180   genjpsi->Init();  // Generating pT and Y parametrsation for the 4pi
181   ratiojpsi = sigmajpsi * brjpsi * fNumberOfCollisions / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
182   AliInfo(Form("Jpsi production cross-section in pp with shadowing %5.3g barns",sigmajpsi));
183   AliInfo(Form("Jpsi production probability per collisions in acceptance via the muonic channel %5.3g",ratiojpsi));
184   // Generation in the kinematical limits of AliGenMUONCocktail
185   genjpsi->SetPtRange(ptMin, ptMax);  
186   genjpsi->SetYRange(yMin, yMax);
187   genjpsi->SetPhiRange(phiMin, phiMax);
188   genjpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
189   // Adding Generator 
190   AddGenerator(genjpsi, "Jpsi", ratiojpsi);
191   fTotalRate+=ratiojpsi;
192
193  // Generating Psi prime Physics
194  // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
195   AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF scaled", "PsiP");
196   genpsiP->SetPtRange(0,100.);// 4pi generation
197   genpsiP->SetYRange(-8.,8);
198   genpsiP->SetPhiRange(0.,360.);
199   genpsiP->SetForceDecay(kDiMuon);
200   genpsiP->SetTrackingFlag(1);
201   // Calculation of the paritcle multiplicity per event in the muonic channel
202   Double_t ratiopsiP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
203   Double_t sigmapsiP     = 4.68e-6 * charmoniumshadowing;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
204   Double_t brpsiP        = 0.0103;           // Branching Ratio for PsiP
205   genpsiP->Init();  // Generating pT and Y parametrsation for the 4pi
206   ratiopsiP = sigmapsiP * brpsiP * fNumberOfCollisions / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
207   AliInfo(Form("Psi prime production cross-section in pp with shadowing %5.3g barns",sigmapsiP));
208   AliInfo(Form("Psi prime production probability per collisions in acceptance via the muonic channel %5.3g",ratiopsiP));
209   // Generation in the kinematical limits of AliGenMUONCocktail
210   genpsiP->SetPtRange(ptMin, ptMax);  
211   genpsiP->SetYRange(yMin, yMax);
212   genpsiP->SetPhiRange(phiMin, phiMax);
213   genpsiP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
214   // Adding Generator 
215   AddGenerator(genpsiP, "PsiP", ratiopsiP);
216   fTotalRate+=ratiopsiP;
217
218   // Generating Upsilon Physics
219   // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
220   AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF scaled", "Upsilon");  
221   genupsilon->SetPtRange(0,100.);  
222   genupsilon->SetYRange(-8.,8);
223   genupsilon->SetPhiRange(0.,360.);
224   genupsilon->SetForceDecay(kDiMuon);
225   genupsilon->SetTrackingFlag(1);
226   Double_t ratioupsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
227   Double_t sigmaupsilon     = 0.501e-6 * beautoniumshadowing;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
228   Double_t brupsilon        = 0.0248;  // Branching Ratio for Upsilon
229   genupsilon->Init();  // Generating pT and Y parametrsation for the 4pi
230   ratioupsilon = sigmaupsilon * brupsilon * fNumberOfCollisions / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
231   AliInfo(Form("Upsilon 1S production cross-section in pp with shadowing %5.3g barns",sigmaupsilon));
232   AliInfo(Form("Upsilon 1S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilon));
233   genupsilon->SetPtRange(ptMin, ptMax);  
234   genupsilon->SetYRange(yMin, yMax);
235   genupsilon->SetPhiRange(phiMin, phiMax);
236   genupsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
237   AddGenerator(genupsilon,"Upsilon", ratioupsilon);
238   fTotalRate+=ratioupsilon;
239
240   // Generating UpsilonP Physics
241   // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
242   AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF Scaled", "UpsilonP");  
243   genupsilonP->SetPtRange(0,100.);  
244   genupsilonP->SetYRange(-8.,8);
245   genupsilonP->SetPhiRange(0.,360.);
246   genupsilonP->SetForceDecay(kDiMuon);
247   genupsilonP->SetTrackingFlag(1);
248   Double_t ratioupsilonP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
249   Double_t sigmaupsilonP     = 0.246e-6 * beautoniumshadowing;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
250   Double_t brupsilonP        = 0.0131;  // Branching Ratio for UpsilonP
251   genupsilonP->Init();  // Generating pT and Y parametrsation for the 4pi
252   ratioupsilonP = sigmaupsilonP * brupsilonP * fNumberOfCollisions / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
253   AliInfo(Form("Upsilon 2S production cross-section in pp with shadowing %5.3g barns",sigmaupsilonP));
254   AliInfo(Form("Upsilon 2S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilonP));
255   genupsilonP->SetPtRange(ptMin, ptMax);  
256   genupsilonP->SetYRange(yMin, yMax);
257   genupsilonP->SetPhiRange(phiMin, phiMax);
258   genupsilonP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
259   AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP);
260   fTotalRate+=ratioupsilonP;
261
262   // Generating UpsilonPP Physics
263   // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
264   AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF Scaled", "UpsilonPP");  
265   genupsilonPP->SetPtRange(0,100.);  
266   genupsilonPP->SetYRange(-8.,8);
267   genupsilonPP->SetPhiRange(0.,360.);
268   genupsilonPP->SetForceDecay(kDiMuon);
269   genupsilonPP->SetTrackingFlag(1);
270   Double_t ratioupsilonPP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
271   Double_t sigmaupsilonPP     = 0.100e-6 * beautoniumshadowing;  //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
272   Double_t brupsilonPP        = 0.0181;  // Branching Ratio for UpsilonPP
273   genupsilonPP->Init();  // Generating pT and Y parametrsation for the 4pi
274   ratioupsilonPP = sigmaupsilonPP * brupsilonPP * fNumberOfCollisions / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
275   AliInfo(Form("Upsilon 3S production cross-section in pp with shadowing %5.3g barns",sigmaupsilonPP));
276   AliInfo(Form("Upsilon 3S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilonPP));
277   genupsilonPP->SetPtRange(ptMin, ptMax);  
278   genupsilonPP->SetYRange(yMin, yMax);
279   genupsilonPP->SetPhiRange(phiMin, phiMax);
280   genupsilonPP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
281   AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP);
282   fTotalRate+=ratioupsilonPP;
283
284 // Generating non-correlated Charm Physics 
285   AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "pp", "Charm");  
286   gencharm->SetPtRange(0,100.);  
287   gencharm->SetYRange(-8.,8);
288   gencharm->SetPhiRange(0.,360.);
289   gencharm->SetForceDecay(kSemiMuonic);
290   gencharm->SetTrackingFlag(1);
291   Double_t ratiocharm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
292   Double_t sigmacharm     = 2. * 6.64e-3 * charmshadowing ;   // 
293   Double_t brcharm        = 0.12;  // Branching Ratio for Charm
294   gencharm->Init();  // Generating pT and Y parametrsation for the 4pi
295   ratiocharm = sigmacharm * brcharm * fNumberOfCollisions / sigmaReaction * gencharm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
296   AliInfo(Form("Charm production cross-section in pp with shadowing %5.3g barns",sigmacharm));
297   AliInfo(Form("Charm production probability per collisions in acceptance via the semi-muonic channel %5.3g",ratiocharm));
298   gencharm->SetPtRange(ptMin, ptMax);  
299   gencharm->SetYRange(yMin, yMax);
300   gencharm->SetPhiRange(phiMin, phiMax);
301   gencharm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
302   AddGenerator(gencharm,"Charm", ratiocharm);
303   fTotalRate+=ratiocharm;
304
305 // Generating non-correlated Beauty Physics 
306   AliGenParam * genbeauty = new AliGenParam(1, AliGenMUONlib::kBeauty, "pp", "Beauty");  
307   genbeauty->SetPtRange(0,100.);  
308   genbeauty->SetYRange(-8.,8);
309   genbeauty->SetPhiRange(0.,360.);
310   genbeauty->SetForceDecay(kSemiMuonic);
311   genbeauty->SetTrackingFlag(1);
312   Double_t ratiobeauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
313   Double_t sigmabeauty     = 2. * 0.210e-3 * beautyshadowing;   // 
314   Double_t brbeauty        = 0.15;  // Branching Ratio for Beauty
315   genbeauty->Init();  // Generating pT and Y parametrsation for the 4pi
316   ratiobeauty = sigmabeauty * brbeauty * fNumberOfCollisions / sigmaReaction * 
317     genbeauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
318   AliInfo(Form("Beauty production cross-section in pp with shadowing %5.3g barns",sigmabeauty));
319   AliInfo(Form("Beauty production probability per collisions in acceptance via the semi-muonic channel %5.3g",ratiobeauty));
320   genbeauty->SetPtRange(ptMin, ptMax);  
321   genbeauty->SetYRange(yMin, yMax);
322   genbeauty->SetPhiRange(phiMin, phiMax);
323   genbeauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
324   AddGenerator(genbeauty,"Beauty", ratiobeauty);
325   fTotalRate+=ratiobeauty;
326
327   // Only if hadronic muons are included in the cocktail
328   if(fHadronicMuons) { 
329     // Generating Pion Physics
330     // The scaling with Npart and Ncoll has been obtained to reproduced tha values presented by Valeri lors de presentatation 
331     // a Clermont Ferrand http://pcrochet.home.cern.ch/pcrochet/files/valerie.pdf
332     //  b range(fm)   Ncoll Npart N_mu pT>0.4 GeV/c
333     //    0 -  3      1982   381    3.62
334     //    3 -  6      1388   297    3.07
335     //    6 -  9       674   177    1.76
336     //    9 - 12       188    71    0.655
337     //   12 - 16        15    10    0.086
338     // We found the hadronic muons scales quite well with the number of participants  
339     AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "default", "Pion");  
340     genpion->SetPtRange(0,100.);  
341     genpion->SetYRange(-8.,8);
342     genpion->SetPhiRange(0.,360.);
343     genpion->SetForceDecay(kPiToMu);
344     genpion->SetTrackingFlag(1);
345     Double_t ratiopion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
346     Double_t sigmapion     = 1.80e-2; // Just for reproducing Valeries's data
347     Double_t brpion        = 0.9999;  // Branching Ratio for Pion 
348     genpion->Init();  // Generating pT and Y parametrsation for the 4pi
349     ratiopion = sigmapion * brpion * (0.93*fNumberOfParticipants+0.07*fNumberOfCollisions) / sigmaReaction * genpion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
350     AliInfo(Form("Pseudo-Pion production cross-section in pp with shadowing %5.3g barns",sigmapion));
351     AliInfo(Form("Pion production probability per collisions in acceptance via the muonic channel %5.3g",ratiopion));
352     genpion->SetPtRange(ptMin, ptMax);  
353     genpion->SetYRange(yMin, yMax);
354     genpion->SetPhiRange(phiMin, phiMax);
355     genpion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
356     AddGenerator(genpion,"Pion", ratiopion);
357     fTotalRate+=ratiopion;
358     
359     // Generating Kaon Physics
360     AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "default", "Kaon");  
361     genkaon->SetPtRange(0,100.);  
362     genkaon->SetYRange(-8.,8);
363     genkaon->SetPhiRange(0.,360.);
364     genkaon->SetForceDecay(kKaToMu);
365     genkaon->SetTrackingFlag(1);
366     Double_t ratiokaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
367     Double_t sigmakaon     = 2.40e-4;   // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06 
368     Double_t brkaon        = 0.6351 ;  // Branching Ratio for Kaon 
369     genkaon->Init();  // Generating pT and Y parametrsation for the 4pi
370     ratiokaon = sigmakaon * brkaon * (0.93*fNumberOfParticipants+0.07*fNumberOfCollisions)/ sigmaReaction * genkaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
371     AliInfo(Form("Pseudo-kaon production cross-section in pp with shadowing %5.3g barns",sigmakaon));
372     AliInfo(Form("Kaon production probability per collisions in acceptance via the muonic channel %5.3g",ratiokaon));
373     genkaon->SetPtRange(ptMin, ptMax);  
374     genkaon->SetYRange(yMin, yMax);
375     genkaon->SetPhiRange(phiMin, phiMax);
376     genkaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
377     AddGenerator(genkaon,"Kaon", ratiokaon);
378     fTotalRate+=ratiokaon;
379   }
380 }
381
382 //_________________________________________________________________________
383 void AliGenMUONCocktail::Generate()
384 {
385   //
386 // Generate event 
387     TIter next(fEntries);
388     AliGenCocktailEntry *entry = 0;
389     AliGenCocktailEntry *preventry = 0;
390     AliGenerator* gen = 0;
391     TObjArray *partArray = gAlice->GetMCApp()->Particles();
392
393 //  Generate the vertex position used by all generators    
394     if(fVertexSmear == kPerEvent) Vertex();
395
396     // Loop on primordialTrigger
397     Bool_t primordialTrigger = kFALSE;
398     while(!primordialTrigger) {
399                 //Reseting stack
400                 AliRunLoader * runloader = gAlice->GetRunLoader();
401                 if (runloader)
402                         if (runloader->Stack())
403                                 runloader->Stack()->Reset();
404                 // Loop over generators and generate events
405                 Int_t igen=0;
406                 Int_t npart =0;
407                 while((entry = (AliGenCocktailEntry*)next())) {
408                         gen = entry->Generator();
409                         gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
410                         if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
411                                 igen++; 
412                                 if (igen == 1) entry->SetFirst(0);
413                                 else  entry->SetFirst((partArray->GetEntriesFast())+1);
414                                 // if ( (fHadronicMuons == kFALSE) && ( (gen->GetName() == "Pions") || (gen->GetName() == "Kaons") ) )
415                                 //  { AliInfo(Form("This generator %s is finally not generated. This is option for hadronic muons.",gen->GetName() ) ); }
416                                 // else {
417                                 gen->SetNumberParticles(npart);
418                                 gen->Generate();
419                                 entry->SetLast(partArray->GetEntriesFast());
420                                 preventry = entry;
421                           // }
422                         }
423                 }  
424                 next.Reset();
425                 // Testing primordial trigger : Muon  pair in the MUON spectrometer acceptance and pTCut
426                 Int_t iPart;
427                 fNGenerated++;
428                 Int_t numberOfMuons=0;
429                 //      printf(">>>fNGenerated is %d\n",fNGenerated);
430                 
431                 TObjArray GoodMuons; // Used in the Invariant Mass selection cut
432                 
433                 for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){      
434                         
435                         if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13)  &&
436                         (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
437                         (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
438                         (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut)                             ) { 
439                                         gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);   
440                                         GoodMuons.AddLast(gAlice->GetMCApp()->Particle(iPart));
441                                         numberOfMuons++;
442                         }                       
443                 }
444                 
445                 // Test the invariant mass of each pair (if cut on Invariant mass is required)
446                 Bool_t InvMassRangeOK = kTRUE;
447                 if(fInvMassCut && (numberOfMuons>=2) ){
448                         TLorentzVector fV1, fV2, fVtot;
449                         InvMassRangeOK = kFALSE;
450                         for(iPart=0; iPart<GoodMuons.GetEntriesFast(); iPart++){      
451                         TParticle * mu1 = ((TParticle *)GoodMuons.At(iPart));
452
453                                 for(int iPart2=iPart+1; iPart2<GoodMuons.GetEntriesFast(); iPart2++){      
454                                         TParticle * mu2 = ((TParticle *)GoodMuons.At(iPart2));
455                                                 
456                                                 fV1.SetPxPyPzE(mu1->Px() ,mu1->Py() ,mu1->Pz() ,mu1->Energy() );
457                                                 fV2.SetPxPyPzE(mu2->Px() ,mu2->Py() ,mu2->Pz() ,mu2->Energy() );
458                                                 fVtot = fV1 + fV2;
459                                                 
460                                                 if(fVtot.M()>fInvMassMinCut && fVtot.M()<fInvMassMaxCut) {
461                                                         InvMassRangeOK = kTRUE;
462                                                         break;
463                                                 }
464                                 }
465                                 if(InvMassRangeOK) break; // Invariant Mass Cut pass as soon as one pair satisfy the criterion
466                         }       
467                 }
468                 
469                 
470                 if ((numberOfMuons >= fMuonMultiplicity) &&  InvMassRangeOK ) primordialTrigger = kTRUE;
471         }
472     fNSucceded++;
473
474     AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
475 }
476
477
478
479
480
481