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 | |
38 | ClassImp(AliGenMUONCocktail) |
39 | |
40 | |
41 | //________________________________________________________________________ |
42 | AliGenMUONCocktail::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 | //_________________________________________________________________________ |
57 | AliGenMUONCocktail::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 | //_________________________________________________________________________ |
72 | AliGenMUONCocktail::~AliGenMUONCocktail() |
73 | { |
74 | // Destructor |
75 | delete fEntries; |
76 | } |
77 | |
78 | //_________________________________________________________________________ |
79 | void 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 | //_________________________________________________________________________ |
224 | void 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 | |