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