]> 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
Adding MUON output to the TrackRefs files. To be used for the evaluation/comparisons...
[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
a2f2f511 31#include <TStyle.h>
5b3a5a5d 32#include <TH1F.h>
a2f2f511 33#include <TH2F.h>
5b3a5a5d 34#include <TF1.h>
35#include <TF2.h>
36#include <TCanvas.h>
37#include <TRandom.h>
38#include <TFile.h>
39
40ClassImp(AliFastGlauber)
41
041f7f97 42TF1* AliFastGlauber::fgWSb = NULL;
43TF2* AliFastGlauber::fgWSbz = NULL;
44TF1* AliFastGlauber::fgWSz = NULL;
45TF1* AliFastGlauber::fgWSta = NULL;
46TF2* AliFastGlauber::fgWStarfi = NULL;
47TF2* AliFastGlauber::fgWAlmond = NULL;
a2f2f511 48//TF2* AliFastGlauber::fWAlmondCurrent = NULL;
041f7f97 49TF1* AliFastGlauber::fgWStaa = NULL;
50TF1* AliFastGlauber::fgWSgeo = NULL;
51TF1* AliFastGlauber::fgWSbinary = NULL;
52TF1* AliFastGlauber::fgWSN = NULL;
53TF1* AliFastGlauber::fgWPathLength0 = NULL;
54TF1* AliFastGlauber::fgWPathLength = NULL;
55TF1* AliFastGlauber::fgWEnergyDensity = NULL;
56TF1* AliFastGlauber::fgWIntRadius = NULL;
57Float_t AliFastGlauber::fgBMax = 0.;
5b3a5a5d 58
59AliFastGlauber::AliFastGlauber()
60{
61// Default Constructor
62//
63// Defaults for Pb
64 SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
65 SetHardCrossSection();
66 SetMaxImpact();
a2f2f511 67 SetLengthDefinition();
5b3a5a5d 68}
69
70void AliFastGlauber::Init(Int_t mode)
71{
72// Initialisation
a2f2f511 73// mode = 0; all functions are calculated
74// mode = 1; overlap function is read from file (for Pb-Pb only)
75// mode = 2; interaction almond functions are read from file
76// (for Pb-Pb only); USE THIS FOR PATH LENGTH CALC.!
77//
78
5b3a5a5d 79//
80// Wood-Saxon
81//
041f7f97 82 fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
83 fgWSb->SetParameter(0, fWSr0);
84 fgWSb->SetParameter(1, fWSd);
85 fgWSb->SetParameter(2, fWSw);
86 fgWSb->SetParameter(3, fWSn);
5b3a5a5d 87
041f7f97 88 fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
89 fgWSbz->SetParameter(0, fWSr0);
90 fgWSbz->SetParameter(1, fWSd);
91 fgWSbz->SetParameter(2, fWSw);
92 fgWSbz->SetParameter(3, fWSn);
5b3a5a5d 93
041f7f97 94 fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
95 fgWSz->SetParameter(0, fWSr0);
96 fgWSz->SetParameter(1, fWSd);
97 fgWSz->SetParameter(2, fWSw);
98 fgWSz->SetParameter(3, fWSn);
5b3a5a5d 99
100//
101// Thickness
102//
041f7f97 103 fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
5b3a5a5d 104
105//
106// Overlap Kernel
107//
041f7f97 108 fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
109 fgWStarfi->SetParameter(0, 0.);
110 fgWStarfi->SetNpx(200);
111 fgWStarfi->SetNpy(20);
5b3a5a5d 112//
f3a04204 113// Almond shaped interaction region
114//
041f7f97 115 fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
116 fgWAlmond->SetParameter(0, 0.);
117 fgWAlmond->SetNpx(200);
a2f2f511 118 fgWAlmond->SetNpy(200);
119
120 if(mode==2) {
121 printf(" Reading interaction almonds from file\n");
122 Char_t almondName[100];
123 TFile* ff = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
124 for(Int_t k=0; k<40; k++) {
125 sprintf(almondName,"WAlmondFixedB%d",k);
126 fWAlmondCurrent = (TF2*)ff->Get(almondName);
127 new(&fWAlmondFixedB[k]) TF2(*fWAlmondCurrent);
128 }
129 }
f3a04204 130//
131// Path Length as a function of Phi
132//
041f7f97 133 fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
134 fgWPathLength0->SetParameter(0, 0.);
8de7e046 135// Pathlength definition
041f7f97 136 fgWPathLength0->SetParameter(1, 0.);
8de7e046 137
041f7f97 138 fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
f3a04204 139// Impact Parameter
041f7f97 140 fgWPathLength->SetParameter(0, 0.);
f3a04204 141// Number of interactions used for average
041f7f97 142 fgWPathLength->SetParameter(1, 1000.);
8de7e046 143// Pathlength definition
041f7f97 144 fgWPathLength->SetParameter(2, 0);
f3a04204 145
041f7f97 146 fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
147 fgWIntRadius->SetParameter(0, 0.);
2a103154 148
149
f3a04204 150//
5b3a5a5d 151// Overlap
152//
153 if (! mode) {
041f7f97 154 fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 0);
155 fgWStaa->SetNpx(100);
5b3a5a5d 156 } else {
157 TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
041f7f97 158 fgWStaa = (TF1*) f->Get("WStaa");
5b3a5a5d 159 }
160
2a103154 161//
041f7f97 162 fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
163 fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
2a103154 164
5b3a5a5d 165//
166// Geometrical Cross-Section
167//
041f7f97 168 fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 0);
169 fgWSgeo->SetNpx(100);
5b3a5a5d 170//
171// Hard cross section (~ binary collisions)
172//
041f7f97 173 fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
174 fgWSbinary->SetParameter(0, fSigmaHard); // mb
175 fgWSbinary->SetNpx(100);
5b3a5a5d 176//
177// Hard collisions per event
178//
041f7f97 179 fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
180 fgWSN->SetNpx(100);
5b3a5a5d 181}
182
183void AliFastGlauber::DrawWSb()
184{
185//
186// Draw Wood-Saxon Nuclear Density Function
187//
188 TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
189 c1->cd();
041f7f97 190 fgWSb->Draw();
5b3a5a5d 191}
192
193void AliFastGlauber::DrawOverlap()
194{
195//
196// Draw Overlap Function
197//
198 TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
199 c2->cd();
041f7f97 200 fgWStaa->Draw();
5b3a5a5d 201}
202
203void AliFastGlauber::DrawThickness()
204{
205//
206// Draw Thickness Function
207//
208 TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
209 c3->cd();
041f7f97 210 fgWSta->Draw();
5b3a5a5d 211}
212
213void AliFastGlauber::DrawGeo()
214{
215//
216// Draw Geometrical Cross-Section
217//
218 TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
219 c3->cd();
041f7f97 220 fgWSgeo->Draw();
5b3a5a5d 221}
222
223void AliFastGlauber::DrawBinary()
224{
225//
a2f2f511 226// Draw Binary Cross-Section
5b3a5a5d 227//
228 TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
229 c4->cd();
041f7f97 230 fgWSbinary->Draw();
5b3a5a5d 231}
232
233void AliFastGlauber::DrawN()
234{
235//
236// Draw Binaries per event
237//
238 TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
239 c5->cd();
041f7f97 240 fgWSN->Draw();
5b3a5a5d 241}
242
f3a04204 243void AliFastGlauber::DrawKernel(Double_t b)
5b3a5a5d 244{
245//
246// Draw Kernel
247//
248 TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
249 c6->cd();
041f7f97 250 fgWStarfi->SetParameter(0, b);
251 fgWStarfi->Draw();
5b3a5a5d 252}
253
f3a04204 254void AliFastGlauber::DrawAlmond(Double_t b)
255{
256//
2a103154 257// Draw Interaction Almond
f3a04204 258//
259 TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
260 c7->cd();
041f7f97 261 fgWAlmond->SetParameter(0, b);
262 fgWAlmond->Draw();
f3a04204 263}
264
8de7e046 265void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
f3a04204 266{
267//
2a103154 268// Draw Path Length
f3a04204 269//
270 TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
271 c8->cd();
041f7f97 272 fgWPathLength0->SetParameter(0, b);
273 fgWPathLength0->SetParameter(1, Double_t(iopt));
274 fgWPathLength0->SetMinimum(0.);
275 fgWPathLength0->SetMaximum(10.);
276 fgWPathLength0->Draw();
f3a04204 277}
278
8de7e046 279void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
f3a04204 280{
281//
2a103154 282// Draw Path Length
f3a04204 283//
284 TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
285 c9->cd();
041f7f97 286 fgWAlmond->SetParameter(0, b);
f3a04204 287
041f7f97 288 fgWPathLength->SetParameter(0, b);
289 fgWPathLength->SetParameter(1, Double_t (ni));
290 fgWPathLength->SetParameter(2, Double_t (iopt));
291 fgWPathLength->SetMinimum(0.);
292 fgWPathLength->SetMaximum(10.);
293 fgWPathLength->Draw();
f3a04204 294}
295
296void AliFastGlauber::DrawIntRadius(Double_t b)
297{
298//
2a103154 299// Draw Interaction Radius
f3a04204 300//
301 TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
302 c10->cd();
041f7f97 303 fgWIntRadius->SetParameter(0, b);
304 fgWIntRadius->SetMinimum(0.);
305 fgWIntRadius->Draw();
f3a04204 306}
307
2a103154 308void AliFastGlauber::DrawEnergyDensity()
309{
310//
311// Draw energy density
312//
313 TCanvas *c11 = new TCanvas("c11","Energy Density",400, 10, 600, 700);
314 c11->cd();
041f7f97 315 fgWEnergyDensity->SetMinimum(0.);
316 fgWEnergyDensity->Draw();
2a103154 317}
318
5b3a5a5d 319Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
320{
321//
322// Woods-Saxon Parameterisation
323// as a function of radius
324//
325 Double_t xx = x[0];
326 Double_t r0 = par[0];
327 Double_t d = par[1];
328 Double_t w = par[2];
329 Double_t n = par[3];
330
331 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
2a103154 332
5b3a5a5d 333 return y;
334}
335
336Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
337{
338//
339// Wood Saxon Parameterisation
340// as a function of z and b
341//
342 Double_t bb = x[0];
343 Double_t zz = x[1];
344 Double_t r0 = par[0];
345 Double_t d = par[1];
346 Double_t w = par[2];
347 Double_t n = par[3];
348 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
349 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
2a103154 350
5b3a5a5d 351 return y;
352}
353
354Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
355{
356//
357// Wood Saxon Parameterisation
358// as a function of z for fixed b
359//
360 Double_t bb = par[4];
361 Double_t zz = x[0];
362 Double_t r0 = par[0];
363 Double_t d = par[1];
364 Double_t w = par[2];
365 Double_t n = par[3];
366 Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
367 Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
2a103154 368
5b3a5a5d 369 return y;
370}
371
f86dad79 372Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
5b3a5a5d 373{
374//
375// Thickness function
376//
377 Double_t b = x[0];
041f7f97 378 fgWSz->SetParameter(4, b);
379 Double_t y = 2. * fgWSz->Integral(0., fgBMax);
5b3a5a5d 380 return y;
381}
382
383
384
385Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
386{
387//
388// Kernel for overlap function
389//
390 Double_t b = par[0];
391 Double_t r1 = x[0];
392 Double_t phi = x[1];
393 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
041f7f97 394 Double_t y = r1 * fgWSta->Eval(r1) * fgWSta->Eval(r2);
5b3a5a5d 395 return y;
396}
397
398
f3a04204 399Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
400{
401//
402// Almond shaped interaction region
403//
404 Double_t b = par[0];
405 Double_t xx = x[0] + b/2.;
406 Double_t yy = x[1];
407 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
408 Double_t phi = TMath::ATan2(yy,xx);
409
410 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
411//
412// Interaction probability calculated as product of thicknesses
413//
041f7f97 414 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
f3a04204 415 return y;
416}
417
418Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
419{
420//
421// Average radius at which interaction takes place
422//
423// Radius in the Almond
424 Double_t r = x[0];
425// Impact parameter
426 Double_t b = par[0];
041f7f97 427 fgWAlmond->SetParameter(0, b);
f3a04204 428// Steps in phi
429 Double_t dphi = 2. * TMath::Pi() / 100.;
430// Average over phi
431 Double_t phi = 0.;
432 Double_t y = 0.;
433
434 for (Int_t i = 0; i < 100; i++) {
435 Double_t xx = r * TMath::Cos(phi);
436 Double_t yy = r * TMath::Sin(phi);
041f7f97 437 y += fgWAlmond->Eval(xx,yy);
f3a04204 438 phi += dphi;
439 } // phi loop
440// Result multiplied by Jacobian (2 pi r)
441 return (2. * TMath::Pi() * y * r / 100.);
442}
443
444Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
445{
446//
447// Path Length as a function of phi for interaction point fixed at (0,0)
448//
449//
450// Steps in r
041f7f97 451 const Int_t kNp = 100;
452 const Double_t kDr = fgBMax/Double_t(kNp);
f3a04204 453// Impact parameter
454 Double_t b = par[0];
8de7e046 455// Path Length definition
456 Int_t iopt = Int_t(par[1]);
457
f3a04204 458// Phi direction in Almond
459 Double_t phi0 = x[0];
460 Double_t r = 0.;
461 Double_t rw = 0.;
462 Double_t w = 0.;
463// Step along radial direction phi
041f7f97 464 for (Int_t i = 0; i < kNp; i++) {
f3a04204 465//
466// Transform into target frame
467//
468 Double_t xx = r * TMath::Cos(phi0) + b / 2.;
469 Double_t yy = r * TMath::Sin(phi0);
470 Double_t phi = TMath::ATan2(yy, xx);
471
472 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
473// Radius in projectile frame
474 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
041f7f97 475 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
f3a04204 476
477 rw += y * r;
478 w += y;
041f7f97 479 r += kDr;
f3a04204 480 } // radial steps
481//
b6d0ed6b 482// My length definition (is exact for hard disk)
8de7e046 483 if (!iopt) {
484 return (2. * rw / w);
485 } else {
041f7f97 486 return TMath::Sqrt(2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01));
8de7e046 487 }
f3a04204 488}
489
490Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
491{
492//
493// Path Length as a function of phi
494// Interaction point from random distribution
495//
496//
497// r-steps
498//
041f7f97 499 const Int_t kNp = 100;
500 const Double_t kDr = fgBMax/Double_t(kNp);
f3a04204 501// Number of interactions
041f7f97 502 const Int_t kNpi = Int_t (par[1]);
f3a04204 503
504//
505// Impact parameter
506 Double_t b = par[0];
8de7e046 507// Path Length definition
508 Int_t iopt = Int_t(par[2]);
f3a04204 509// Phi direction
510 Double_t phi0 = x[0];
511
512 printf("phi0 %f \n", phi0);
513
514// Path length
515 Double_t l = 0.;
516
041f7f97 517 for (Int_t in = 0; in < kNpi; in ++) {
f3a04204 518 Double_t rw = 0.;
519 Double_t w = 0.;
520
521 // Interaction point
522 Double_t x0, y0;
041f7f97 523 fgWAlmond->GetRandom2(x0, y0);
f3a04204 524// Initial radius
525 Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
041f7f97 526 Int_t nps = Int_t ((fgBMax - r0)/kDr) - 1;
f3a04204 527
528 Double_t r = 0.;
529// Radial steps
530 for (Int_t i = 0; (i < nps ); i++) {
531
532// Transform into target frame
533 Double_t xx = x0 + r * TMath::Cos(phi0) + b / 2.;
534 Double_t yy = y0 + r * TMath::Sin(phi0);
535 Double_t phi = TMath::ATan2(yy, xx);
536 Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
537// Radius in projectile frame
538 Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
041f7f97 539 Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
f3a04204 540
541 rw += y * r;
542 w += y;
041f7f97 543 r += kDr;
f3a04204 544 } // steps
545// Average over interactions
8de7e046 546 if (!iopt) {
547 l += (2. * rw / w);
548 } else {
041f7f97 549 l+= 2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01);
8de7e046 550 }
f3a04204 551 } // interactions
9db56006 552 if (!iopt)
041f7f97 553 return (l / Double_t(kNpi));
9db56006 554 else
041f7f97 555 return (TMath::Sqrt(l / Double_t(kNpi)));
f3a04204 556}
557
f86dad79 558Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
5b3a5a5d 559{
560//
561// Overlap function
562//
563 Double_t b = x[0];
041f7f97 564 fgWStarfi->SetParameter(0, b);
5b3a5a5d 565/*
566 Double_t al[2];
567 Double_t bl[2];
568 al[0] = 0.;
569 al[1] = 0.;
570 bl[0] = 6.6;
571 bl[1] = TMath::Pi();
572 Double_t err;
573
041f7f97 574 Double_t y = 2. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
5b3a5a5d 575 printf("WStaa: %f %f %f\n", b, y, err);
576*/
577//
578// MC Integration
579//
580 Double_t y = 0;
581 for (Int_t i = 0; i < 100000; i++)
582 {
583 Double_t phi = TMath::Pi() * gRandom->Rndm();
041f7f97 584 Double_t b1 = fgBMax * gRandom->Rndm();
585 y += fgWStarfi->Eval(b1, phi);
5b3a5a5d 586 }
041f7f97 587 y *= 2. * 0.1 * 208. * 208. * TMath::Pi() * fgBMax / 100000.;
5b3a5a5d 588 return y;
589}
590
f86dad79 591Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
5b3a5a5d 592{
593//
594// Geometrical Cross-Section
595//
596 Double_t b = x[0];
041f7f97 597 Double_t taa = fgWStaa->Eval(b);
598 const Double_t kSigma = 55.6; // mbarn
5b3a5a5d 599
041f7f97 600 Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- kSigma * taa)); // fm
5b3a5a5d 601 return y;
602}
603
604
605Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
606{
607//
c2140715 608// Number of binary collisions
5b3a5a5d 609//
610 Double_t b = x[0];
611 Double_t sigma = par[0];
041f7f97 612 Double_t taa = fgWStaa->Eval(b);
5b3a5a5d 613
614 Double_t y = 2. * TMath::Pi() * b * sigma * taa; // fm
615 return y;
616}
617
f86dad79 618Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
5b3a5a5d 619{
620//
c2140715 621// Number of hard processes per event
5b3a5a5d 622//
623 Double_t b = x[0];
041f7f97 624 Double_t y = fgWSbinary->Eval(b)/fgWSgeo->Eval(b);
5b3a5a5d 625 return y;
626}
627
2a103154 628Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
629{
630//
631// Initial energy density as a function of the impact parameter
632//
633 Double_t b = x[0];
634 Double_t rA = par[0];
b6d0ed6b 635//
636// Attention: area of transverse reaction zone in hard-sphere approximation !
637 Double_t saa = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA * rA
638 - b * TMath::Sqrt(rA * rA - b * b/ 4.);
041f7f97 639 Double_t taa = fgWStaa->Eval(b);
2a103154 640
641 return (taa/saa);
642}
643
5b3a5a5d 644void AliFastGlauber::SimulateTrigger(Int_t n)
645{
646 //
647 // Simulates Trigger
648 //
649 TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
650 TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
651 TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
652 TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
653
654 mbtH->SetXTitle("b [fm]");
655 hdtH->SetXTitle("b [fm]");
656 mbmH->SetXTitle("Multiplicity");
657 hdmH->SetXTitle("Multiplicity");
658
659 TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
660 c0->Divide(2,1);
661 TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
662 c1->Divide(1,2);
663
664 //
665 //
666 Init(1);
667 for (Int_t iev = 0; iev < n; iev++)
668 {
669 Float_t b, p, mult;
670 GetRandom(b, p, mult);
671 mbtH->Fill(b,1.);
672 hdtH->Fill(b, p);
673 mbmH->Fill(mult, 1.);
674 hdmH->Fill(mult, p);
675
676 c0->cd(1);
677 mbtH->Draw();
678 c0->cd(2);
679 hdtH->Draw();
680 c0->Update();
681
682 c1->cd(1);
683 mbmH->Draw();
684 c1->cd(2);
685 hdmH->Draw();
686 c1->Update();
687 }
688}
689
690void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
691{
692 //
693 // Gives back a random impact parameter, hard trigger probability and multiplicity
694 //
041f7f97 695 b = fgWSgeo->GetRandom();
696 Float_t mu = fgWSN->Eval(b);
5b3a5a5d 697 p = 1.-TMath::Exp(-mu);
041f7f97 698 mult = 6000./fgWSN->Eval(1.) * mu;
5b3a5a5d 699}
700
c2140715 701void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
702{
703 //
704 // Gives back a random impact parameter bin, and hard trigger decission
705 //
041f7f97 706 Float_t b = fgWSgeo->GetRandom();
707 Float_t mu = fgWSN->Eval(b) * fSigmaHard;
c2140715 708 Float_t p = 1.-TMath::Exp(-mu);
709 if (b < 5.) {
710 bin = 1;
711 } else if (b < 8.6) {
712 bin = 2;
713 } else if (b < 11.2) {
714 bin = 3;
715 } else if (b < 13.2) {
716 bin = 4;
204e011f 717 } else if (b < 15.0) {
c2140715 718 bin = 5;
204e011f 719 } else {
720 bin = 6;
c2140715 721 }
722
723 hard = kFALSE;
724
725 Float_t r = gRandom->Rndm();
726
727 if (r < p) hard = kTRUE;
728}
729
730
5b3a5a5d 731Float_t AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
732{
733 //
734 // Gives back a random impact parameter in the range bmin .. bmax
735 //
736
737 Float_t b = -1.;
738 while(b < bmin || b > bmax)
041f7f97 739 b = fgWSgeo->GetRandom();
5b3a5a5d 740 return b;
741}
742
743Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
744{
745 //
746 // Return cross-section integrated from b1 to b2
747 //
748
041f7f97 749 return fgWSgeo->Integral(b1, b2)/100.;
5b3a5a5d 750}
751
752Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
753{
754 //
755 // Return raction of hard cross-section integrated from b1 to b2
756 //
757
041f7f97 758 return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
5b3a5a5d 759}
760
761
762Double_t AliFastGlauber::Binaries(Double_t b)
763{
764 //
765 // Return number of binary collisions normalized to 1 at b=0
766 //
767
041f7f97 768 return fgWSN->Eval(b)/fgWSN->Eval(0.001);
5b3a5a5d 769}
a2f2f511 770
771//=================== Added by A. Dainese 11/02/04 ===========================
772void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
773{
774 //
775 // Set limits of centrality class as fractions of the geomtrical
776 // cross section
777 //
778 if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
779 printf(" Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
780 return;
781 }
782
783 Double_t bLow=0.,bUp=0.;
784 Double_t xsecFr=0.;
785
786 while(xsecFr<xsecFrLow) {
787 xsecFr = fgWSgeo->Integral(0.,bLow)/fgWSgeo->Integral(0.,100.);
788 bLow += 0.1;
789 }
790 bUp = bLow;
791 while(xsecFr<xsecFrUp) {
792 xsecFr = fgWSgeo->Integral(0.,bUp)/fgWSgeo->Integral(0.,100.);
793 bUp += 0.1;
794 }
795
796 printf(" Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm\n",
797 xsecFrLow,xsecFrUp,bLow,bUp);
798
799 fgWSbinary->SetRange(bLow,bUp);
800
801 return;
802}
803
804void AliFastGlauber::GetRandomBHard(Double_t& b)
805{
806 //
807 // Get random impact parameter according to distribution of
808 // hard (binary) cross-section, in the range defined by the centrality class
809 //
810 b = fgWSbinary->GetRandom();
811 Int_t bin = 2*(Int_t)b;
812 if( (b-(Int_t)b) > 0.5) bin++;
813 fWAlmondCurrent = &fWAlmondFixedB[bin];
814 return;
815}
816
817void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
818{
819 //
820 // Get random position of parton production point according to
821 // product of thickness functions
822 //
823 fWAlmondCurrent->GetRandom2(x,y);
824 return;
825}
826
827void AliFastGlauber::GetRandomPhi(Double_t& phi)
828{
829 //
830 // Get random parton azimuthal propagation direction
831 //
832 phi = 2.*TMath::Pi()*gRandom->Rndm();
833 return;
834}
835
836Double_t AliFastGlauber::CalculateLength(Double_t b,
837 Double_t x0,Double_t y0,Double_t phi0)
838{
839 //
840 // Calculate path length for a parton with production point (x0,y0)
841 // and propagation direction (ux=cos(phi0),uy=sin(phi0))
842 // in a collision with impact parameter b
843 //
844
845 // number of steps in l
846 const Int_t kNp = 100;
847 const Double_t kDl = fgBMax/Double_t(kNp);
848
849 Double_t l,integral,integral1,integral2,r0,xx,yy,phi,r1,r2,ell;
850 Double_t rhocoll,rhocollHalfMax;
851 Int_t i,nps;
852
853 if(fEllDef==1) {
854 //
855 // Definition 1:
856 //
857 // ell = 2 * \int_0^\infty dl*l*\rho_coll(x0+l*ux,y0+l*uy) /
858 // \int_0^\infty dl*\rho_coll(x0+l*ux,y0+l*uy)
859 //
860
861
862 l = 0.;
863 integral1 = 0.;
864 integral2 = 0.;
865
866 // Initial radius
867 r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
868 nps = Int_t ((fgBMax - r0)/kDl) - 1;
869
870 // Radial steps
871 for (i = 0; i < nps; i++) {
872
873 // Transform into target frame
874 xx = x0 + l * TMath::Cos(phi0) + b / 2.;
875 yy = y0 + l * TMath::Sin(phi0);
876 phi = TMath::ATan2(yy, xx);
877 r1 = TMath::Sqrt(xx * xx + yy * yy);
878 // Radius in projectile frame
879 r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
880 rhocoll = fgWSta->Eval(r1) * fgWSta->Eval(r2);
881
882 integral1 += rhocoll * l * kDl;
883 integral2 += rhocoll * kDl;
884 l += kDl;
885 } // steps
886
887 ell = (2. * integral1 / integral2);
888 return ell;
889
890 } else if(fEllDef==2) {
891 //
892 // Definition 2:
893 //
894 // ell = \int_0^\infty dl*
895 // \Theta(\rho_coll(x0+l*ux,y0+l*uy)-0.5*\rho_coll(0,0))
896 //
897
898 l = 0.;
899 integral = 0.;
900
901 // Initial radius
902 r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
903 nps = Int_t ((fgBMax - r0)/kDl) - 1;
904
905 rhocollHalfMax = 0.5*fWAlmondCurrent->Eval(0.,0.);
906
907 // Radial steps
908 for (i = 0; i < nps; i++) {
909
910 // Transform into target frame
911 xx = x0 + l * TMath::Cos(phi0) + b / 2.;
912 yy = y0 + l * TMath::Sin(phi0);
913 phi = TMath::ATan2(yy, xx);
914 r1 = TMath::Sqrt(xx * xx + yy * yy);
915 // Radius in projectile frame
916 r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
917 rhocoll = fgWSta->Eval(r1) * fgWSta->Eval(r2);
918
919 if(rhocoll>rhocollHalfMax) integral += kDl;
920
921 l += kDl;
922 } // steps
923
924 ell = integral;
925 return ell;
926
927 } else {
928 printf(" Wrong length definition setting!\n");
929 return -1.;
930 }
931}
932
933void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
934{
935 //
936 // Return length from random b, x0, y0, phi0
937 //
938 Double_t x0,y0,phi0;
939 if(b<0.) GetRandomBHard(b);
940 GetRandomXY(x0,y0);
941 GetRandomPhi(phi0);
942
943 ell = CalculateLength(b,x0,y0,phi0);
944
945 return;
946}
947
948void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
949 Double_t b)
950{
951 //
952 // Return 2 lengths back to back from random b, x0, y0, phi0
953 //
954 Double_t x0,y0,phi0;
955 if(b<0.) GetRandomBHard(b);
956 GetRandomXY(x0,y0);
957 GetRandomPhi(phi0);
958 Double_t phi0plusPi = phi0+TMath::Pi();
959
960 ell1 = CalculateLength(b,x0,y0,phi0);
961 ell2 = CalculateLength(b,x0,y0,phi0plusPi);
962
963 return;
964}
965
966void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell,
967 Double_t b)
968{
969 //
970 // Returns lenghts for n partons with azimuthal angles phi[n]
971 // from random b, x0, y0
972 //
973 Double_t x0,y0;
974 if(b<0.) GetRandomBHard(b);
975 GetRandomXY(x0,y0);
976
977 for(Int_t i=0; i<n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
978
979 return;
980}
981
982void AliFastGlauber::PlotBDistr(Int_t n)
983{
984 //
985 // Plot distribution of n impact parameters
986 //
987 Double_t b;
988 TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax);
989 hB->SetXTitle("b [fm]");
990 hB->SetYTitle("dN/db [a.u.]");
991 hB->SetFillColor(3);
992
993 for(Int_t i=0; i<n; i++) {
994 GetRandomBHard(b);
995 hB->Fill(b);
996 }
997
998 TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
999 cB->cd();
1000 hB->Draw();
1001
1002 return;
1003}
1004
1005void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname)
1006{
1007 //
1008 // Plot length distribution
1009 //
1010 Double_t ell;
1011 TH1F *hEll = new TH1F("hEll","Length distribution",16,-0.5,15.5);
1012 hEll->SetXTitle("Transverse path length, L [fm]");
1013 hEll->SetYTitle("Probability");
1014 hEll->SetFillColor(2);
1015
1016 for(Int_t i=0; i<n; i++) {
1017 GetLength(ell);
1018 hEll->Fill(ell);
1019 }
1020 hEll->Scale(1/(Double_t)n);
1021
1022 TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
1023 cL->cd();
1024 hEll->Draw();
1025
1026 if(save) {
1027 TFile *f = new TFile(fname,"recreate");
1028 hEll->Write();
1029 f->Close();
1030 }
1031 return;
1032}
1033
1034void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,Char_t *fname)
1035{
1036 //
1037 // Plot lengths back-to-back distributions
1038 //
1039 Double_t ell1,ell2;
1040 TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15);
1041 hElls->SetXTitle("Transverse path length, L [fm]");
1042 hElls->SetYTitle("Transverse path length, L [fm]");
1043
1044 for(Int_t i=0; i<n; i++) {
1045 GetLengthsBackToBack(ell1,ell2);
1046 hElls->Fill(ell1,ell2);
1047 }
1048 hElls->Scale(1/(Double_t)n);
1049
1050
1051 TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
1052 gStyle->SetPalette(1,0);
1053 cLs->cd();
1054 hElls->Draw("col,Z");
1055
1056 if(save) {
1057 TFile *f = new TFile(fname,"recreate");
1058 hElls->Write();
1059 f->Close();
1060 }
1061 return;
1062}
1063
1064void AliFastGlauber::StoreAlmonds()
1065{
1066 //
1067 // Store in glauberPbPb.root 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
1068 //
1069
1070 Char_t almondName[100];
1071 TFile* ff = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root","update");
1072 for(Int_t k=0; k<40; k++) {
1073 sprintf(almondName,"WAlmondFixedB%d",k);
1074 Double_t b = 0.25+k*0.5;
1075 printf(" b = %f\n",b);
1076 fgWAlmond->SetParameter(0,b);
1077
1078 fgWAlmond->Write(almondName);
1079 }
1080 ff->Close();
1081
1082 return;
1083}
1084
1085void AliFastGlauber::PlotAlmonds()
1086{
1087 //
1088 // Plot almonds for some impact parameters
1089 //
1090
1091 TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
1092 gStyle->SetPalette(1,0);
1093 c->Divide(2,2);
1094 c->cd(1);
1095 fWAlmondFixedB[0].Draw("cont1");
1096 c->cd(2);
1097 fWAlmondFixedB[10].Draw("cont1");
1098 c->cd(3);
1099 fWAlmondFixedB[20].Draw("cont1");
1100 c->cd(4);
1101 fWAlmondFixedB[30].Draw("cont1");
1102
1103 return;
1104}