GetRandom returning impact parameter bin and flag for hard process added.
[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$
c2140715 18Revision 1.3 2003/04/21 09:35:53 morsch
19Protection against division by 0 in Binaries().
20
a7ef5179 21Revision 1.2 2003/04/14 14:23:44 morsch
22Correction in Binaries().
23
270181a5 24Revision 1.1 2003/04/10 08:48:13 morsch
25First 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
39ClassImp(AliFastGlauber)
40
41TF1* AliFastGlauber::fWSb = NULL;
42TF2* AliFastGlauber::fWSbz = NULL;
43TF1* AliFastGlauber::fWSz = NULL;
44TF1* AliFastGlauber::fWSta = NULL;
45TF2* AliFastGlauber::fWStarfi = NULL;
46TF1* AliFastGlauber::fWStaa = NULL;
47TF1* AliFastGlauber::fWSgeo = NULL;
48TF1* AliFastGlauber::fWSbinary = NULL;
49TF1* AliFastGlauber::fWSN = NULL;
50Float_t AliFastGlauber::fbMax = 0.;
51
52AliFastGlauber::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
62void 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
127void 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
137void 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
147void 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
157void 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
167void 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
177void 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
187void 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
197Double_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
213Double_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
230Double_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
248Double_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
261Double_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
275Double_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
308Double_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
322Double_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
335Double_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
350void 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
396void 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 407void 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 435Float_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
447Double_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
456Double_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
466Double_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}