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