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