GetRandom returning impact parameter bin and flag for hard process added.
[u/mrichter/AliRoot.git] / FASTSIM / AliFastGlauber.cxx
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 Revision 1.3  2003/04/21 09:35:53  morsch
19 Protection against division by 0 in Binaries().
20
21 Revision 1.2  2003/04/14 14:23:44  morsch
22 Correction in Binaries().
23
24 Revision 1.1  2003/04/10 08:48:13  morsch
25 First commit.
26
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
39 ClassImp(AliFastGlauber)
40
41 TF1*    AliFastGlauber::fWSb      = NULL;     
42 TF2*    AliFastGlauber::fWSbz     = NULL;    
43 TF1*    AliFastGlauber::fWSz      = NULL;     
44 TF1*    AliFastGlauber::fWSta     = NULL;    
45 TF2*    AliFastGlauber::fWStarfi  = NULL; 
46 TF1*    AliFastGlauber::fWStaa    = NULL;   
47 TF1*    AliFastGlauber::fWSgeo    = NULL;   
48 TF1*    AliFastGlauber::fWSbinary = NULL;   
49 TF1*    AliFastGlauber::fWSN      = NULL;   
50 Float_t AliFastGlauber::fbMax     = 0.;
51
52 AliFastGlauber::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
62 void 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
127 void 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
137 void 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
147 void 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
157 void 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
167 void 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
177 void 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
187 void 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
197 Double_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
213 Double_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
230 Double_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
248 Double_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
261 Double_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
275 Double_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
308 Double_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
322 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
323 {
324 //
325 //  Number of binary collisions
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
335 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* par)
336 {
337 //
338 //  Number of hard processes per event
339 //
340     Double_t b     = x[0];
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     }
347     return y;
348 }
349
350 void 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
396 void 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);
404         mult = 6000./fWSN->Eval(0.) * mu;
405 }
406
407 void 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
435 Float_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
447 Double_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
456 Double_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
466 Double_t AliFastGlauber::Binaries(Double_t b)
467 {
468     //
469     // Return number of binary collisions normalized to 1 at b=0
470     //
471     
472     return fWSN->Eval(b)/fWSN->Eval(0.);
473 }