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