Introducing the interaction time into the aliroot generators. In case of gaussian...
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
CommitLineData
103ac317 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//
ba762ab4 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
103ac317 22// The free parameeters are :
23// pp reaction cross-section
24// production cross-sections in pp collisions
b44c3901 25// July 07:added heavy quark production from AliGenCorrHF and heavy quark
26// production switched off in Pythia
f81a4aca 27// Aug. 07: added trigger cut on total momentum
00e6c5ee 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)
ba8bf3f5 31//-----------------
00e6c5ee 32// 2009: added polarization (L. Bianchi)
ba8bf3f5 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
ba762ab4 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
103ac317 46#include <TObjArray.h>
47#include <TParticle.h>
48#include <TF1.h>
aaa95f22 49#include <TVirtualMC.h>
103ac317 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"
aaa95f22 59#include "AliDecayer.h"
103ac317 60#include "AliLog.h"
b44c3901 61#include "AliGenCorrHF.h"
00e6c5ee 62#include "AliDecayerPolarized.h"
103ac317 63
64ClassImp(AliGenMUONCocktailpp)
65
66//________________________________________________________________________
67AliGenMUONCocktailpp::AliGenMUONCocktailpp()
1c56e311 68 :AliGenCocktail(),
aaa95f22 69 fDecayer(0),
70 fDecayModeResonance(kAll),
71 fDecayModePythia(kAll),
1c56e311 72 fMuonMultiplicity(0),
73 fMuonPtCut(0.),
f81a4aca 74 fMuonPCut(0.),
18bc0c8e 75 fMuonThetaMinCut(0.),
f81a4aca 76 fMuonThetaMaxCut(180.),
aaa95f22 77 fMuonOriginCut(-999.),
18bc0c8e 78 fNSucceded(0),
9ff768ee 79 fNGenerated(0),
00e6c5ee 80
81 fJpsiPol(0),
ba8bf3f5 82 fChic1Pol(0),
83 fChic2Pol(0),
00e6c5ee 84 fPsiPPol(0),
85 fUpsPol(0),
86 fUpsPPol(0),
87 fUpsPPPol(0),
88 fPolFrame(0),
89
ba762ab4 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
ba8bf3f5 117// x-sections for pp @ 7 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
118// bottomnium as for 10 TeV
ba762ab4 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
ba8bf3f5 131//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
132// scaled down according to ccbar and bbar cross-sections
ba762ab4 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
ba8bf3f5 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
ba762ab4 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
103ac317 160}
93a2041b 161
103ac317 162//_________________________________________________________________________
163AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
ba762ab4 164{
103ac317 165// Destructor
ba762ab4 166
167}
168
169//_________________________________________________________________________
170void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy)
103ac317 171{
ba762ab4 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];
103ac317 184}
185
186//_________________________________________________________________________
00e6c5ee 187void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
188 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
ba762ab4 189// setter for resonances polarization
00e6c5ee 190 if (strcmp(PolFrame,"kColSop")==0){
ba762ab4 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;
00e6c5ee 197 } else if (strcmp(PolFrame,"kHelicity")==0){
ba762ab4 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
00e6c5ee 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//_________________________________________________________________________
8058ead1 217void AliGenMUONCocktailpp::CreateCocktail()
103ac317 218{
ba762ab4 219// create and add resonances and open HF to the coctail
220 Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
221
103ac317 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
18bc0c8e 232// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
233// corrected from feed down of higher resonances
103ac317 234
9ff768ee 235 Double_t sigmajpsi = fSigmaJPsi;
ba8bf3f5 236 Double_t sigmachic1 = fSigmaChic1;
237 Double_t sigmachic2 = fSigmaChic2;
9ff768ee 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;
3285b419 244
00e6c5ee 245// Cross sections corrected with the BR in mu+mu-
246// (only in case of use of AliDecayerPolarized)
247
ba762ab4 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;}
00e6c5ee 255
aaa95f22 256 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
257
ba762ab4 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
c478a271 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");
ba762ab4 282
283 if(cmsEnergy == 10){
00e6c5ee 284 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
ba8bf3f5 285 genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 10", "Chic1");
ba8bf3f5 286 genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 10", "Chic2");
00e6c5ee 287 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
00e6c5ee 288 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
00e6c5ee 289
00e6c5ee 290 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
00e6c5ee 291 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
ba762ab4 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
00e6c5ee 310 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
418e091f 311 } else {
312 AliError("Initialisation failed");
313 return;
00e6c5ee 314 }
418e091f 315
00e6c5ee 316
ba762ab4 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);
b44c3901 324
325//------------------------------------------------------------------
326// Generator of charm
ba762ab4 327 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
8058ead1 328 gencharm->SetMomentumRange(0,9999);
329 gencharm->SetForceDecay(kAll);
ba762ab4 330 Double_t ratioccbar = sigmaccbar/fSigmaReaction;
331 if (!gMC) gencharm->SetDecayer(fDecayer);
8058ead1 332 gencharm->Init();
9ff768ee 333 if (!fSigmaSilent) {
00e6c5ee 334 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
3285b419 335 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
336 }
8058ead1 337 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
103ac317 338//------------------------------------------------------------------
b44c3901 339// Generator of beauty
ba762ab4 340 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
8058ead1 341 genbeauty->SetMomentumRange(0,9999);
342 genbeauty->SetForceDecay(kAll);
ba762ab4 343 Double_t ratiobbbar = sigmabbbar/fSigmaReaction;
344 if (!gMC) genbeauty->SetDecayer(fDecayer);
8058ead1 345 genbeauty->Init();
9ff768ee 346 if (!fSigmaSilent) {
00e6c5ee 347 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
3285b419 348 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
349 }
ba762ab4 350 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
b44c3901 351
8058ead1 352//-------------------------------------------------------------------
103ac317 353// Pythia generator
8058ead1 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);
ba8bf3f5 371
8058ead1 372}
103ac317 373
ba762ab4 374//-------------------------------------------------------------------
375void 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 }
418e091f 415 if (decReso) {
416 decReso->SetForceDecay(kAll);
417 decReso->Init();
418 genReso->SetDecayer(decReso);
419 }
ba762ab4 420 }
421
422 AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
423}
424
425//-------------------------------------------------------------------
8058ead1 426void AliGenMUONCocktailpp::Init()
427{
8058ead1 428 // Initialisation
aaa95f22 429 TIter next(fEntries);
430 AliGenCocktailEntry *entry;
431 if (fStack) {
432 while((entry = (AliGenCocktailEntry*)next())) {
433 entry->Generator()->SetStack(fStack);
434 }
435 }
103ac317 436}
437
438//_________________________________________________________________________
439void AliGenMUONCocktailpp::Generate()
440{
441// Generate event
442 TIter next(fEntries);
443 AliGenCocktailEntry *entry = 0;
444 AliGenCocktailEntry *preventry = 0;
445 AliGenerator* gen = 0;
446
01c89ce1 447 if (fHeader) delete fHeader;
448 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
449
d94af0c1 450 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
103ac317 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
33c3c91a 461 AliRunLoader * runloader = AliRunLoader::Instance();
103ac317 462 if (runloader)
463 if (runloader->Stack())
34da8494 464 runloader->Stack()->Clean();
103ac317 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));
21391258 473 gen->SetTime(fTime);
103ac317 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
ba762ab4 491
103ac317 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++;
b44c3901 496 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22 497 for(iPart=0; iPart<maxPart; iPart++){
103ac317 498
aaa95f22 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) &&
f81a4aca 504 (part->Pt()>fMuonPtCut) &&
505 (part->P()>fMuonPCut)) {
aaa95f22 506 numberOfMuons++;
103ac317 507 }
aaa95f22 508 }
509 }
01c89ce1 510 if (numberOfMuons >= fMuonMultiplicity) {
511 primordialTrigger = kTRUE;
512 fHeader->SetNProduced(maxPart);
513 }
514
103ac317 515 }
516 fNSucceded++;
517
01c89ce1 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);
21391258 523 fHeader->SetInteractionTime(fTime);
01c89ce1 524
525 gAlice->SetGenEventHeader(fHeader);
526
aaa95f22 527// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 528 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
529}
530
531