]> git.uio.no Git - u/mrichter/AliRoot.git/blame - FASTSIM/AliFastGlauber.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Transient pointer to AliRunDigitizer
[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$ */
041f7f97 17//
18// Utility class to make simple Glauber type calculations for collision geometries:
19// Impact parameter, production points, reaction plane dependence
20// The SimulateTrigger method can be used for simple MB and hard-process
21// (binary scaling) trigger studies.
22// Some basic quantities can be visualized directly.
23// The default set-up for PbPb collisions can be read from a file calling Init(1).
24//
25//
26// Author: andreas.morsch@cern.ch
5b3a5a5d 27
28// from AliRoot
29#include "AliFastGlauber.h"
30// from root
31#include <TH1F.h>
32#include <TF1.h>
33#include <TF2.h>
34#include <TCanvas.h>
35#include <TRandom.h>
36#include <TFile.h>
37
38ClassImp(AliFastGlauber)
39
041f7f97 40TF1* AliFastGlauber::fgWSb = NULL;
41TF2* AliFastGlauber::fgWSbz = NULL;
42TF1* AliFastGlauber::fgWSz = NULL;
43TF1* AliFastGlauber::fgWSta = NULL;
44TF2* AliFastGlauber::fgWStarfi = NULL;
45TF2* AliFastGlauber::fgWAlmond = NULL;
46TF1* AliFastGlauber::fgWStaa = NULL;
47TF1* AliFastGlauber::fgWSgeo = NULL;
48TF1* AliFastGlauber::fgWSbinary = NULL;
49TF1* AliFastGlauber::fgWSN = NULL;
50TF1* AliFastGlauber::fgWPathLength0 = NULL;
51TF1* AliFastGlauber::fgWPathLength = NULL;
52TF1* AliFastGlauber::fgWEnergyDensity = NULL;
53TF1* AliFastGlauber::fgWIntRadius = NULL;
54Float_t AliFastGlauber::fgBMax = 0.;
5b3a5a5d 55
56AliFastGlauber::AliFastGlauber()
57{
58// Default Constructor
59//
60// Defaults for Pb
61 SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
62 SetHardCrossSection();
63 SetMaxImpact();
64}
65
66void AliFastGlauber::Init(Int_t mode)
67{
68// Initialisation
69//
70// Wood-Saxon
71//
041f7f97 72 fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
73 fgWSb->SetParameter(0, fWSr0);
74 fgWSb->SetParameter(1, fWSd);
75 fgWSb->SetParameter(2, fWSw);
76 fgWSb->SetParameter(3, fWSn);
5b3a5a5d 77
041f7f97 78 fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
79 fgWSbz->SetParameter(0, fWSr0);
80 fgWSbz->SetParameter(1, fWSd);
81 fgWSbz->SetParameter(2, fWSw);
82 fgWSbz->SetParameter(3, fWSn);
5b3a5a5d 83
041f7f97 84 fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
85 fgWSz->SetParameter(0, fWSr0);
86 fgWSz->SetParameter(1, fWSd);
87 fgWSz->SetParameter(2, fWSw);
88 fgWSz->SetParameter(3, fWSn);
5b3a5a5d 89
90//
91// Thickness
92//
041f7f97 93 fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
5b3a5a5d 94
95//
96// Overlap Kernel
97//
041f7f97 98 fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
99 fgWStarfi->SetParameter(0, 0.);
100 fgWStarfi->SetNpx(200);
101 fgWStarfi->SetNpy(20);
5b3a5a5d 102//
f3a04204 103// Almond shaped interaction region
104//
041f7f97 105 fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
106 fgWAlmond->SetParameter(0, 0.);
107 fgWAlmond->SetNpx(200);
108 fgWAlmond->SetNpy(200);
f3a04204 109//
110// Path Length as a function of Phi
111//
041f7f97 112 fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
113 fgWPathLength0->SetParameter(0, 0.);
8de7e046 114// Pathlength definition
041f7f97 115 fgWPathLength0->SetParameter(1, 0.);
8de7e046 116
041f7f97 117 fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
f3a04204 118// Impact Parameter
041f7f97 119 fgWPathLength->SetParameter(0, 0.);
f3a04204 120// Number of interactions used for average
041f7f97 121 fgWPathLength->SetParameter(1, 1000.);
8de7e046 122// Pathlength definition
041f7f97 123 fgWPathLength->SetParameter(2, 0);
f3a04204 124
041f7f97 125 fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
126 fgWIntRadius->SetParameter(0, 0.);
2a103154 127
128
f3a04204 129//
5b3a5a5d 130// Overlap
131//
132 if (! mode) {
041f7f97 133 fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 0);
134 fgWStaa->SetNpx(100);
5b3a5a5d 135 } else {
136 TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
041f7f97 137 fgWStaa = (TF1*) f->Get("WStaa");
5b3a5a5d 138 }
139
2a103154 140//
041f7f97 141 fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
142 fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
2a103154 143
5b3a5a5d 144//
145// Geometrical Cross-Section
146//
041f7f97 147 fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 0);
148 fgWSgeo->SetNpx(100);
5b3a5a5d 149//
150// Hard cross section (~ binary collisions)
151//
041f7f97 152 fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
153 fgWSbinary->SetParameter(0, fSigmaHard); // mb
154 fgWSbinary->SetNpx(100);
5b3a5a5d 155//
156// Hard collisions per event
157//
041f7f97 158 fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
159 fgWSN->SetNpx(100);
5b3a5a5d 160}
161
162void AliFastGlauber::DrawWSb()
163{
164//
165// Draw Wood-Saxon Nuclear Density Function
166//
167 TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
168 c1->cd();
041f7f97 169 fgWSb->Draw();
5b3a5a5d 170}
171
172void AliFastGlauber::DrawOverlap()
173{
174//
175// Draw Overlap Function
176//
177 TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
178 c2->cd();
041f7f97 179 fgWStaa->Draw();
5b3a5a5d 180}
181
182void AliFastGlauber::DrawThickness()
183{
184//
185// Draw Thickness Function
186//
187 TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
188 c3->cd();
041f7f97 189 fgWSta->Draw();
5b3a5a5d 190}
191
192void AliFastGlauber::DrawGeo()
193{
194//
195// Draw Geometrical Cross-Section
196//
197 TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
198 c3->cd();
041f7f97 199 fgWSgeo->Draw();
5b3a5a5d 200}
201
202void AliFastGlauber::DrawBinary()
203{
204//
205// Draw Wood-Saxon Nuclear Density Function
206//
207 TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
208 c4->cd();
041f7f97 209 fgWSbinary->Draw();
5b3a5a5d 210}
211
212void AliFastGlauber::DrawN()
213{
214//
215// Draw Binaries per event
216//
217 TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
218 c5->cd();
041f7f97 219 fgWSN->Draw();
5b3a5a5d 220}
221
f3a04204 222void AliFastGlauber::DrawKernel(Double_t b)
5b3a5a5d 223{
224//
225// Draw Kernel
226//
227 TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
228 c6->cd();
041f7f97 229 fgWStarfi->SetParameter(0, b);
230 fgWStarfi->Draw();
5b3a5a5d 231}
232
f3a04204 233void AliFastGlauber::DrawAlmond(Double_t b)
234{
235//
2a103154 236// Draw Interaction Almond
f3a04204 237//
238 TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
239 c7->cd();
041f7f97 240 fgWAlmond->SetParameter(0, b);
241 fgWAlmond->Draw();
f3a04204 242}
243
8de7e046 244void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
f3a04204 245{
246//
2a103154 247// Draw Path Length
f3a04204 248//
249 TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
250 c8->cd();
041f7f97 251 fgWPathLength0->SetParameter(0, b);
252 fgWPathLength0->SetParameter(1, Double_t(iopt));
253 fgWPathLength0->SetMinimum(0.);
254 fgWPathLength0->SetMaximum(10.);
255 fgWPathLength0->Draw();
f3a04204 256}
257
8de7e046 258void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
f3a04204 259{
260//
2a103154 261// Draw Path Length
f3a04204 262//
263 TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
264 c9->cd();
041f7f97 265 fgWAlmond->SetParameter(0, b);
f3a04204 266
041f7f97 267 fgWPathLength->SetParameter(0, b);
268 fgWPathLength->SetParameter(1, Double_t (ni));
269 fgWPathLength->SetParameter(2, Double_t (iopt));
270 fgWPathLength->SetMinimum(0.);
271 fgWPathLength->SetMaximum(10.);
272 fgWPathLength->Draw();
f3a04204 273}
274
275void AliFastGlauber::DrawIntRadius(Double_t b)
276{
277//
2a103154 278// Draw Interaction Radius
f3a04204 279//
280 TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
281 c10->cd();
041f7f97 282 fgWIntRadius->SetParameter(0, b);
283 fgWIntRadius->SetMinimum(0.);
284 fgWIntRadius->Draw();
f3a04204 285}
286
2a103154 287void AliFastGlauber::DrawEnergyDensity()
288{
289//
290// Draw energy density
291//
292 TCanvas *c11 = new TCanvas("c11","Energy Density",400, 10, 600, 700);
293 c11->cd();
041f7f97 294 fgWEnergyDensity->SetMinimum(0.);
295 fgWEnergyDensity->Draw();
2a103154 296}
297
5b3a5a5d 298Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
299{
300//
301// Woods-Saxon Parameterisation
302// as a function of radius
303//
304 Double_t xx = x[0];
305 Double_t r0 = par[0];
306 Double_t d = par[1];
307 Double_t w = par[2];
308 Double_t n = par[3];
309
310 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
2a103154 311
5b3a5a5d 312 return y;
313}
314
315Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
316{
317//
318// Wood Saxon Parameterisation
319// as a function of z and b
320//
321 Double_t bb = x[0];
322 Double_t zz = x[1];
323 Double_t r0 = par[0];
324 Double_t d = par[1];
325 Double_t w = par[2];
326 Double_t n = par[3];
327 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
328 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
2a103154 329
5b3a5a5d 330 return y;
331}
332
333Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
334{
335//
336// Wood Saxon Parameterisation
337// as a function of z for fixed b
338//
339 Double_t bb = par[4];
340 Double_t zz = x[0];
341 Double_t r0 = par[0];
342 Double_t d = par[1];
343 Double_t w = par[2];
344 Double_t n = par[3];
345 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
346 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
2a103154 347
5b3a5a5d 348 return y;
349}
350
f86dad79 351Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
5b3a5a5d 352{
353//
354// Thickness function
355//
356 Double_t b = x[0];
041f7f97 357 fgWSz->SetParameter(4, b);
358 Double_t y = 2. * fgWSz->Integral(0., fgBMax);
5b3a5a5d 359 return y;
360}
361
362
363
364Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
365{
366//
367// Kernel for overlap function
368//
369 Double_t b = par[0];
370 Double_t r1 = x[0];
371 Double_t phi = x[1];
372 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
041f7f97 373 Double_t y = r1 * fgWSta->Eval(r1) * fgWSta->Eval(r2);
5b3a5a5d 374 return y;
375}
376
377
f3a04204 378Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
379{
380//
381// Almond shaped interaction region
382//
383 Double_t b = par[0];
384 Double_t xx = x[0] + b/2.;
385 Double_t yy = x[1];
386 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
387 Double_t phi = TMath::ATan2(yy,xx);
388
389 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
390//
391// Interaction probability calculated as product of thicknesses
392//
041f7f97 393 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
f3a04204 394 return y;
395}
396
397Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
398{
399//
400// Average radius at which interaction takes place
401//
402// Radius in the Almond
403 Double_t r = x[0];
404// Impact parameter
405 Double_t b = par[0];
041f7f97 406 fgWAlmond->SetParameter(0, b);
f3a04204 407// Steps in phi
408 Double_t dphi = 2. * TMath::Pi() / 100.;
409// Average over phi
410 Double_t phi = 0.;
411 Double_t y = 0.;
412
413 for (Int_t i = 0; i < 100; i++) {
414 Double_t xx = r * TMath::Cos(phi);
415 Double_t yy = r * TMath::Sin(phi);
041f7f97 416 y += fgWAlmond->Eval(xx,yy);
f3a04204 417 phi += dphi;
418 } // phi loop
419// Result multiplied by Jacobian (2 pi r)
420 return (2. * TMath::Pi() * y * r / 100.);
421}
422
423Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
424{
425//
426// Path Length as a function of phi for interaction point fixed at (0,0)
427//
428//
429// Steps in r
041f7f97 430 const Int_t kNp = 100;
431 const Double_t kDr = fgBMax/Double_t(kNp);
f3a04204 432// Impact parameter
433 Double_t b = par[0];
8de7e046 434// Path Length definition
435 Int_t iopt = Int_t(par[1]);
436
f3a04204 437// Phi direction in Almond
438 Double_t phi0 = x[0];
439 Double_t r = 0.;
440 Double_t rw = 0.;
441 Double_t w = 0.;
442// Step along radial direction phi
041f7f97 443 for (Int_t i = 0; i < kNp; i++) {
f3a04204 444//
445// Transform into target frame
446//
447 Double_t xx = r * TMath::Cos(phi0) + b / 2.;
448 Double_t yy = r * TMath::Sin(phi0);
449 Double_t phi = TMath::ATan2(yy, xx);
450
451 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
452// Radius in projectile frame
453 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
041f7f97 454 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
f3a04204 455
456 rw += y * r;
457 w += y;
041f7f97 458 r += kDr;
f3a04204 459 } // radial steps
460//
b6d0ed6b 461// My length definition (is exact for hard disk)
8de7e046 462 if (!iopt) {
463 return (2. * rw / w);
464 } else {
041f7f97 465 return TMath::Sqrt(2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01));
8de7e046 466 }
f3a04204 467}
468
469Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
470{
471//
472// Path Length as a function of phi
473// Interaction point from random distribution
474//
475//
476// r-steps
477//
041f7f97 478 const Int_t kNp = 100;
479 const Double_t kDr = fgBMax/Double_t(kNp);
f3a04204 480// Number of interactions
041f7f97 481 const Int_t kNpi = Int_t (par[1]);
f3a04204 482
483//
484// Impact parameter
485 Double_t b = par[0];
8de7e046 486// Path Length definition
487 Int_t iopt = Int_t(par[2]);
f3a04204 488// Phi direction
489 Double_t phi0 = x[0];
490
491 printf("phi0 %f \n", phi0);
492
493// Path length
494 Double_t l = 0.;
495
041f7f97 496 for (Int_t in = 0; in < kNpi; in ++) {
f3a04204 497 Double_t rw = 0.;
498 Double_t w = 0.;
499
500 // Interaction point
501 Double_t x0, y0;
041f7f97 502 fgWAlmond->GetRandom2(x0, y0);
f3a04204 503// Initial radius
504 Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
041f7f97 505 Int_t nps = Int_t ((fgBMax - r0)/kDr) - 1;
f3a04204 506
507 Double_t r = 0.;
508// Radial steps
509 for (Int_t i = 0; (i < nps ); i++) {
510
511// Transform into target frame
512 Double_t xx = x0 + r * TMath::Cos(phi0) + b / 2.;
513 Double_t yy = y0 + r * TMath::Sin(phi0);
514 Double_t phi = TMath::ATan2(yy, xx);
515 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
516// Radius in projectile frame
517 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
041f7f97 518 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
f3a04204 519
520 rw += y * r;
521 w += y;
041f7f97 522 r += kDr;
f3a04204 523 } // steps
524// Average over interactions
8de7e046 525 if (!iopt) {
526 l += (2. * rw / w);
527 } else {
041f7f97 528 l+= 2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01);
8de7e046 529 }
f3a04204 530 } // interactions
9db56006 531 if (!iopt)
041f7f97 532 return (l / Double_t(kNpi));
9db56006 533 else
041f7f97 534 return (TMath::Sqrt(l / Double_t(kNpi)));
f3a04204 535}
536
f86dad79 537Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
5b3a5a5d 538{
539//
540// Overlap function
541//
542 Double_t b = x[0];
041f7f97 543 fgWStarfi->SetParameter(0, b);
5b3a5a5d 544/*
545 Double_t al[2];
546 Double_t bl[2];
547 al[0] = 0.;
548 al[1] = 0.;
549 bl[0] = 6.6;
550 bl[1] = TMath::Pi();
551 Double_t err;
552
041f7f97 553 Double_t y = 2. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
5b3a5a5d 554 printf("WStaa: %f %f %f\n", b, y, err);
555*/
556//
557// MC Integration
558//
559 Double_t y = 0;
560 for (Int_t i = 0; i < 100000; i++)
561 {
562 Double_t phi = TMath::Pi() * gRandom->Rndm();
041f7f97 563 Double_t b1 = fgBMax * gRandom->Rndm();
564 y += fgWStarfi->Eval(b1, phi);
5b3a5a5d 565 }
041f7f97 566 y *= 2. * 0.1 * 208. * 208. * TMath::Pi() * fgBMax / 100000.;
5b3a5a5d 567 return y;
568}
569
f86dad79 570Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
5b3a5a5d 571{
572//
573// Geometrical Cross-Section
574//
575 Double_t b = x[0];
041f7f97 576 Double_t taa = fgWStaa->Eval(b);
577 const Double_t kSigma = 55.6; // mbarn
5b3a5a5d 578
041f7f97 579 Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- kSigma * taa)); // fm
5b3a5a5d 580 return y;
581}
582
583
584Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
585{
586//
c2140715 587// Number of binary collisions
5b3a5a5d 588//
589 Double_t b = x[0];
590 Double_t sigma = par[0];
041f7f97 591 Double_t taa = fgWStaa->Eval(b);
5b3a5a5d 592
593 Double_t y = 2. * TMath::Pi() * b * sigma * taa; // fm
594 return y;
595}
596
f86dad79 597Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
5b3a5a5d 598{
599//
c2140715 600// Number of hard processes per event
5b3a5a5d 601//
602 Double_t b = x[0];
041f7f97 603 Double_t y = fgWSbinary->Eval(b)/fgWSgeo->Eval(b);
5b3a5a5d 604 return y;
605}
606
2a103154 607Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
608{
609//
610// Initial energy density as a function of the impact parameter
611//
612 Double_t b = x[0];
613 Double_t rA = par[0];
b6d0ed6b 614//
615// Attention: area of transverse reaction zone in hard-sphere approximation !
616 Double_t saa = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA * rA
617 - b * TMath::Sqrt(rA * rA - b * b/ 4.);
041f7f97 618 Double_t taa = fgWStaa->Eval(b);
2a103154 619
620 return (taa/saa);
621}
622
5b3a5a5d 623void AliFastGlauber::SimulateTrigger(Int_t n)
624{
625 //
626 // Simulates Trigger
627 //
628 TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
629 TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
630 TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
631 TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
632
633 mbtH->SetXTitle("b [fm]");
634 hdtH->SetXTitle("b [fm]");
635 mbmH->SetXTitle("Multiplicity");
636 hdmH->SetXTitle("Multiplicity");
637
638 TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
639 c0->Divide(2,1);
640 TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
641 c1->Divide(1,2);
642
643 //
644 //
645 Init(1);
646 for (Int_t iev = 0; iev < n; iev++)
647 {
648 Float_t b, p, mult;
649 GetRandom(b, p, mult);
650 mbtH->Fill(b,1.);
651 hdtH->Fill(b, p);
652 mbmH->Fill(mult, 1.);
653 hdmH->Fill(mult, p);
654
655 c0->cd(1);
656 mbtH->Draw();
657 c0->cd(2);
658 hdtH->Draw();
659 c0->Update();
660
661 c1->cd(1);
662 mbmH->Draw();
663 c1->cd(2);
664 hdmH->Draw();
665 c1->Update();
666 }
667}
668
669void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
670{
671 //
672 // Gives back a random impact parameter, hard trigger probability and multiplicity
673 //
041f7f97 674 b = fgWSgeo->GetRandom();
675 Float_t mu = fgWSN->Eval(b);
5b3a5a5d 676 p = 1.-TMath::Exp(-mu);
041f7f97 677 mult = 6000./fgWSN->Eval(1.) * mu;
5b3a5a5d 678}
679
c2140715 680void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
681{
682 //
683 // Gives back a random impact parameter bin, and hard trigger decission
684 //
041f7f97 685 Float_t b = fgWSgeo->GetRandom();
686 Float_t mu = fgWSN->Eval(b) * fSigmaHard;
c2140715 687 Float_t p = 1.-TMath::Exp(-mu);
688 if (b < 5.) {
689 bin = 1;
690 } else if (b < 8.6) {
691 bin = 2;
692 } else if (b < 11.2) {
693 bin = 3;
694 } else if (b < 13.2) {
695 bin = 4;
204e011f 696 } else if (b < 15.0) {
c2140715 697 bin = 5;
204e011f 698 } else {
699 bin = 6;
c2140715 700 }
701
702 hard = kFALSE;
703
704 Float_t r = gRandom->Rndm();
705
706 if (r < p) hard = kTRUE;
707}
708
709
5b3a5a5d 710Float_t AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
711{
712 //
713 // Gives back a random impact parameter in the range bmin .. bmax
714 //
715
716 Float_t b = -1.;
717 while(b < bmin || b > bmax)
041f7f97 718 b = fgWSgeo->GetRandom();
5b3a5a5d 719 return b;
720}
721
722Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
723{
724 //
725 // Return cross-section integrated from b1 to b2
726 //
727
041f7f97 728 return fgWSgeo->Integral(b1, b2)/100.;
5b3a5a5d 729}
730
731Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
732{
733 //
734 // Return raction of hard cross-section integrated from b1 to b2
735 //
736
041f7f97 737 return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
5b3a5a5d 738}
739
740
741Double_t AliFastGlauber::Binaries(Double_t b)
742{
743 //
744 // Return number of binary collisions normalized to 1 at b=0
745 //
746
041f7f97 747 return fgWSN->Eval(b)/fgWSN->Eval(0.001);
5b3a5a5d 748}