New class added.
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.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 // Class to create the coktail for physics with muons for pp collisions
20 // using the The followoing sources: 
21 // jpsi, psiP, upsilon, upsilonP, upsilonPP, open charm and open beauty
22 // The free parameeters are :
23 //      pp reaction cross-section
24 //      production cross-sections in pp collisions 
25 // July 07:added heavy quark production from AliGenCorrHF and heavy quark 
26 //         production switched off in Pythia 
27 // Aug. 07: added trigger cut on total momentum
28 // 2009: added possibility to hide x-sections (B. Vulpescu)
29 // 2009: added possibility to have the cocktail (fast generator and param.) 
30 // for pp @ 10 TeV or pp @ 14 TeV (N. Bastid)
31 //-----------------
32 // 2009:  added polarization (L. Bianchi)
33 //------------------
34 // 11/2009: added chi_c1 & chi_c2 (P.Crochet & N.Bastid). 
35 // Cross-sections for charmonia are now directly taken from the Yellow Report 
36 // (hep-ph/0311048) Tab.9, page 19. See below for details w.r.t. beam energy. 
37 // usage: see example of Config in $ALICE_ROOT/prod/LHC09a10/Config.C
38 //------------------------
39 // 04/2010:
40 // - CMS energy passed via parameter
41 // i.e. gener->SetCMSEnergy(AliGenMUONCocktailpp::kCMS07TeV) in Config.C 
42 // - resonances now added to the cocktail via AddReso2Generator 
43 // - cleaning 
44 // B.Vulpescu & P.Crochet         
45  
46 #include <TObjArray.h>
47 #include <TParticle.h>
48 #include <TF1.h>
49 #include <TVirtualMC.h>
50 #include "AliGenCocktailEventHeader.h"
51
52 #include "AliGenCocktailEntry.h"
53 #include "AliGenMUONCocktailpp.h"
54 #include "AliGenMUONlib.h"
55 #include "AliGenParam.h"
56 #include "AliMC.h"
57 #include "AliRun.h"
58 #include "AliStack.h"
59 #include "AliDecayer.h"
60 #include "AliLog.h"
61 #include "AliGenCorrHF.h"
62 #include "AliDecayerPolarized.h"
63
64 ClassImp(AliGenMUONCocktailpp)  
65   
66 //________________________________________________________________________
67 AliGenMUONCocktailpp::AliGenMUONCocktailpp()
68     :AliGenCocktail(),
69      fDecayer(0),
70      fDecayModeResonance(kAll),
71      fDecayModePythia(kAll),
72      fMuonMultiplicity(0),
73      fMuonPtCut(0.),
74      fMuonPCut(0.),
75      fMuonThetaMinCut(0.),
76      fMuonThetaMaxCut(180.),
77      fMuonOriginCut(-999.),
78      fNSucceded(0),
79      fNGenerated(0),
80
81      fJpsiPol(0), 
82      fChic1Pol(0), 
83      fChic2Pol(0), 
84      fPsiPPol(0), 
85      fUpsPol(0), 
86      fUpsPPol(0), 
87      fUpsPPPol(0),
88      fPolFrame(0),
89
90      fCMSEnergyTeV(0),
91      fCMSEnergyTeVArray(),
92      fSigmaReaction(0),  
93      fSigmaReactionArray(),  
94      fSigmaJPsi(0),      
95      fSigmaJPsiArray(),      
96      fSigmaChic1(0),     
97      fSigmaChic1Array(),     
98      fSigmaChic2(0),     
99      fSigmaChic2Array(),     
100      fSigmaPsiP(0),      
101      fSigmaPsiPArray(),      
102      fSigmaUpsilon(0),   
103      fSigmaUpsilonArray(),   
104      fSigmaUpsilonP(0),  
105      fSigmaUpsilonPArray(),  
106      fSigmaUpsilonPP(0), 
107      fSigmaUpsilonPPArray(), 
108      fSigmaCCbar(0),     
109      fSigmaCCbarArray(),     
110      fSigmaBBbar(0),
111      fSigmaBBbarArray(),
112
113      fSigmaSilent(kFALSE)
114 {
115 // Constructor
116
117 // x-sections for pp @ 7 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
118 // bottomnium as for 10 TeV
119     fCMSEnergyTeVArray[0] =   7.00;
120     fSigmaReactionArray[0] =  0.0695;
121     fSigmaJPsiArray[0] =      21.8e-6;
122     fSigmaChic1Array[0] =     21.1e-6;
123     fSigmaChic2Array[0] =     34.9e-6;
124     fSigmaPsiPArray[0] =      4.93e-6;
125     fSigmaUpsilonArray[0] =   0.463e-6;
126     fSigmaUpsilonPArray[0] =  0.154e-6;
127     fSigmaUpsilonPPArray[0] = 0.0886e-6;
128     fSigmaCCbarArray[0] =     6.91e-3;
129     fSigmaBBbarArray[0] =     0.232e-3;
130     
131 //x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
132 // scaled down according to ccbar and bbar cross-sections
133     fCMSEnergyTeVArray[1] =   10.00;
134     fSigmaReactionArray[1] =  0.0695;
135     fSigmaJPsiArray[1] =      26.06e-6;
136     fSigmaChic1Array[1] =     25.18e-6;
137     fSigmaChic2Array[1] =     41.58e-6;
138     fSigmaPsiPArray[1] =      5.88e-6;
139     fSigmaUpsilonArray[1] =   0.658e-6;
140     fSigmaUpsilonPArray[1] =  0.218e-6;
141     fSigmaUpsilonPPArray[1] = 0.122e-6;
142     fSigmaCCbarArray[1] =     8.9e-3;
143     fSigmaBBbarArray[1] =     0.33e-3;
144     
145 //x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
146 // bottomonium from hep-ph/0311048 Tab.9, page 19 taken inton account that 
147 // feed-down from chib is included
148     fCMSEnergyTeVArray[2] =   14.00;
149     fSigmaReactionArray[2] =  0.070;
150     fSigmaJPsiArray[2] =      32.9e-6;
151     fSigmaChic1Array[2] =     31.8e-6;
152     fSigmaChic2Array[2] =     52.5e-6;
153     fSigmaPsiPArray[2] =      7.43e-6;
154     fSigmaUpsilonArray[2] =   0.989e-6;
155     fSigmaUpsilonPArray[2] =  0.502e-6;
156     fSigmaUpsilonPPArray[2] = 0.228e-6;
157     fSigmaCCbarArray[2] =     11.2e-3;
158     fSigmaBBbarArray[2] =     0.51e-3;
159     
160 }
161
162 //_________________________________________________________________________
163 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
164 {
165 // Destructor
166
167 }
168
169 //_________________________________________________________________________
170 void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy) 
171 {
172 // setter for CMSEnergy and corresponding cross-sections
173   fCMSEnergyTeV   = fCMSEnergyTeVArray[cmsEnergy];
174   fSigmaReaction  = fSigmaReactionArray[cmsEnergy];  
175   fSigmaJPsi      = fSigmaJPsiArray[cmsEnergy];      
176   fSigmaChic1     = fSigmaChic1Array[cmsEnergy];     
177   fSigmaChic2     = fSigmaChic2Array[cmsEnergy];     
178   fSigmaPsiP      = fSigmaPsiPArray[cmsEnergy];      
179   fSigmaUpsilon   = fSigmaUpsilonArray[cmsEnergy];   
180   fSigmaUpsilonP  = fSigmaUpsilonPArray[cmsEnergy];  
181   fSigmaUpsilonPP = fSigmaUpsilonPPArray[cmsEnergy]; 
182   fSigmaCCbar     = fSigmaCCbarArray[cmsEnergy];     
183   fSigmaBBbar     = fSigmaBBbarArray[cmsEnergy];
184 }
185
186 //_________________________________________________________________________
187 void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
188                                 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
189 // setter for resonances polarization   
190    if (strcmp(PolFrame,"kColSop")==0){
191        fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
192        fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
193        fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
194        fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
195        fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
196        fPolFrame = 0;
197    } else if (strcmp(PolFrame,"kHelicity")==0){
198        fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
199        fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
200        fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
201        fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
202        fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
203        fPolFrame = 1;
204
205    } else {
206        AliInfo(Form("The polarization frame is not valid"));
207        AliInfo(Form("No polarization will be set"));
208        fJpsiPol=0.;
209        fPsiPPol=0.;
210        fUpsPol=0.;
211        fUpsPPol=0.;
212        fUpsPPPol=0.;
213    }
214 }
215
216 //_________________________________________________________________________
217 void AliGenMUONCocktailpp::CreateCocktail()
218 {
219 // create and add resonances and open HF to the coctail
220     Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
221
222 // These limits are only used for renormalization of quarkonia cross section
223 // Pythia events are generated in 4pi  
224     Double_t ptMin  = fPtMin;
225     Double_t ptMax  = fPtMax;
226     Double_t yMin   = fYMin;;
227     Double_t yMax   = fYMax;;
228     Double_t phiMin = fPhiMin*180./TMath::Pi();
229     Double_t phiMax = fPhiMax*180./TMath::Pi();
230     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));
231     
232 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and 
233 // corrected from feed down of higher resonances 
234
235     Double_t sigmajpsi      = fSigmaJPsi;  
236     Double_t sigmachic1     = fSigmaChic1;  
237     Double_t sigmachic2     = fSigmaChic2;  
238     Double_t sigmapsiP      = fSigmaPsiP;  
239     Double_t sigmaupsilon   = fSigmaUpsilon;  
240     Double_t sigmaupsilonP  = fSigmaUpsilonP;  
241     Double_t sigmaupsilonPP = fSigmaUpsilonPP;
242     Double_t sigmaccbar     = fSigmaCCbar;
243     Double_t sigmabbbar     = fSigmaBBbar;
244
245 // Cross sections corrected with the BR in mu+mu-
246 // (only in case of use of AliDecayerPolarized)
247
248     if(TMath::Abs(fJpsiPol) > 1.e-30) {sigmajpsi      = fSigmaJPsi*0.0593;}
249     if(TMath::Abs(fChic1Pol) > 1.e-30) {sigmachic1     = fSigmaChic1*0.;} // tb consistent
250     if(TMath::Abs(fChic2Pol) > 1.e-30) {sigmachic2     = fSigmaChic2*0.;} // tb consistent 
251     if(TMath::Abs(fPsiPPol) > 1.e-30) {sigmapsiP      = fSigmaPsiP*0.0075;}
252     if(TMath::Abs(fUpsPol) > 1.e-30) {sigmaupsilon   = fSigmaUpsilon*0.0248;}
253     if(TMath::Abs(fUpsPPol) > 1.e-30) {sigmaupsilonP  = fSigmaUpsilonP*0.0193;}
254     if(TMath::Abs(fUpsPPPol) > 1.e-30) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
255
256     AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
257
258 // Create and add resonances to the generator
259     AliGenParam * genjpsi=0;
260     AliGenParam * genchic1=0;
261     AliGenParam * genchic2=0;
262     AliGenParam * genpsiP=0;
263     AliGenParam * genupsilon=0;
264     AliGenParam * genupsilonP=0;
265     AliGenParam * genupsilonPP=0;
266     
267     Char_t nameJpsi[10];    
268     Char_t nameChic1[10];    
269     Char_t nameChic2[10];    
270     Char_t namePsiP[10];
271     Char_t nameUps[10];    
272     Char_t nameUpsP[10];    
273     Char_t nameUpsPP[10];    
274
275     snprintf(nameJpsi,10, "Jpsi");    
276     snprintf(nameChic1,10, "Chic1");
277     snprintf(nameChic2,10, "Chic2");
278     snprintf(namePsiP,10, "PsiP");
279     snprintf(nameUps,10, "Ups");
280     snprintf(nameUpsP,10, "UpsP");
281     snprintf(nameUpsPP,10, "UpsPP");
282
283     if(cmsEnergy == 10){
284         genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
285         genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 10", "Chic1");
286         genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 10", "Chic2");
287         genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
288         genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
289
290         genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
291         genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
292     } else if (cmsEnergy == 7){
293         genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 7", "Jpsi");
294         genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 7", "Chic1");
295         genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 7", "Chic2");
296         genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 7", "PsiP");
297
298         genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 7", "Upsilon");
299         genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 7", "UpsilonP");
300         genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 7", "UpsilonPP");
301     } else if (cmsEnergy == 14){
302         genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi");
303         genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp ", "Chic1");
304         genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp ", "Chic2");
305         genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
306
307         genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
308         genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
309
310         genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");    
311     } else {
312         AliError("Initialisation failed");
313         return;
314     }
315     
316
317     AddReso2Generator(nameJpsi,genjpsi,sigmajpsi,fJpsiPol);
318     AddReso2Generator(nameChic1,genchic2,sigmachic1,fChic2Pol);
319     AddReso2Generator(nameChic2,genpsiP,sigmapsiP,fPsiPPol);    
320     AddReso2Generator(namePsiP,genchic1,sigmachic1,fChic1Pol);    
321     AddReso2Generator(nameUps,genupsilon,sigmaupsilon,fUpsPol);    
322     AddReso2Generator(nameUpsP,genupsilonP,sigmaupsilonP,fUpsPPol);    
323     AddReso2Generator(nameUpsPP,genupsilonPP,sigmaupsilonPP,fUpsPPPol);    
324
325 //------------------------------------------------------------------
326 // Generator of charm
327     AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
328     gencharm->SetMomentumRange(0,9999);
329     gencharm->SetForceDecay(kAll);
330     Double_t ratioccbar = sigmaccbar/fSigmaReaction;
331     if (!gMC) gencharm->SetDecayer(fDecayer);  
332     gencharm->Init();
333     if (!fSigmaSilent) {
334       AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
335       AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
336     }
337     AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
338 //------------------------------------------------------------------
339 // Generator of beauty
340     AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
341     genbeauty->SetMomentumRange(0,9999);
342     genbeauty->SetForceDecay(kAll);
343     Double_t ratiobbbar = sigmabbbar/fSigmaReaction;
344     if (!gMC) genbeauty->SetDecayer(fDecayer);  
345     genbeauty->Init();
346     if (!fSigmaSilent) {
347       AliInfo(Form("b-bbar prod. cross-section in pp  %5.3g b",sigmabbbar));
348       AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
349     }
350     AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
351
352 //-------------------------------------------------------------------
353 // Pythia generator
354 //
355 // This has to go into the Config.C
356 //
357 //    AliGenPythia *pythia = new AliGenPythia(1);
358 //    pythia->SetProcess(kPyMbMSEL1);
359 //    pythia->SetStrucFunc(kCTEQ5L);
360 //    pythia->SetEnergyCMS(14000.);
361 //    AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
362 //    Decay_t dt = gener->GetDecayModePythia();
363 //    pythia->SetForceDecay(dt);
364 //    pythia->SetPtRange(0.,100.);
365 //    pythia->SetYRange(-8.,8.);
366 //    pythia->SetPhiRange(0.,360.);
367 //    pythia->SetPtHard(2.76,-1.0);
368 //    pythia->SwitchHFOff();
369 //    pythia->Init(); 
370 //    AddGenerator(pythia,"Pythia",1);
371
372 }
373
374 //-------------------------------------------------------------------
375 void AliGenMUONCocktailpp::AddReso2Generator(Char_t* nameReso, 
376                                              AliGenParam* const genReso,
377                                              Double_t sigmaReso,
378                                              Double_t polReso)
379 {
380 // add resonances to the cocktail
381     Double_t phiMin = fPhiMin*180./TMath::Pi();
382     Double_t phiMax = fPhiMax*180./TMath::Pi();
383
384     // first step: generation in 4pi
385     genReso->SetPtRange(0.,100.); 
386     genReso->SetYRange(-8.,8.);
387     genReso->SetPhiRange(0.,360.);
388     genReso->SetForceDecay(fDecayModeResonance);
389     if (!gMC) genReso->SetDecayer(fDecayer);  
390     genReso->Init();  // generation in 4pi
391 // Ratios with respect to the reaction cross-section in the 
392 // kinematics limit of the MUONCocktail
393     Double_t ratioReso = sigmaReso / fSigmaReaction * genReso->GetRelativeArea(fPtMin,fPtMax,fYMin,fYMax,phiMin,phiMax);
394     if (!fSigmaSilent) {
395       AliInfo(Form("%s prod. cross-section in pp %5.3g b",nameReso,sigmaReso));
396       AliInfo(Form("%s prod. probability per collision in acceptance %5.3g",nameReso,ratioReso));
397     }
398 // second step: generation in selected kinematical range
399     genReso->SetPtRange(fPtMin, fPtMax);  
400     genReso->SetYRange(fYMin, fYMax);
401     genReso->SetPhiRange(phiMin, phiMax);
402     genReso->Init(); // generation in selected kinematical range
403
404      TVirtualMCDecayer *decReso = 0;
405      if(TMath::Abs(polReso) > 1.e-30){
406       AliInfo(Form("******Setting polarized decayer for %s''",nameReso));
407       if(fPolFrame==0) {
408         decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
409         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",polReso));
410           }
411       if(fPolFrame==1) {
412         decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
413         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",polReso));
414           }
415       if (decReso) {
416           decReso->SetForceDecay(kAll);
417           decReso->Init();
418           genReso->SetDecayer(decReso);
419       }
420      }
421     
422     AddGenerator(genReso,nameReso,ratioReso); // Adding Generator    
423 }
424
425 //-------------------------------------------------------------------
426 void AliGenMUONCocktailpp::Init()
427 {
428     // Initialisation
429     TIter next(fEntries);
430     AliGenCocktailEntry *entry;
431     if (fStack) {
432         while((entry = (AliGenCocktailEntry*)next())) {
433             entry->Generator()->SetStack(fStack);
434         }
435     }
436 }
437
438 //_________________________________________________________________________
439 void AliGenMUONCocktailpp::Generate()
440 {
441 // Generate event 
442     TIter next(fEntries);
443     AliGenCocktailEntry *entry = 0;
444     AliGenCocktailEntry *preventry = 0;
445     AliGenerator* gen = 0;
446
447     if (fHeader) delete fHeader;
448     fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
449
450     const TObjArray *partArray = gAlice->GetMCApp()->Particles();
451     
452 // Generate the vertex position used by all generators    
453     if(fVertexSmear == kPerEvent) Vertex();
454
455 // Loop on primordialTrigger: 
456 // minimum muon multiplicity above a pt cut in a theta acceptance region
457
458     Bool_t primordialTrigger = kFALSE;
459     while(!primordialTrigger) {
460         //Reseting stack
461         AliRunLoader * runloader = AliRunLoader::Instance();
462         if (runloader)
463             if (runloader->Stack())
464                 runloader->Stack()->Clean();
465         // Loop over generators and generate events
466         Int_t igen = 0;
467         Int_t npart = 0;        
468         const char* genName = 0;
469         while((entry = (AliGenCocktailEntry*)next())) {
470             gen = entry->Generator();
471             genName = entry->GetName();
472             gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
473             gen->SetTime(fTime);
474
475             npart = (strcmp(genName,"Pythia") == 0) ? 1 :
476                 gRandom->Poisson(entry->Rate());
477
478             if (npart > 0) {
479                 igen++; 
480                 if (igen == 1) entry->SetFirst(0);              
481                 else  entry->SetFirst((partArray->GetEntriesFast())+1);
482                 
483                 gen->SetNumberParticles(npart);         
484                 gen->Generate();
485                 entry->SetLast(partArray->GetEntriesFast());
486                 preventry = entry;
487             }
488         }  
489         next.Reset();
490
491
492 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut 
493 // in the muon spectrometer acceptance
494         Int_t iPart;
495         fNGenerated++;
496         Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
497         for(iPart=0; iPart<maxPart; iPart++){      
498             
499           TParticle *part = gAlice->GetMCApp()->Particle(iPart);
500           if ( TMath::Abs(part->GetPdgCode()) == 13 ){  
501             if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
502                (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
503                (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
504                (part->Pt()>fMuonPtCut) &&
505                (part->P()>fMuonPCut)) {
506               numberOfMuons++;
507             }
508           }
509         }       
510         if (numberOfMuons >= fMuonMultiplicity) {
511           primordialTrigger = kTRUE;
512           fHeader->SetNProduced(maxPart);
513         }
514
515     }
516     fNSucceded++;
517
518     TArrayF eventVertex;
519     eventVertex.Set(3);
520     for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
521
522     fHeader->SetPrimaryVertex(eventVertex);
523     fHeader->SetInteractionTime(fTime);
524
525     gAlice->SetGenEventHeader(fHeader);
526
527 //     AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
528     AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
529 }
530
531