Changes needed in PHOS (Y.Schutz)
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktail.cxx
CommitLineData
84954c47 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
9538aedd 18//
84954c47 19// Classe to create the MUON coktail for physics in the Alice muon spectrometer
9538aedd 20// The followoing muons sources are included in this cocktail:
21// jpsi, upsilon, non-correlated open and beauty, and muons from pion and kaons.
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
26// Hard probes are supposed to scale with Ncoll and hadronic production with (0.8Ncoll+0.2*Npart)
27// There is a primordial trigger wiche requires :
28// a minimum muon multiplicity above a pT cut in a theta acceptance cone
29//
30// Gines Martinez, jan 2004, Nantes martinez@in2p3.fr
84954c47 31
32
33//
34
9538aedd 35//#include <TList.h>
84954c47 36#include <TObjArray.h>
37#include <TF1.h>
38#include <TParticle.h>
39
40#include "AliGenParam.h"
41#include "AliGenMUONlib.h"
42#include "AliGenMUONCocktail.h"
43#include "AliGenCocktailEntry.h"
44#include "AliCollisionGeometry.h"
45#include "AliRun.h"
46#include "AliMC.h"
47#include "AliStack.h"
48
49ClassImp(AliGenMUONCocktail)
50
51
52//________________________________________________________________________
53AliGenMUONCocktail::AliGenMUONCocktail()
54 :AliGenCocktail()
55{
56// Constructor
57 fTotalRate =0;
58 fNSucceded=0;
59 fNGenerated=0;
60 fMuonMultiplicity=1;
61 fMuonPtCut= 1.;
62 fMuonThetaMinCut=171.;
63 fMuonThetaMaxCut=178.;
64 fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
65 //
66}
67//_________________________________________________________________________
68AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
69 AliGenCocktail(cocktail)
70{
71// Copy constructor
72 fTotalRate =0;
73 fNSucceded=0;
74 fNGenerated=0;
75 fMuonMultiplicity=1;
76 fMuonPtCut= 1.;
77 fMuonThetaMinCut=171.;
78 fMuonThetaMaxCut=178.;
79 fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
80 //
81}
82//_________________________________________________________________________
83AliGenMUONCocktail::~AliGenMUONCocktail()
84{
85// Destructor
84954c47 86}
87
88//_________________________________________________________________________
89void AliGenMUONCocktail::Init()
90{
91 // Defining MUON physics cocktail
92
93 // Kinematical limits for particle generation
94 Float_t ptMin = fPtMin;
95 Float_t ptMax = fPtMax;
96 Float_t yMin = fYMin;;
97 Float_t yMax = fYMax;;
98 Float_t phiMin = fPhiMin*180./TMath::Pi();
99 Float_t phiMax = fPhiMax*180./TMath::Pi();
100 printf(">>> Kinematical range pT:%f:%f y:%f:%f Phi:%f:%f\n",ptMin,ptMax,yMin,yMax,phiMin,phiMax);
101
9538aedd 102 Float_t sigmaReaction = 0.072; // MINB pp at LHC energies 72 mb
88e5db43 103
84954c47 104 // Generating J/Psi Physics
9538aedd 105 AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi");
84954c47 106 // 4pi Generation
9538aedd 107 genjpsi->SetPtRange(0,100.);
108 genjpsi->SetYRange(-8.,8);
109 genjpsi->SetPhiRange(0.,360.);
110 genjpsi->SetForceDecay(kDiMuon);
111 genjpsi->SetTrackingFlag(1);
84954c47 112 // Calculation of the paritcle multiplicity per event in the muonic channel
9538aedd 113 Float_t ratiojpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
114 Float_t sigmajpsi = 31.0e-6 * 0.437; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing
115 Float_t brjpsi = 0.0588; // Branching Ratio for JPsi
116 genjpsi->Init(); // Generating pT and Y parametrsation for the 4pi
117 ratiojpsi = sigmajpsi * brjpsi * fNumberOfCollisions / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
118 printf(">>> ratio jpsi %g et %g Ncol %g sigma %g\n",ratiojpsi,genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigmajpsi );
84954c47 119 // Generation in the kinematical limits of AliGenMUONCocktail
9538aedd 120 genjpsi->SetPtRange(ptMin, ptMax);
121 genjpsi->SetYRange(yMin, yMax);
122 genjpsi->SetPhiRange(phiMin, phiMax);
123 genjpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
84954c47 124 // Adding Generator
9538aedd 125 AddGenerator(genjpsi, "Jpsi", ratiojpsi);
126 fTotalRate+=ratiojpsi;
84954c47 127
88e5db43 128 // Generating Psi prime Physics
129 AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "Vogt", "PsiP");
130 // 4pi Generation
131 genpsiP->SetPtRange(0,100.);
132 genpsiP->SetYRange(-8.,8);
133 genpsiP->SetPhiRange(0.,360.);
134 genpsiP->SetForceDecay(kDiMuon);
135 genpsiP->SetTrackingFlag(1);
136 // Calculation of the paritcle multiplicity per event in the muonic channel
137 Float_t ratiopsiP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
138 Float_t sigmapsiP = 4.68e-6 * 0.437; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing
139 Float_t brpsiP = 0.0103; // Branching Ratio for PsiP
140 genpsiP->Init(); // Generating pT and Y parametrsation for the 4pi
141 ratiopsiP = sigmapsiP * brpsiP * fNumberOfCollisions / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
142 printf(">>> ratio psiP %g et %g Ncol %g sigma %g\n",ratiopsiP,genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigmapsiP );
143 // Generation in the kinematical limits of AliGenMUONCocktail
144 genpsiP->SetPtRange(ptMin, ptMax);
145 genpsiP->SetYRange(yMin, yMax);
146 genpsiP->SetPhiRange(phiMin, phiMax);
147 genpsiP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
148 // Adding Generator
149 AddGenerator(genpsiP, "PsiP", ratiopsiP);
150 fTotalRate+=ratiopsiP;
151
152
84954c47 153// Generating Upsilon Physics
9538aedd 154 AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon");
155 genupsilon->SetPtRange(0,100.);
156 genupsilon->SetYRange(-8.,8);
157 genupsilon->SetPhiRange(0.,360.);
158 genupsilon->SetForceDecay(kDiMuon);
159 genupsilon->SetTrackingFlag(1);
160 Float_t ratioupsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
161 Float_t sigmaupsilon = 0.501e-6 * 0.674; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing
162 Float_t brupsilon = 0.0248; // Branching Ratio for Upsilon
163 genupsilon->Init(); // Generating pT and Y parametrsation for the 4pi
164 ratioupsilon = sigmaupsilon * brupsilon * fNumberOfCollisions / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
165 printf(">>> ratio upsilon %g et %g\n",ratioupsilon, genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
166 genupsilon->SetPtRange(ptMin, ptMax);
167 genupsilon->SetYRange(yMin, yMax);
168 genupsilon->SetPhiRange(phiMin, phiMax);
169 genupsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
170 AddGenerator(genupsilon,"Upsilon", ratioupsilon);
171 fTotalRate+=ratioupsilon;
84954c47 172
88e5db43 173// Generating UpsilonP Physics
174 AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "Vogt", "UpsilonP");
175 genupsilonP->SetPtRange(0,100.);
176 genupsilonP->SetYRange(-8.,8);
177 genupsilonP->SetPhiRange(0.,360.);
178 genupsilonP->SetForceDecay(kDiMuon);
179 genupsilonP->SetTrackingFlag(1);
180 Float_t ratioupsilonP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
181 Float_t sigmaupsilonP = 0.246e-6 * 0.674; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing
182 Float_t brupsilonP = 0.0131; // Branching Ratio for UpsilonP
183 genupsilonP->Init(); // Generating pT and Y parametrsation for the 4pi
184 ratioupsilonP = sigmaupsilonP * brupsilonP * fNumberOfCollisions / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
185 printf(">>> ratio upsilonP %g et %g\n",ratioupsilonP, genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
186 genupsilonP->SetPtRange(ptMin, ptMax);
187 genupsilonP->SetYRange(yMin, yMax);
188 genupsilonP->SetPhiRange(phiMin, phiMax);
189 genupsilonP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
190 AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP);
191 fTotalRate+=ratioupsilonP;
192
193
194// Generating UpsilonPP Physics
195 AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "Vogt", "UpsilonPP");
196 genupsilonPP->SetPtRange(0,100.);
197 genupsilonPP->SetYRange(-8.,8);
198 genupsilonPP->SetPhiRange(0.,360.);
199 genupsilonPP->SetForceDecay(kDiMuon);
200 genupsilonPP->SetTrackingFlag(1);
201 Float_t ratioupsilonPP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
202 Float_t sigmaupsilonPP = 0.100e-6 * 0.674; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing
203 Float_t brupsilonPP = 0.0181; // Branching Ratio for UpsilonPP
204 genupsilonPP->Init(); // Generating pT and Y parametrsation for the 4pi
205 ratioupsilonPP = sigmaupsilonPP * brupsilonPP * fNumberOfCollisions / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
206 printf(">>> ratio upsilonPP %g et %g\n",ratioupsilonPP, genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
207 genupsilonPP->SetPtRange(ptMin, ptMax);
208 genupsilonPP->SetYRange(yMin, yMax);
209 genupsilonPP->SetPhiRange(phiMin, phiMax);
210 genupsilonPP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
211 AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP);
212 fTotalRate+=ratioupsilonPP;
213
214
84954c47 215// Generating Charm Physics
9538aedd 216 AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm");
217 gencharm->SetPtRange(0,100.);
218 gencharm->SetYRange(-8.,8);
219 gencharm->SetPhiRange(0.,360.);
220 gencharm->SetForceDecay(kSemiMuonic);
221 gencharm->SetTrackingFlag(1);
222 Float_t ratiocharm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
223 Float_t sigmacharm = 2. * 6.64e-3 * 0.65 ; //
224 Float_t brcharm = 0.12; // Branching Ratio for Charm
225 gencharm->Init(); // Generating pT and Y parametrsation for the 4pi
226 ratiocharm = sigmacharm * brcharm * fNumberOfCollisions / sigmaReaction *
227 gencharm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
228 gencharm->SetPtRange(ptMin, ptMax);
229 gencharm->SetYRange(yMin, yMax);
230 gencharm->SetPhiRange(phiMin, phiMax);
231 gencharm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
84954c47 232
9538aedd 233 printf(">>> ratio charm %f\n",ratiocharm);
234 AddGenerator(gencharm,"Charm", ratiocharm);
235 fTotalRate+=ratiocharm;
84954c47 236
237// Generating Beauty Physics "Correlated Pairs"
9538aedd 238 AliGenParam * genbeauty = new AliGenParam(2, AliGenMUONlib::kBeauty, "Vogt", "Beauty");
239 genbeauty->SetPtRange(0,100.);
240 genbeauty->SetYRange(-8.,8);
241 genbeauty->SetPhiRange(0.,360.);
242 genbeauty->SetForceDecay(kSemiMuonic);
243 genbeauty->SetTrackingFlag(1);
244 Float_t ratiobeauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
245 Float_t sigmabeauty = 2. * 0.210e-3 * 0.84; //
246 Float_t brbeauty = 0.15; // Branching Ratio for Beauty
247 genbeauty->Init(); // Generating pT and Y parametrsation for the 4pi
248 ratiobeauty = sigmabeauty * brbeauty * fNumberOfCollisions / sigmaReaction *
249 genbeauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
250 genbeauty->SetPtRange(ptMin, ptMax);
251 genbeauty->SetYRange(yMin, yMax);
252 genbeauty->SetPhiRange(phiMin, phiMax);
253 genbeauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
84954c47 254
9538aedd 255 printf(">>> ratio beauty %f\n",ratiobeauty);
256 AddGenerator(genbeauty,"Beauty", ratiobeauty);
257 fTotalRate+=ratiobeauty;
84954c47 258
259// Generating Pion Physics
9538aedd 260 AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion");
261 genpion->SetPtRange(0,100.);
262 genpion->SetYRange(-8.,8);
263 genpion->SetPhiRange(0.,360.);
264 genpion->SetForceDecay(kPiToMu);
265 genpion->SetTrackingFlag(1);
266 Float_t ratiopion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
267 Float_t sigmapion = 0.93e-2; // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06
268 Float_t brpion = 0.9999; // Branching Ratio for Pion
269 genpion->Init(); // Generating pT and Y parametrsation for the 4pi
270 ratiopion = sigmapion * brpion * (0.80*fNumberOfParticipants+0.2*fNumberOfCollisions) / sigmaReaction * genpion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
271 genpion->SetPtRange(ptMin, ptMax);
272 genpion->SetYRange(yMin, yMax);
273 genpion->SetPhiRange(phiMin, phiMax);
274 genpion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
275 printf(">>> ratio pion %f\n",ratiopion);
276 AddGenerator(genpion,"Pion", ratiopion);
277 fTotalRate+=ratiopion;
84954c47 278
279// Generating Kaon Physics
9538aedd 280 AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon");
281 genkaon->SetPtRange(0,100.);
282 genkaon->SetYRange(-8.,8);
283 genkaon->SetPhiRange(0.,360.);
284 genkaon->SetForceDecay(kKaToMu);
285 genkaon->SetTrackingFlag(1);
286 Float_t ratiokaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
287 Float_t sigmakaon = 1.23e-4; // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06
288 Float_t brkaon = 0.6351 ; // Branching Ratio for Kaon
289 genkaon->Init(); // Generating pT and Y parametrsation for the 4pi
290 ratiokaon = sigmakaon * brkaon * (0.80*fNumberOfParticipants+0.2*fNumberOfCollisions)/ sigmaReaction * genkaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
291 genkaon->SetPtRange(ptMin, ptMax);
292 genkaon->SetYRange(yMin, yMax);
293 genkaon->SetPhiRange(phiMin, phiMax);
294 genkaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
295 printf(">>> ratio kaon %f\n",ratiokaon);
296 AddGenerator(genkaon,"Kaon", ratiokaon);
297 fTotalRate+=ratiokaon;
84954c47 298}
299
300//_________________________________________________________________________
301void AliGenMUONCocktail::Generate()
302{
303 //
304// Generate event
305 TIter next(fEntries);
306 AliGenCocktailEntry *entry = 0;
307 AliGenCocktailEntry *preventry = 0;
308 AliGenerator* gen = 0;
309
310 TObjArray *partArray = gAlice->GetMCApp()->Particles();
311
312//
313// Generate the vertex position used by all generators
314//
315 if(fVertexSmear == kPerEvent) Vertex();
9538aedd 316 Bool_t primordialTrigger = kFALSE;
84954c47 317
9538aedd 318 while(!primordialTrigger) {
84954c47 319
320 //Reseting stack
321 AliRunLoader * runloader = gAlice->GetRunLoader();
322 if (runloader)
323 if (runloader->Stack())
324 runloader->Stack()->Reset();
325 //
326 // Loop over generators and generate events
327 Int_t igen=0;
328 Int_t npart =0;
329
330 while((entry = (AliGenCocktailEntry*)next())) {
331 gen = entry->Generator();
332 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
333 if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
334 igen++;
335 if (igen ==1) entry->SetFirst(0);
336 else entry->SetFirst((partArray->GetEntriesFast())+1);
337 gen->SetNumberParticles(npart);
338 gen->Generate();
339 entry->SetLast(partArray->GetEntriesFast());
340 preventry = entry;
341 }
342 }
343 next.Reset();
344
345 // Tesitng primordial trigger : Muon pair in the MUON spectrometer acceptance 171,178 and pTCut
346 Int_t iPart;
347 fNGenerated++;
348 Int_t numberOfMuons=0;
349 // printf(">>>fNGenerated is %d\n",fNGenerated);
350 for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){
351 // gAlice->GetMCApp()->Particle(iPart)->Print();
352 if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13) &&
353 (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
354 (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
355 (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut) ) {
356 gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);
357 numberOfMuons++;
358 }
359 }
360 // printf(">>> Number of Muons is %d \n", numberOfMuons);
9538aedd 361 if (numberOfMuons >= fMuonMultiplicity ) primordialTrigger = kTRUE;
84954c47 362 }
363 //printf(">>> Trigger Accepted!!! \n");
364 fNSucceded++;
365 // Float_t Ratio = (Float_t) fNSucceded/fNGenerated;
366 // printf("Generated %d, Succeded %d and Ratio %f\n",fNGenerated, fNSucceded,Ratio);
367}
368
369
370
371
372
373