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