Coding violations (Eric)
[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
18// Classe to create the MUON coktail for physics in the Alice muon spectrometer
19// Gines Martinez, jan 2004, Nantes martinez@in2p3.fr
20
21
22//
23
24#include <TList.h>
25#include <TObjArray.h>
26#include <TF1.h>
27#include <TParticle.h>
28
29#include "AliGenParam.h"
30#include "AliGenMUONlib.h"
31#include "AliGenMUONCocktail.h"
32#include "AliGenCocktailEntry.h"
33#include "AliCollisionGeometry.h"
34#include "AliRun.h"
35#include "AliMC.h"
36#include "AliStack.h"
37
38ClassImp(AliGenMUONCocktail)
39
40
41//________________________________________________________________________
42AliGenMUONCocktail::AliGenMUONCocktail()
43 :AliGenCocktail()
44{
45// Constructor
46 fTotalRate =0;
47 fNSucceded=0;
48 fNGenerated=0;
49 fMuonMultiplicity=1;
50 fMuonPtCut= 1.;
51 fMuonThetaMinCut=171.;
52 fMuonThetaMaxCut=178.;
53 fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
54 //
55}
56//_________________________________________________________________________
57AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
58 AliGenCocktail(cocktail)
59{
60// Copy constructor
61 fTotalRate =0;
62 fNSucceded=0;
63 fNGenerated=0;
64 fMuonMultiplicity=1;
65 fMuonPtCut= 1.;
66 fMuonThetaMinCut=171.;
67 fMuonThetaMaxCut=178.;
68 fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
69 //
70}
71//_________________________________________________________________________
72AliGenMUONCocktail::~AliGenMUONCocktail()
73{
74// Destructor
75 delete fEntries;
76}
77
78//_________________________________________________________________________
79void AliGenMUONCocktail::Init()
80{
81 // Defining MUON physics cocktail
82
83 // Kinematical limits for particle generation
84 Float_t ptMin = fPtMin;
85 Float_t ptMax = fPtMax;
86 Float_t yMin = fYMin;;
87 Float_t yMax = fYMax;;
88 Float_t phiMin = fPhiMin*180./TMath::Pi();
89 Float_t phiMax = fPhiMax*180./TMath::Pi();
90 printf(">>> Kinematical range pT:%f:%f y:%f:%f Phi:%f:%f\n",ptMin,ptMax,yMin,yMax,phiMin,phiMax);
91
92 Float_t sigma_reaction = 0.072; // MINB pp at LHC energies 72 mb
93
94 // Generating J/Psi Physics
95 AliGenParam * gen_jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi");
96 // 4pi Generation
97 gen_jpsi->SetPtRange(0,100.);
98 gen_jpsi->SetYRange(-8.,8);
99 gen_jpsi->SetPhiRange(0.,360.);
100 gen_jpsi->SetForceDecay(kDiMuon);
101 gen_jpsi->SetTrackingFlag(1);
102 // Calculation of the paritcle multiplicity per event in the muonic channel
103 Float_t ratio_jpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
104 Float_t sigma_jpsi = 31.0e-6 * 0.437; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing
105 Float_t br_jpsi = 0.0588; // Branching Ratio for JPsi
106 gen_jpsi->Init(); // Generating pT and Y parametrsation for the 4pi
107 ratio_jpsi = sigma_jpsi * br_jpsi * fNumberOfCollisions / sigma_reaction * gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
108 printf(">>> ratio jpsi %g et %g Ncol %g sigma %g\n",ratio_jpsi,gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigma_jpsi );
109 // Generation in the kinematical limits of AliGenMUONCocktail
110 gen_jpsi->SetPtRange(ptMin, ptMax);
111 gen_jpsi->SetYRange(yMin, yMax);
112 gen_jpsi->SetPhiRange(phiMin, phiMax);
113 gen_jpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
114 // Adding Generator
115 AddGenerator(gen_jpsi, "Jpsi", ratio_jpsi);
116 fTotalRate+=ratio_jpsi;
117
118// Generating Upsilon Physics
119 AliGenParam * gen_upsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon");
120 gen_upsilon->SetPtRange(0,100.);
121 gen_upsilon->SetYRange(-8.,8);
122 gen_upsilon->SetPhiRange(0.,360.);
123 gen_upsilon->SetForceDecay(kDiMuon);
124 gen_upsilon->SetTrackingFlag(1);
125 Float_t ratio_upsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
126 Float_t sigma_upsilon = 0.501e-6 * 0.674; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing
127 Float_t br_upsilon = 0.0248; // Branching Ratio for Upsilon
128 gen_upsilon->Init(); // Generating pT and Y parametrsation for the 4pi
129 ratio_upsilon = sigma_upsilon * br_upsilon * fNumberOfCollisions / sigma_reaction * gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
130 printf(">>> ratio upsilon %g et %g\n",ratio_upsilon, gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
131 gen_upsilon->SetPtRange(ptMin, ptMax);
132 gen_upsilon->SetYRange(yMin, yMax);
133 gen_upsilon->SetPhiRange(phiMin, phiMax);
134 gen_upsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
135 AddGenerator(gen_upsilon,"Upsilon", ratio_upsilon);
136 fTotalRate+=ratio_upsilon;
137
138// Generating Charm Physics
139 AliGenParam * gen_charm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm");
140 gen_charm->SetPtRange(0,100.);
141 gen_charm->SetYRange(-8.,8);
142 gen_charm->SetPhiRange(0.,360.);
143 gen_charm->SetForceDecay(kSemiMuonic);
144 gen_charm->SetTrackingFlag(1);
145 Float_t ratio_charm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
146 Float_t sigma_charm = 2. * 6.64e-3 * 0.65 ; //
147 Float_t br_charm = 0.12; // Branching Ratio for Charm
148 gen_charm->Init(); // Generating pT and Y parametrsation for the 4pi
149 ratio_charm = sigma_charm * br_charm * fNumberOfCollisions / sigma_reaction *
150 gen_charm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
151 gen_charm->SetPtRange(ptMin, ptMax);
152 gen_charm->SetYRange(yMin, yMax);
153 gen_charm->SetPhiRange(phiMin, phiMax);
154 gen_charm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
155
156 printf(">>> ratio charm %f\n",ratio_charm);
157 AddGenerator(gen_charm,"Charm", ratio_charm);
158 fTotalRate+=ratio_charm;
159
160// Generating Beauty Physics "Correlated Pairs"
161 AliGenParam * gen_beauty = new AliGenParam(2, AliGenMUONlib::kBeauty, "Vogt", "Beauty");
162 gen_beauty->SetPtRange(0,100.);
163 gen_beauty->SetYRange(-8.,8);
164 gen_beauty->SetPhiRange(0.,360.);
165 gen_beauty->SetForceDecay(kSemiMuonic);
166 gen_beauty->SetTrackingFlag(1);
167 Float_t ratio_beauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
168 Float_t sigma_beauty = 2. * 0.210e-3 * 0.84; //
169 Float_t br_beauty = 0.15; // Branching Ratio for Beauty
170 gen_beauty->Init(); // Generating pT and Y parametrsation for the 4pi
171 ratio_beauty = sigma_beauty * br_beauty * fNumberOfCollisions / sigma_reaction *
172 gen_beauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
173 gen_beauty->SetPtRange(ptMin, ptMax);
174 gen_beauty->SetYRange(yMin, yMax);
175 gen_beauty->SetPhiRange(phiMin, phiMax);
176 gen_beauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
177
178 printf(">>> ratio beauty %f\n",ratio_beauty);
179 AddGenerator(gen_beauty,"Beauty", ratio_beauty);
180 fTotalRate+=ratio_beauty;
181
182// Generating Pion Physics
183 AliGenParam * gen_pion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion");
184 gen_pion->SetPtRange(0,100.);
185 gen_pion->SetYRange(-8.,8);
186 gen_pion->SetPhiRange(0.,360.);
187 gen_pion->SetForceDecay(kPiToMu);
188 gen_pion->SetTrackingFlag(1);
189 Float_t ratio_pion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
190 Float_t sigma_pion = 1.6e-2; // A ojo TO be studied in detail.
191 Float_t br_pion = 0.9999; // Branching Ratio for Pion
192 gen_pion->Init(); // Generating pT and Y parametrsation for the 4pi
193 ratio_pion = sigma_pion * br_pion * fNumberOfParticipants / sigma_reaction * gen_pion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
194 gen_pion->SetPtRange(ptMin, ptMax);
195 gen_pion->SetYRange(yMin, yMax);
196 gen_pion->SetPhiRange(phiMin, phiMax);
197 gen_pion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
198 printf(">>> ratio pion %f\n",ratio_pion);
199 AddGenerator(gen_pion,"Pion", ratio_pion);
200 fTotalRate+=ratio_pion;
201
202// Generating Kaon Physics
203 AliGenParam * gen_kaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon");
204 gen_kaon->SetPtRange(0,100.);
205 gen_kaon->SetYRange(-8.,8);
206 gen_kaon->SetPhiRange(0.,360.);
207 gen_kaon->SetForceDecay(kKaToMu);
208 gen_kaon->SetTrackingFlag(1);
209 Float_t ratio_kaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
210 Float_t sigma_kaon = 1.8e-4; // OJO
211 Float_t br_kaon = 0.6351 ; // Branching Ratio for Kaon
212 gen_kaon->Init(); // Generating pT and Y parametrsation for the 4pi
213 ratio_kaon = sigma_kaon * br_kaon * fNumberOfParticipants/ sigma_reaction * gen_kaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
214 gen_kaon->SetPtRange(ptMin, ptMax);
215 gen_kaon->SetYRange(yMin, yMax);
216 gen_kaon->SetPhiRange(phiMin, phiMax);
217 gen_kaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
218 printf(">>> ratio kaon %f\n",ratio_kaon);
219 AddGenerator(gen_kaon,"Kaon", ratio_kaon);
220 fTotalRate+=ratio_kaon;
221}
222
223//_________________________________________________________________________
224void AliGenMUONCocktail::Generate()
225{
226 //
227// Generate event
228 TIter next(fEntries);
229 AliGenCocktailEntry *entry = 0;
230 AliGenCocktailEntry *preventry = 0;
231 AliGenerator* gen = 0;
232
233 TObjArray *partArray = gAlice->GetMCApp()->Particles();
234
235//
236// Generate the vertex position used by all generators
237//
238 if(fVertexSmear == kPerEvent) Vertex();
239 Bool_t PrimordialTrigger = kFALSE;
240
241 while(!PrimordialTrigger) {
242
243 //Reseting stack
244 AliRunLoader * runloader = gAlice->GetRunLoader();
245 if (runloader)
246 if (runloader->Stack())
247 runloader->Stack()->Reset();
248 //
249 // Loop over generators and generate events
250 Int_t igen=0;
251 Int_t npart =0;
252
253 while((entry = (AliGenCocktailEntry*)next())) {
254 gen = entry->Generator();
255 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
256 if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
257 igen++;
258 if (igen ==1) entry->SetFirst(0);
259 else entry->SetFirst((partArray->GetEntriesFast())+1);
260 gen->SetNumberParticles(npart);
261 gen->Generate();
262 entry->SetLast(partArray->GetEntriesFast());
263 preventry = entry;
264 }
265 }
266 next.Reset();
267
268 // Tesitng primordial trigger : Muon pair in the MUON spectrometer acceptance 171,178 and pTCut
269 Int_t iPart;
270 fNGenerated++;
271 Int_t numberOfMuons=0;
272 // printf(">>>fNGenerated is %d\n",fNGenerated);
273 for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){
274 // gAlice->GetMCApp()->Particle(iPart)->Print();
275 if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13) &&
276 (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
277 (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
278 (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut) ) {
279 gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);
280 numberOfMuons++;
281 }
282 }
283 // printf(">>> Number of Muons is %d \n", numberOfMuons);
284 if (numberOfMuons >= fMuonMultiplicity ) PrimordialTrigger = kTRUE;
285 }
286 //printf(">>> Trigger Accepted!!! \n");
287 fNSucceded++;
288 // Float_t Ratio = (Float_t) fNSucceded/fNGenerated;
289 // printf("Generated %d, Succeded %d and Ratio %f\n",fNGenerated, fNSucceded,Ratio);
290}
291
292
293
294
295
296