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