coverity fixes
[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
275 sprintf(nameJpsi,"Jpsi");
276 sprintf(nameChic1,"Chic1");
277 sprintf(nameChic2,"Chic2");
278 sprintf(namePsiP,"PsiP");
279 sprintf(nameUps,"Ups");
280 sprintf(nameUpsP,"UpsP");
281 sprintf(nameUpsPP,"UpsPP");
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");
311 }
312
ba762ab4 313 AddReso2Generator(nameJpsi,genjpsi,sigmajpsi,fJpsiPol);
314 AddReso2Generator(nameChic1,genchic2,sigmachic1,fChic2Pol);
315 AddReso2Generator(nameChic2,genpsiP,sigmapsiP,fPsiPPol);
316 AddReso2Generator(namePsiP,genchic1,sigmachic1,fChic1Pol);
317 AddReso2Generator(nameUps,genupsilon,sigmaupsilon,fUpsPol);
318 AddReso2Generator(nameUpsP,genupsilonP,sigmaupsilonP,fUpsPPol);
319 AddReso2Generator(nameUpsPP,genupsilonPP,sigmaupsilonPP,fUpsPPPol);
b44c3901 320
321//------------------------------------------------------------------
322// Generator of charm
ba762ab4 323 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
8058ead1 324 gencharm->SetMomentumRange(0,9999);
325 gencharm->SetForceDecay(kAll);
ba762ab4 326 Double_t ratioccbar = sigmaccbar/fSigmaReaction;
327 if (!gMC) gencharm->SetDecayer(fDecayer);
8058ead1 328 gencharm->Init();
9ff768ee 329 if (!fSigmaSilent) {
00e6c5ee 330 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
3285b419 331 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
332 }
8058ead1 333 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
103ac317 334//------------------------------------------------------------------
b44c3901 335// Generator of beauty
ba762ab4 336 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
8058ead1 337 genbeauty->SetMomentumRange(0,9999);
338 genbeauty->SetForceDecay(kAll);
ba762ab4 339 Double_t ratiobbbar = sigmabbbar/fSigmaReaction;
340 if (!gMC) genbeauty->SetDecayer(fDecayer);
8058ead1 341 genbeauty->Init();
9ff768ee 342 if (!fSigmaSilent) {
00e6c5ee 343 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
3285b419 344 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
345 }
ba762ab4 346 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
b44c3901 347
8058ead1 348//-------------------------------------------------------------------
103ac317 349// Pythia generator
8058ead1 350//
351// This has to go into the Config.C
352//
353// AliGenPythia *pythia = new AliGenPythia(1);
354// pythia->SetProcess(kPyMbMSEL1);
355// pythia->SetStrucFunc(kCTEQ5L);
356// pythia->SetEnergyCMS(14000.);
357// AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
358// Decay_t dt = gener->GetDecayModePythia();
359// pythia->SetForceDecay(dt);
360// pythia->SetPtRange(0.,100.);
361// pythia->SetYRange(-8.,8.);
362// pythia->SetPhiRange(0.,360.);
363// pythia->SetPtHard(2.76,-1.0);
364// pythia->SwitchHFOff();
365// pythia->Init();
366// AddGenerator(pythia,"Pythia",1);
ba8bf3f5 367
8058ead1 368}
103ac317 369
ba762ab4 370//-------------------------------------------------------------------
371void AliGenMUONCocktailpp::AddReso2Generator(Char_t* nameReso,
372 AliGenParam* const genReso,
373 Double_t sigmaReso,
374 Double_t polReso)
375{
376// add resonances to the cocktail
377 Double_t phiMin = fPhiMin*180./TMath::Pi();
378 Double_t phiMax = fPhiMax*180./TMath::Pi();
379
380 // first step: generation in 4pi
381 genReso->SetPtRange(0.,100.);
382 genReso->SetYRange(-8.,8.);
383 genReso->SetPhiRange(0.,360.);
384 genReso->SetForceDecay(fDecayModeResonance);
385 if (!gMC) genReso->SetDecayer(fDecayer);
386 genReso->Init(); // generation in 4pi
387// Ratios with respect to the reaction cross-section in the
388// kinematics limit of the MUONCocktail
389 Double_t ratioReso = sigmaReso / fSigmaReaction * genReso->GetRelativeArea(fPtMin,fPtMax,fYMin,fYMax,phiMin,phiMax);
390 if (!fSigmaSilent) {
391 AliInfo(Form("%s prod. cross-section in pp %5.3g b",nameReso,sigmaReso));
392 AliInfo(Form("%s prod. probability per collision in acceptance %5.3g",nameReso,ratioReso));
393 }
394// second step: generation in selected kinematical range
395 genReso->SetPtRange(fPtMin, fPtMax);
396 genReso->SetYRange(fYMin, fYMax);
397 genReso->SetPhiRange(phiMin, phiMax);
398 genReso->Init(); // generation in selected kinematical range
399
400 TVirtualMCDecayer *decReso = 0;
401 if(TMath::Abs(polReso) > 1.e-30){
402 AliInfo(Form("******Setting polarized decayer for %s''",nameReso));
403 if(fPolFrame==0) {
404 decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
405 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",polReso));
406 }
407 if(fPolFrame==1) {
408 decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
409 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",polReso));
410 }
411 decReso->SetForceDecay(kAll);
412 decReso->Init();
413 genReso->SetDecayer(decReso);
414 }
415
416 AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
417}
418
419//-------------------------------------------------------------------
8058ead1 420void AliGenMUONCocktailpp::Init()
421{
8058ead1 422 // Initialisation
aaa95f22 423 TIter next(fEntries);
424 AliGenCocktailEntry *entry;
425 if (fStack) {
426 while((entry = (AliGenCocktailEntry*)next())) {
427 entry->Generator()->SetStack(fStack);
428 }
429 }
103ac317 430}
431
432//_________________________________________________________________________
433void AliGenMUONCocktailpp::Generate()
434{
435// Generate event
436 TIter next(fEntries);
437 AliGenCocktailEntry *entry = 0;
438 AliGenCocktailEntry *preventry = 0;
439 AliGenerator* gen = 0;
440
01c89ce1 441 if (fHeader) delete fHeader;
442 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
443
d94af0c1 444 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
103ac317 445
446// Generate the vertex position used by all generators
447 if(fVertexSmear == kPerEvent) Vertex();
448
449// Loop on primordialTrigger:
450// minimum muon multiplicity above a pt cut in a theta acceptance region
451
452 Bool_t primordialTrigger = kFALSE;
453 while(!primordialTrigger) {
454 //Reseting stack
33c3c91a 455 AliRunLoader * runloader = AliRunLoader::Instance();
103ac317 456 if (runloader)
457 if (runloader->Stack())
34da8494 458 runloader->Stack()->Clean();
103ac317 459 // Loop over generators and generate events
460 Int_t igen = 0;
461 Int_t npart = 0;
462 const char* genName = 0;
463 while((entry = (AliGenCocktailEntry*)next())) {
464 gen = entry->Generator();
465 genName = entry->GetName();
466 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
467
468 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
469 gRandom->Poisson(entry->Rate());
470
471 if (npart > 0) {
472 igen++;
473 if (igen == 1) entry->SetFirst(0);
474 else entry->SetFirst((partArray->GetEntriesFast())+1);
475
476 gen->SetNumberParticles(npart);
477 gen->Generate();
478 entry->SetLast(partArray->GetEntriesFast());
479 preventry = entry;
480 }
481 }
482 next.Reset();
483
ba762ab4 484
103ac317 485// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
486// in the muon spectrometer acceptance
487 Int_t iPart;
488 fNGenerated++;
b44c3901 489 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22 490 for(iPart=0; iPart<maxPart; iPart++){
103ac317 491
aaa95f22 492 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
493 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
494 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
495 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
496 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
f81a4aca 497 (part->Pt()>fMuonPtCut) &&
498 (part->P()>fMuonPCut)) {
aaa95f22 499 numberOfMuons++;
103ac317 500 }
aaa95f22 501 }
502 }
01c89ce1 503 if (numberOfMuons >= fMuonMultiplicity) {
504 primordialTrigger = kTRUE;
505 fHeader->SetNProduced(maxPart);
506 }
507
103ac317 508 }
509 fNSucceded++;
510
01c89ce1 511 TArrayF eventVertex;
512 eventVertex.Set(3);
513 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
514
515 fHeader->SetPrimaryVertex(eventVertex);
516
517 gAlice->SetGenEventHeader(fHeader);
518
aaa95f22 519// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 520 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
521}
522
523