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