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