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