5b3a5a5d |
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 | /* |
17 | $Log$ |
c2140715 |
18 | Revision 1.3 2003/04/21 09:35:53 morsch |
19 | Protection against division by 0 in Binaries(). |
20 | |
a7ef5179 |
21 | Revision 1.2 2003/04/14 14:23:44 morsch |
22 | Correction in Binaries(). |
23 | |
270181a5 |
24 | Revision 1.1 2003/04/10 08:48:13 morsch |
25 | First commit. |
26 | |
5b3a5a5d |
27 | */ |
28 | |
29 | // from AliRoot |
30 | #include "AliFastGlauber.h" |
31 | // from root |
32 | #include <TH1F.h> |
33 | #include <TF1.h> |
34 | #include <TF2.h> |
35 | #include <TCanvas.h> |
36 | #include <TRandom.h> |
37 | #include <TFile.h> |
38 | |
39 | ClassImp(AliFastGlauber) |
40 | |
41 | TF1* AliFastGlauber::fWSb = NULL; |
42 | TF2* AliFastGlauber::fWSbz = NULL; |
43 | TF1* AliFastGlauber::fWSz = NULL; |
44 | TF1* AliFastGlauber::fWSta = NULL; |
45 | TF2* AliFastGlauber::fWStarfi = NULL; |
46 | TF1* AliFastGlauber::fWStaa = NULL; |
47 | TF1* AliFastGlauber::fWSgeo = NULL; |
48 | TF1* AliFastGlauber::fWSbinary = NULL; |
49 | TF1* AliFastGlauber::fWSN = NULL; |
50 | Float_t AliFastGlauber::fbMax = 0.; |
51 | |
52 | AliFastGlauber::AliFastGlauber() |
53 | { |
54 | // Default Constructor |
55 | // |
56 | // Defaults for Pb |
57 | SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4); |
58 | SetHardCrossSection(); |
59 | SetMaxImpact(); |
60 | } |
61 | |
62 | void AliFastGlauber::Init(Int_t mode) |
63 | { |
64 | // Initialisation |
65 | // |
66 | // Wood-Saxon |
67 | // |
68 | fWSb = new TF1("WSb", WSb, 0, fbMax, 4); |
69 | fWSb->SetParameter(0, fWSr0); |
70 | fWSb->SetParameter(1, fWSd); |
71 | fWSb->SetParameter(2, fWSw); |
72 | fWSb->SetParameter(3, fWSn); |
73 | |
74 | fWSbz = new TF2("WSbz", WSbz, 0, fbMax, 4); |
75 | fWSbz->SetParameter(0, fWSr0); |
76 | fWSbz->SetParameter(1, fWSd); |
77 | fWSbz->SetParameter(2, fWSw); |
78 | fWSbz->SetParameter(3, fWSn); |
79 | |
80 | fWSz = new TF1("WSz", WSz, 0, fbMax, 5); |
81 | fWSz->SetParameter(0, fWSr0); |
82 | fWSz->SetParameter(1, fWSd); |
83 | fWSz->SetParameter(2, fWSw); |
84 | fWSz->SetParameter(3, fWSn); |
85 | |
86 | // |
87 | // Thickness |
88 | // |
89 | fWSta = new TF1("WSta", WSta, 0., fbMax, 0); |
90 | |
91 | // |
92 | // Overlap Kernel |
93 | // |
94 | fWStarfi = new TF2("WStarfi", WStarfi, 0., fbMax, 0., TMath::Pi(), 1); |
95 | fWStarfi->SetParameter(0, 0.); |
96 | fWStarfi->SetNpx(200); |
97 | fWStarfi->SetNpy(20); |
98 | // |
99 | // Overlap |
100 | // |
101 | if (! mode) { |
102 | fWStaa = new TF1("WStaa", WStaa, 0., fbMax, 0); |
103 | fWStaa->SetNpx(100); |
104 | } else { |
105 | TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root"); |
106 | fWStaa = (TF1*) f->Get("WStaa"); |
107 | } |
108 | |
109 | // |
110 | // Geometrical Cross-Section |
111 | // |
112 | fWSgeo = new TF1("WSgeo", WSgeo, 0., fbMax, 0); |
113 | fWSgeo->SetNpx(100); |
114 | // |
115 | // Hard cross section (~ binary collisions) |
116 | // |
117 | fWSbinary = new TF1("WSbinary", WSbinary, 0., fbMax, 1); |
118 | fWSbinary->SetParameter(0, fSigmaHard); // mb |
119 | fWSbinary->SetNpx(100); |
120 | // |
121 | // Hard collisions per event |
122 | // |
123 | fWSN = new TF1("WSN", WSN, 0., fbMax, 1); |
124 | fWSN->SetNpx(100); |
125 | } |
126 | |
127 | void AliFastGlauber::DrawWSb() |
128 | { |
129 | // |
130 | // Draw Wood-Saxon Nuclear Density Function |
131 | // |
132 | TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700); |
133 | c1->cd(); |
134 | fWSb->Draw(); |
135 | } |
136 | |
137 | void AliFastGlauber::DrawOverlap() |
138 | { |
139 | // |
140 | // Draw Overlap Function |
141 | // |
142 | TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700); |
143 | c2->cd(); |
144 | fWStaa->Draw(); |
145 | } |
146 | |
147 | void AliFastGlauber::DrawThickness() |
148 | { |
149 | // |
150 | // Draw Thickness Function |
151 | // |
152 | TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700); |
153 | c3->cd(); |
154 | fWSta->Draw(); |
155 | } |
156 | |
157 | void AliFastGlauber::DrawGeo() |
158 | { |
159 | // |
160 | // Draw Geometrical Cross-Section |
161 | // |
162 | TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700); |
163 | c3->cd(); |
164 | fWSgeo->Draw(); |
165 | } |
166 | |
167 | void AliFastGlauber::DrawBinary() |
168 | { |
169 | // |
170 | // Draw Wood-Saxon Nuclear Density Function |
171 | // |
172 | TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700); |
173 | c4->cd(); |
174 | fWSbinary->Draw(); |
175 | } |
176 | |
177 | void AliFastGlauber::DrawN() |
178 | { |
179 | // |
180 | // Draw Binaries per event |
181 | // |
182 | TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700); |
183 | c5->cd(); |
184 | fWSN->Draw(); |
185 | } |
186 | |
187 | void AliFastGlauber::DrawKernel() |
188 | { |
189 | // |
190 | // Draw Kernel |
191 | // |
192 | TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700); |
193 | c6->cd(); |
194 | fWStarfi->Draw(); |
195 | } |
196 | |
197 | Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par) |
198 | { |
199 | // |
200 | // Woods-Saxon Parameterisation |
201 | // as a function of radius |
202 | // |
203 | Double_t xx = x[0]; |
204 | Double_t r0 = par[0]; |
205 | Double_t d = par[1]; |
206 | Double_t w = par[2]; |
207 | Double_t n = par[3]; |
208 | |
209 | Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d)); |
210 | return y; |
211 | } |
212 | |
213 | Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par) |
214 | { |
215 | // |
216 | // Wood Saxon Parameterisation |
217 | // as a function of z and b |
218 | // |
219 | Double_t bb = x[0]; |
220 | Double_t zz = x[1]; |
221 | Double_t r0 = par[0]; |
222 | Double_t d = par[1]; |
223 | Double_t w = par[2]; |
224 | Double_t n = par[3]; |
225 | Double_t xx = TMath::Sqrt(bb*bb+zz*zz); |
226 | Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d)); |
227 | return y; |
228 | } |
229 | |
230 | Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par) |
231 | { |
232 | // |
233 | // Wood Saxon Parameterisation |
234 | // as a function of z for fixed b |
235 | // |
236 | Double_t bb = par[4]; |
237 | Double_t zz = x[0]; |
238 | Double_t r0 = par[0]; |
239 | Double_t d = par[1]; |
240 | Double_t w = par[2]; |
241 | Double_t n = par[3]; |
242 | Double_t xx = TMath::Sqrt(bb*bb+zz*zz); |
243 | Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d)); |
244 | // if (y < 1.e-6) y = 0; |
245 | return y; |
246 | } |
247 | |
248 | Double_t AliFastGlauber::WSta(Double_t* x, Double_t* par) |
249 | { |
250 | // |
251 | // Thickness function |
252 | // |
253 | Double_t b = x[0]; |
254 | fWSz->SetParameter(4, b); |
255 | Double_t y = 2. * fWSz->Integral(0., fbMax); |
256 | return y; |
257 | } |
258 | |
259 | |
260 | |
261 | Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par) |
262 | { |
263 | // |
264 | // Kernel for overlap function |
265 | // |
266 | Double_t b = par[0]; |
267 | Double_t r1 = x[0]; |
268 | Double_t phi = x[1]; |
269 | Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); |
270 | Double_t y = r1 * fWSta->Eval(r1) * fWSta->Eval(r2); |
271 | return y; |
272 | } |
273 | |
274 | |
275 | Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par) |
276 | { |
277 | // |
278 | // Overlap function |
279 | // |
280 | Double_t b = x[0]; |
281 | fWStarfi->SetParameter(0, b); |
282 | /* |
283 | Double_t al[2]; |
284 | Double_t bl[2]; |
285 | al[0] = 0.; |
286 | al[1] = 0.; |
287 | bl[0] = 6.6; |
288 | bl[1] = TMath::Pi(); |
289 | Double_t err; |
290 | |
291 | Double_t y = 2. * fWStarfi->IntegralMultiple(2, al, bl, 0.001, err); |
292 | printf("WStaa: %f %f %f\n", b, y, err); |
293 | */ |
294 | // |
295 | // MC Integration |
296 | // |
297 | Double_t y = 0; |
298 | for (Int_t i = 0; i < 100000; i++) |
299 | { |
300 | Double_t phi = TMath::Pi() * gRandom->Rndm(); |
301 | Double_t b1 = fbMax * gRandom->Rndm(); |
302 | y += fWStarfi->Eval(b1, phi); |
303 | } |
304 | y *= 2. * 0.1 * 208. * 208. * TMath::Pi() * fbMax / 100000.; |
305 | return y; |
306 | } |
307 | |
308 | Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par) |
309 | { |
310 | // |
311 | // Geometrical Cross-Section |
312 | // |
313 | Double_t b = x[0]; |
314 | Double_t taa = fWStaa->Eval(b); |
315 | const Double_t sigma = 55.6; // mbarn |
316 | |
317 | Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigma * taa)); // fm |
318 | return y; |
319 | } |
320 | |
321 | |
322 | Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par) |
323 | { |
324 | // |
c2140715 |
325 | // Number of binary collisions |
5b3a5a5d |
326 | // |
327 | Double_t b = x[0]; |
328 | Double_t sigma = par[0]; |
329 | Double_t taa = fWStaa->Eval(b); |
330 | |
331 | Double_t y = 2. * TMath::Pi() * b * sigma * taa; // fm |
332 | return y; |
333 | } |
334 | |
335 | Double_t AliFastGlauber::WSN(Double_t* x, Double_t* par) |
336 | { |
337 | // |
c2140715 |
338 | // Number of hard processes per event |
5b3a5a5d |
339 | // |
340 | Double_t b = x[0]; |
a7ef5179 |
341 | Double_t y; |
342 | if (b != 0.) { |
343 | y = fWSbinary->Eval(b)/fWSgeo->Eval(b); |
344 | } else { |
345 | y = fWSbinary->Eval(1.e-4)/fWSgeo->Eval(1.e-4); |
346 | } |
5b3a5a5d |
347 | return y; |
348 | } |
349 | |
350 | void AliFastGlauber::SimulateTrigger(Int_t n) |
351 | { |
352 | // |
353 | // Simulates Trigger |
354 | // |
355 | TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.); |
356 | TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.); |
357 | TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.); |
358 | TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.); |
359 | |
360 | mbtH->SetXTitle("b [fm]"); |
361 | hdtH->SetXTitle("b [fm]"); |
362 | mbmH->SetXTitle("Multiplicity"); |
363 | hdmH->SetXTitle("Multiplicity"); |
364 | |
365 | TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700); |
366 | c0->Divide(2,1); |
367 | TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700); |
368 | c1->Divide(1,2); |
369 | |
370 | // |
371 | // |
372 | Init(1); |
373 | for (Int_t iev = 0; iev < n; iev++) |
374 | { |
375 | Float_t b, p, mult; |
376 | GetRandom(b, p, mult); |
377 | mbtH->Fill(b,1.); |
378 | hdtH->Fill(b, p); |
379 | mbmH->Fill(mult, 1.); |
380 | hdmH->Fill(mult, p); |
381 | |
382 | c0->cd(1); |
383 | mbtH->Draw(); |
384 | c0->cd(2); |
385 | hdtH->Draw(); |
386 | c0->Update(); |
387 | |
388 | c1->cd(1); |
389 | mbmH->Draw(); |
390 | c1->cd(2); |
391 | hdmH->Draw(); |
392 | c1->Update(); |
393 | } |
394 | } |
395 | |
396 | void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult) |
397 | { |
398 | // |
399 | // Gives back a random impact parameter, hard trigger probability and multiplicity |
400 | // |
401 | b = fWSgeo->GetRandom(); |
402 | Float_t mu = fWSN->Eval(b); |
403 | p = 1.-TMath::Exp(-mu); |
a7ef5179 |
404 | mult = 6000./fWSN->Eval(0.) * mu; |
5b3a5a5d |
405 | } |
406 | |
c2140715 |
407 | void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard) |
408 | { |
409 | // |
410 | // Gives back a random impact parameter bin, and hard trigger decission |
411 | // |
412 | Float_t b = fWSgeo->GetRandom(); |
413 | Float_t mu = fWSN->Eval(b) * fSigmaHard; |
414 | Float_t p = 1.-TMath::Exp(-mu); |
415 | if (b < 5.) { |
416 | bin = 1; |
417 | } else if (b < 8.6) { |
418 | bin = 2; |
419 | } else if (b < 11.2) { |
420 | bin = 3; |
421 | } else if (b < 13.2) { |
422 | bin = 4; |
423 | } else { |
424 | bin = 5; |
425 | } |
426 | |
427 | hard = kFALSE; |
428 | |
429 | Float_t r = gRandom->Rndm(); |
430 | |
431 | if (r < p) hard = kTRUE; |
432 | } |
433 | |
434 | |
5b3a5a5d |
435 | Float_t AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax) |
436 | { |
437 | // |
438 | // Gives back a random impact parameter in the range bmin .. bmax |
439 | // |
440 | |
441 | Float_t b = -1.; |
442 | while(b < bmin || b > bmax) |
443 | b = fWSgeo->GetRandom(); |
444 | return b; |
445 | } |
446 | |
447 | Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2) |
448 | { |
449 | // |
450 | // Return cross-section integrated from b1 to b2 |
451 | // |
452 | |
453 | return fWSgeo->Integral(b1, b2)/100.; |
454 | } |
455 | |
456 | Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) |
457 | { |
458 | // |
459 | // Return raction of hard cross-section integrated from b1 to b2 |
460 | // |
461 | |
462 | return fWSbinary->Integral(b1, b2)/fWSbinary->Integral(0., 100.); |
463 | } |
464 | |
465 | |
466 | Double_t AliFastGlauber::Binaries(Double_t b) |
467 | { |
468 | // |
469 | // Return number of binary collisions normalized to 1 at b=0 |
470 | // |
471 | |
270181a5 |
472 | return fWSN->Eval(b)/fWSN->Eval(0.); |
5b3a5a5d |
473 | } |