Adding include files where needed for latest ROOT
[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//
65aa45f2 18// Utility class to make simple Glauber type calculations
a42548b0 19// for SYMMETRIC collision geometries (AA):
041f7f97 20// Impact parameter, production points, reaction plane dependence
65aa45f2 21//
041f7f97 22// The SimulateTrigger method can be used for simple MB and hard-process
23// (binary scaling) trigger studies.
65aa45f2 24//
041f7f97 25// Some basic quantities can be visualized directly.
041f7f97 26//
65aa45f2 27// The default set-up for PbPb or AUAu collisions can be read from a file
28// calling Init(1) or Init(2) if you want to read Almonds too.
29//
30// ***** If you change settings dont forget to call init afterwards, *****
31// ***** in order to update the formulas with the new parameters. *****
041f7f97 32//
33// Author: andreas.morsch@cern.ch
65aa45f2 34//=================== Added by A. Dainese 11/02/04 ===========================
35// Calculate path length for a parton with production point (x0,y0)
36// and propagation direction (ux=cos(phi0),uy=sin(phi0))
37// in a collision with impact parameter b and functions that make use
38// of it.
39//=================== Added by A. Dainese 05/03/04 ===========================
40// Calculation of line integrals I0 and I1
41// integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
42// integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
43// mostly for use in the Quenching class
44//=================== Added by C. Loizdes 27/03/04 ===========================
45// Handling of AuAu collisions
46// More get/set functions
47// Comments, units and clearing of code
48//
5b3a5a5d 49
50// from AliRoot
51#include "AliFastGlauber.h"
52// from root
9edefa04 53#include <Riostream.h>
54#include <TCanvas.h>
5b3a5a5d 55#include <TF1.h>
56#include <TF2.h>
5b3a5a5d 57#include <TFile.h>
9edefa04 58#include <TH1F.h>
59#include <TH2F.h>
60#include <TLegend.h>
61#include <TMath.h>
65aa45f2 62#include <TROOT.h>
9edefa04 63#include <TRandom.h>
64#include <TStyle.h>
5b3a5a5d 65
66ClassImp(AliFastGlauber)
67
65aa45f2 68Float_t AliFastGlauber::fgBMax = 0.;
041f7f97 69TF1* AliFastGlauber::fgWSb = NULL;
70TF2* AliFastGlauber::fgWSbz = NULL;
71TF1* AliFastGlauber::fgWSz = NULL;
72TF1* AliFastGlauber::fgWSta = NULL;
73TF2* AliFastGlauber::fgWStarfi = NULL;
74TF2* AliFastGlauber::fgWAlmond = NULL;
75TF1* AliFastGlauber::fgWStaa = NULL;
76TF1* AliFastGlauber::fgWSgeo = NULL;
77TF1* AliFastGlauber::fgWSbinary = NULL;
78TF1* AliFastGlauber::fgWSN = NULL;
79TF1* AliFastGlauber::fgWPathLength0 = NULL;
80TF1* AliFastGlauber::fgWPathLength = NULL;
81TF1* AliFastGlauber::fgWEnergyDensity = NULL;
82TF1* AliFastGlauber::fgWIntRadius = NULL;
1bc228f5 83TF2* AliFastGlauber::fgWKParticipants = NULL;
84TF1* AliFastGlauber::fgWParticipants = NULL;
65aa45f2 85TF2* AliFastGlauber::fgWAlmondCurrent = NULL;
7f2f270b 86TF2* AliFastGlauber::fgWAlmondFixedB[40];
65aa45f2 87const Int_t AliFastGlauber::fgkMCInts = 100000;
88Int_t AliFastGlauber::fgCounter = 0;
5b3a5a5d 89
e6e76983 90AliFastGlauber::AliFastGlauber():
91 fWSr0(0.),
92 fWSd(0.),
93 fWSw(0.),
94 fWSn(0.),
95 fSigmaHard(0.),
96 fSigmaNN(0.),
97 fA(0),
98 fBmin(0.),
99 fBmax(0.),
100 fEllDef(0),
101 fName()
5b3a5a5d 102{
a42548b0 103 // Default Constructor
65aa45f2 104 fgCounter++;
105 if(fgCounter>1)
e6e76983 106 Error("AliFastGlauber","More than one instance (%d) is not supported, check your code!",fgCounter);
65aa45f2 107
108 // Defaults for Pb
109 SetMaxImpact();
110 SetLengthDefinition();
111 SetPbPbLHC();
112}
113
a42548b0 114AliFastGlauber::AliFastGlauber(const AliFastGlauber & gl)
e6e76983 115 :TObject(gl),
116 fWSr0(0.),
117 fWSd(0.),
118 fWSw(0.),
119 fWSn(0.),
120 fSigmaHard(0.),
121 fSigmaNN(0.),
122 fA(0),
123 fBmin(0.),
124 fBmax(0.),
125 fEllDef(0),
126 fName()
a42548b0 127{
128// Copy constructor
129 gl.Copy(*this);
130}
131
65aa45f2 132AliFastGlauber::~AliFastGlauber()
133{
a42548b0 134// Destructor
65aa45f2 135 fgCounter--;
7f2f270b 136 for(Int_t k=0; k<40; k++) delete fgWAlmondFixedB[k];
65aa45f2 137}
138
139void AliFastGlauber::SetAuAuRhic()
140{
141 //Set all parameters for RHIC
142 SetWoodSaxonParametersAu();
143 SetHardCrossSection();
710a8d90 144 SetNNCrossSection(42);
65aa45f2 145 SetNucleus(197);
146 SetFileName("$(ALICE_ROOT)/FASTSIM/data/glauberAuAu.root");
147}
148
149void AliFastGlauber::SetPbPbLHC()
150{
151 //Set all parameters for LHC
152 SetWoodSaxonParametersPb();
153 SetHardCrossSection();
154 SetNNCrossSection();
155 SetNucleus();
156 SetFileName();
5b3a5a5d 157}
158
159void AliFastGlauber::Init(Int_t mode)
160{
65aa45f2 161 // Initialisation
162 // mode = 0; all functions are calculated
163 // mode = 1; overlap function is read from file (for Pb-Pb only)
164 // mode = 2; interaction almond functions are read from file
165 // USE THIS FOR PATH LENGTH CALC.!
166 //
a2f2f511 167
65aa45f2 168 //
169 Reset();
170 //
5b3a5a5d 171
65aa45f2 172 //
173 // Wood-Saxon
174 //
175 fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
176 fgWSb->SetParameter(0, fWSr0);
177 fgWSb->SetParameter(1, fWSd);
178 fgWSb->SetParameter(2, fWSw);
179 fgWSb->SetParameter(3, fWSn);
180
fac5662b 181 fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 0, fgBMax, 4);
65aa45f2 182 fgWSbz->SetParameter(0, fWSr0);
183 fgWSbz->SetParameter(1, fWSd);
184 fgWSbz->SetParameter(2, fWSw);
185 fgWSbz->SetParameter(3, fWSn);
186
187 fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
188 fgWSz->SetParameter(0, fWSr0);
189 fgWSz->SetParameter(1, fWSd);
190 fgWSz->SetParameter(2, fWSw);
191 fgWSz->SetParameter(3, fWSn);
192
193 //
194 // Thickness
195 //
196 fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
5b3a5a5d 197
65aa45f2 198 //
199 // Overlap Kernel
200 //
201 fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
202 fgWStarfi->SetParameter(0, 0.);
203 fgWStarfi->SetNpx(200);
204 fgWStarfi->SetNpy(20);
8de7e046 205
65aa45f2 206 //
207 // Participants Kernel
208 //
209 fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 3);
210 fgWKParticipants->SetParameter(0, 0.);
211 fgWKParticipants->SetParameter(1, fSigmaNN);
212 fgWKParticipants->SetParameter(2, fA);
213 fgWKParticipants->SetNpx(200);
214 fgWKParticipants->SetNpy(20);
f3a04204 215
65aa45f2 216 //
217 // Overlap and Participants
218 //
219 if (!mode) {
220 fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 1);
221 fgWStaa->SetNpx(100);
222 fgWStaa->SetParameter(0,fA);
223 fgWStaa->SetNpx(100);
224 fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 2);
225 fgWParticipants->SetParameter(0, fSigmaNN);
226 fgWParticipants->SetParameter(1, fA);
227 fgWParticipants->SetNpx(100);
228 } else {
229 Info("Init","Reading overlap function from file %s",fName.Data());
230 TFile* f = new TFile(fName.Data());
231 if(!f->IsOpen()){
232 Fatal("Init", "Could not open file %s",fName.Data());
233 }
234 fgWStaa = (TF1*) f->Get("WStaa");
235 fgWParticipants = (TF1*) f->Get("WParticipants");
236 delete f;
237 }
2a103154 238
65aa45f2 239 //
240 // Energy Density
241 //
242 fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
243 fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
244
245 //
246 // Geometrical Cross-Section
247 //
248 fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 1);
249 fgWSgeo->SetParameter(0,fSigmaNN); //mbarn
250 fgWSgeo->SetNpx(100);
2a103154 251
65aa45f2 252 //
253 // Hard cross section (binary collisions)
254 //
255 fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
256 fgWSbinary->SetParameter(0, fSigmaHard); //mbarn
257 fgWSbinary->SetNpx(100);
1bc228f5 258
65aa45f2 259 //
260 // Hard collisions per event
261 //
262 fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
263 fgWSN->SetNpx(100);
264
265 //
266 // Almond shaped interaction region
267 //
268 fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
269 fgWAlmond->SetParameter(0, 0.);
270 fgWAlmond->SetNpx(200);
271 fgWAlmond->SetNpy(200);
272
273 if(mode==2) {
274 Info("Init","Reading interaction almonds from file: %s",fName.Data());
275 Char_t almondName[100];
276 TFile* ff = new TFile(fName.Data());
277 for(Int_t k=0; k<40; k++) {
278 sprintf(almondName,"WAlmondFixedB%d",k);
279 fgWAlmondCurrent = (TF2*)ff->Get(almondName);
7f2f270b 280 fgWAlmondFixedB[k] = fgWAlmondCurrent;
5b3a5a5d 281 }
65aa45f2 282 delete ff;
283 }
284
285 fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
286 fgWIntRadius->SetParameter(0, 0.);
287
288 //
289 // Path Length as a function of Phi
290 //
291 fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
292 fgWPathLength0->SetParameter(0, 0.);
293 fgWPathLength0->SetParameter(1, 0.); //Pathlength definition
294
295 fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
296 fgWPathLength->SetParameter(0, 0.); //Impact Parameter
297 fgWPathLength->SetParameter(1, 1000.); //Number of interactions used for average
298 fgWPathLength->SetParameter(2, 0); //Pathlength definition
299}
300
a42548b0 301void AliFastGlauber::Reset() const
65aa45f2 302{
303 //
304 // Reset dynamic allocated formulas
305 // in case init is called twice
306
307 if(fgWSb) delete fgWSb;
308 if(fgWSbz) delete fgWSbz;
309 if(fgWSz) delete fgWSz;
310 if(fgWSta) delete fgWSta;
311 if(fgWStarfi) delete fgWStarfi;
312 if(fgWAlmond) delete fgWAlmond;
313 if(fgWStaa) delete fgWStaa;
314 if(fgWSgeo) delete fgWSgeo;
315 if(fgWSbinary) delete fgWSbinary;
316 if(fgWSN) delete fgWSN;
317 if(fgWPathLength0) delete fgWPathLength0;
318 if(fgWPathLength) delete fgWPathLength;
319 if(fgWEnergyDensity) delete fgWEnergyDensity;
320 if(fgWIntRadius) delete fgWIntRadius;
321 if(fgWKParticipants) delete fgWKParticipants;
322 if(fgWParticipants) delete fgWParticipants;
5b3a5a5d 323}
324
710a8d90 325void AliFastGlauber::DrawWSb() const
5b3a5a5d 326{
65aa45f2 327 //
328 // Draw Wood-Saxon Nuclear Density Function
329 //
330 TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
331 c1->cd();
332 Double_t max=fgWSb->GetMaximum(0,fgBMax)*1.01;
333 TH2F *h2f=new TH2F("h2fwsb","Wood Saxon: #rho(r) = n (1-#omega(r/r_{0})^2)/(1+exp((r-r_{0})/d)) [fm^{-3}]",2,0,fgBMax,2,0,max);
334 h2f->SetStats(0);
335 h2f->GetXaxis()->SetTitle("r [fm]");
336 h2f->GetYaxis()->SetNoExponent(kTRUE);
337 h2f->GetYaxis()->SetTitle("#rho [fm^{-3}]");
338 h2f->Draw();
339 fgWSb->Draw("same");
340 TLegend *l1a = new TLegend(0.45,0.6,.90,0.8);
341 l1a->SetFillStyle(0);
342 l1a->SetBorderSize(0);
343 Char_t label[100];
344 sprintf(label,"r_{0} = %.2f fm",fWSr0);
345 l1a->AddEntry(fgWSb,label,"");
346 sprintf(label,"d = %.2f fm",fWSd);
347 l1a->AddEntry(fgWSb,label,"");
348 sprintf(label,"n = %.2e fm^{-3}",fWSn);
349 l1a->AddEntry(fgWSb,label,"");
350 sprintf(label,"#omega = %.2f",fWSw);
351 l1a->AddEntry(fgWSb,label,"");
352 l1a->Draw();
353 c1->Update();
5b3a5a5d 354}
355
710a8d90 356void AliFastGlauber::DrawOverlap() const
5b3a5a5d 357{
65aa45f2 358 //
359 // Draw Overlap Function
360 //
361 TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
362 c2->cd();
363 Double_t max=fgWStaa->GetMaximum(0,fgBMax)*1.01;
a42548b0 364 TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0, max);
65aa45f2 365 h2f->SetStats(0);
366 h2f->GetXaxis()->SetTitle("b [fm]");
367 h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]");
368 h2f->Draw();
369 fgWStaa->Draw("same");
5b3a5a5d 370}
371
710a8d90 372void AliFastGlauber::DrawParticipants() const
1bc228f5 373{
65aa45f2 374 //
375 // Draw Number of Participants Npart
376 //
377 TCanvas *c3 = new TCanvas("c3","Participants",400,10,600,700);
378 c3->cd();
379 Double_t max=fgWParticipants->GetMaximum(0,fgBMax)*1.01;
380 TH2F *h2f=new TH2F("h2fpart","Number of Participants",2,0,fgBMax,2,0,max);
381 h2f->SetStats(0);
382 h2f->GetXaxis()->SetTitle("b [fm]");
383 h2f->GetYaxis()->SetTitle("N_{part}");
384 h2f->Draw();
385 fgWParticipants->Draw("same");
386 TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
387 l1a->SetFillStyle(0);
388 l1a->SetBorderSize(0);
389 Char_t label[100];
390 sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
391 l1a->AddEntry(fgWParticipants,label,"");
392 l1a->Draw();
393 c3->Update();
1bc228f5 394}
395
710a8d90 396void AliFastGlauber::DrawThickness() const
5b3a5a5d 397{
65aa45f2 398 //
399 // Draw Thickness Function
400 //
401 TCanvas *c4 = new TCanvas("c4","Thickness",400,10,600,700);
402 c4->cd();
403 Double_t max=fgWSta->GetMaximum(0,fgBMax)*1.01;
404 TH2F *h2f=new TH2F("h2fta","Thickness function: T_{A} [fm^{-2}]",2,0,fgBMax,2,0,max);
405 h2f->SetStats(0);
406 h2f->GetXaxis()->SetTitle("b [fm]");
407 h2f->GetYaxis()->SetTitle("T_{A} [fm^{-2}]");
408 h2f->Draw();
409 fgWSta->Draw("same");
5b3a5a5d 410}
411
710a8d90 412void AliFastGlauber::DrawGeo() const
5b3a5a5d 413{
65aa45f2 414 //
415 // Draw Geometrical Cross-Section
416 //
417 TCanvas *c5 = new TCanvas("c5","Geometrical Cross-Section",400,10,600,700);
418 c5->cd();
419 Double_t max=fgWSgeo->GetMaximum(0,fgBMax)*1.01;
420 TH2F *h2f=new TH2F("h2fgeo","Differential Geometrical Cross-Section: d#sigma^{geo}_{AB}/db [fm]",2,0,fgBMax,2,0,max);
421 h2f->SetStats(0);
422 h2f->GetXaxis()->SetTitle("b [fm]");
423 h2f->GetYaxis()->SetTitle("d#sigma^{geo}_{AB}/db [fm]");
424 h2f->Draw();
425 fgWSgeo->Draw("same");
426 TLegend *l1a = new TLegend(0.10,0.8,.40,0.9);
427 l1a->SetFillStyle(0);
428 l1a->SetBorderSize(0);
429 Char_t label[100];
430 sprintf(label,"#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
431 l1a->AddEntry(fgWSgeo,label,"");
432 l1a->Draw();
433 c5->Update();
5b3a5a5d 434}
435
710a8d90 436void AliFastGlauber::DrawBinary() const
5b3a5a5d 437{
65aa45f2 438 //
439 // Draw Binary Cross-Section
440 //
441 TCanvas *c6 = new TCanvas("c6","Binary Cross-Section",400,10,600,700);
442 c6->cd();
443 Double_t max=fgWSbinary->GetMaximum(0,fgBMax)*1.01;
444 TH2F *h2f=new TH2F("h2fbinary","Differential Binary Cross-Section: #sigma^{hard}_{NN} dT_{AB}/db [fm]",2,0,fgBMax,2,0,max);
445 h2f->SetStats(0);
446 h2f->GetXaxis()->SetTitle("b [fm]");
447 h2f->GetYaxis()->SetTitle("d#sigma^{hard}_{AB}/db [fm]");
448 h2f->Draw();
449 fgWSbinary->Draw("same");
450 TLegend *l1a = new TLegend(0.50,0.8,.90,0.9);
451 l1a->SetFillStyle(0);
452 l1a->SetBorderSize(0);
453 Char_t label[100];
454 sprintf(label,"#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
455 l1a->AddEntry(fgWSb,label,"");
456 l1a->Draw();
457 c6->Update();
5b3a5a5d 458}
459
710a8d90 460void AliFastGlauber::DrawN() const
5b3a5a5d 461{
65aa45f2 462 //
463 // Draw Binaries per event (Ncoll)
464 //
465 TCanvas *c7 = new TCanvas("c7","Binaries per event",400,10,600,700);
466 c7->cd();
467 Double_t max=fgWSN->GetMaximum(0,fgBMax)*1.01;
468 TH2F *h2f=new TH2F("h2fhardcols","Number of hard collisions: T_{AB} #sigma^{hard}_{NN}/#sigma_{AB}^{geo}",2,0,fgBMax,2,0,max);
469 h2f->SetStats(0);
470 h2f->GetXaxis()->SetTitle("b [fm]");
471 h2f->GetYaxis()->SetTitle("N_{coll}");
472 h2f->Draw();
473 fgWSN->Draw("same");
474 TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
475 l1a->SetFillStyle(0);
476 l1a->SetBorderSize(0);
477 Char_t label[100];
478 sprintf(label,"#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
479 l1a->AddEntry(fgWSN,label,"");
480 sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
481 l1a->AddEntry(fgWSN,label,"");
482 l1a->Draw();
483 c7->Update();
5b3a5a5d 484}
485
710a8d90 486void AliFastGlauber::DrawKernel(Double_t b) const
5b3a5a5d 487{
65aa45f2 488 //
489 // Draw Kernel
490 //
491 TCanvas *c8 = new TCanvas("c8","Kernel",400,10,600,700);
492 c8->cd();
493 fgWStarfi->SetParameter(0, b);
494 TH2F *h2f=new TH2F("h2fkernel","Kernel of Overlap function: d^{2}T_{AB}/dr/d#phi [fm^{-3}]",2,0,fgBMax,2,0,TMath::Pi());
495 h2f->SetStats(0);
496 h2f->GetXaxis()->SetTitle("r [fm]");
497 h2f->GetYaxis()->SetTitle("#phi [rad]");
498 h2f->Draw();
499 fgWStarfi->Draw("same");
500 TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
501 l1a->SetFillStyle(0);
502 l1a->SetBorderSize(0);
503 Char_t label[100];
a0e187b1 504 sprintf(label,"b = %.1f fm",b);
65aa45f2 505 l1a->AddEntry(fgWStarfi,label,"");
506 l1a->Draw();
507 c8->Update();
5b3a5a5d 508}
509
710a8d90 510void AliFastGlauber::DrawAlmond(Double_t b) const
f3a04204 511{
65aa45f2 512 //
513 // Draw Interaction Almond
514 //
515 TCanvas *c9 = new TCanvas("c9","Almond",400,10,600,700);
516 c9->cd();
517 fgWAlmond->SetParameter(0, b);
518 TH2F *h2f=new TH2F("h2falmond","Interaction Almond [fm^{-4}]",2,0,fgBMax,2,0,fgBMax);
519 h2f->SetStats(0);
520 h2f->GetXaxis()->SetTitle("x [fm]");
521 h2f->GetYaxis()->SetTitle("y [fm]");
522 h2f->Draw();
523 fgWAlmond->Draw("same");
524 TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
525 l1a->SetFillStyle(0);
526 l1a->SetBorderSize(0);
527 Char_t label[100];
a0e187b1 528 sprintf(label,"b = %.1f fm",b);
65aa45f2 529 l1a->AddEntry(fgWAlmond,label,"");
530 l1a->Draw();
531 c9->Update();
f3a04204 532}
533
710a8d90 534void AliFastGlauber::DrawEnergyDensity() const
f3a04204 535{
65aa45f2 536 //
537 // Draw energy density
538 //
539 TCanvas *c10 = new TCanvas("c10","Energy Density",400, 10, 600, 700);
540 c10->cd();
541 fgWEnergyDensity->SetMinimum(0.);
542 Double_t max=fgWEnergyDensity->GetMaximum(0,fgWEnergyDensity->GetParameter(0))*1.01;
543 TH2F *h2f=new TH2F("h2fenergydens","Energy density",2,0,fgBMax,2,0,max);
544 h2f->SetStats(0);
545 h2f->GetXaxis()->SetTitle("b [fm]");
546 h2f->GetYaxis()->SetTitle("fm^{-4}");
547 h2f->Draw();
548 fgWEnergyDensity->Draw("same");
549 c10->Update();
f3a04204 550}
551
710a8d90 552void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt) const
f3a04204 553{
65aa45f2 554 //
555 // Draw Path Length
556 //
557 TCanvas *c11 = new TCanvas("c11","Path Length",400,10,600,700);
558 c11->cd();
559 fgWPathLength0->SetParameter(0, b);
560 fgWPathLength0->SetParameter(1, Double_t(iopt));
561 fgWPathLength0->SetMinimum(0.);
562 fgWPathLength0->SetMaximum(10.);
563 TH2F *h2f=new TH2F("h2fpathlength0","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
564 h2f->SetStats(0);
565 h2f->GetXaxis()->SetTitle("#phi [rad]");
566 h2f->GetYaxis()->SetTitle("l [fm]");
567 h2f->Draw();
568 fgWPathLength0->Draw("same");
f3a04204 569}
570
710a8d90 571void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt) const
f3a04204 572{
65aa45f2 573 //
574 // Draw Path Length
575 //
576 TCanvas *c12 = new TCanvas("c12","Path Length",400,10,600,700);
577 c12->cd();
578 fgWAlmond->SetParameter(0, b);
579 fgWPathLength->SetParameter(0, b);
580 fgWPathLength->SetParameter(1, Double_t (ni));
581 fgWPathLength->SetParameter(2, Double_t (iopt));
582 fgWPathLength->SetMinimum(0.);
583 fgWPathLength->SetMaximum(10.);
584 TH2F *h2f=new TH2F("h2fpathlength","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
585 h2f->SetStats(0);
586 h2f->GetXaxis()->SetTitle("#phi [rad]");
587 h2f->GetYaxis()->SetTitle("l [fm]");
588 h2f->Draw();
589 fgWPathLength->Draw("same");
f3a04204 590}
591
710a8d90 592void AliFastGlauber::DrawIntRadius(Double_t b) const
2a103154 593{
65aa45f2 594 //
595 // Draw Interaction Radius
596 //
597 TCanvas *c13 = new TCanvas("c13","Interaction Radius",400,10,600,700);
598 c13->cd();
599 fgWIntRadius->SetParameter(0, b);
600 fgWIntRadius->SetMinimum(0);
601 Double_t max=fgWIntRadius->GetMaximum(0,fgBMax)*1.01;
602 TH2F *h2f=new TH2F("h2fintradius","Interaction Density",2,0.,fgBMax,2,0,max);
603 h2f->SetStats(0);
604 h2f->GetXaxis()->SetTitle("r [fm]");
605 h2f->GetYaxis()->SetTitle("[fm^{-3}]");
606 h2f->Draw();
607 fgWIntRadius->Draw("same");
2a103154 608}
609
5b3a5a5d 610Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
611{
65aa45f2 612 //
613 // Woods-Saxon Parameterisation
614 // as a function of radius (xx)
615 //
710a8d90 616 const Double_t kxx = x[0]; //fm
617 const Double_t kr0 = par[0]; //fm
618 const Double_t kd = par[1]; //fm
619 const Double_t kw = par[2]; //no units
620 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
621 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
65aa45f2 622 return y; //fm^-3
5b3a5a5d 623}
624
625Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
626{
65aa45f2 627 //
628 // Wood Saxon Parameterisation
629 // as a function of z and b
630 //
710a8d90 631 const Double_t kbb = x[0]; //fm
632 const Double_t kzz = x[1]; //fm
633 const Double_t kr0 = par[0]; //fm
634 const Double_t kd = par[1]; //fm
635 const Double_t kw = par[2]; //no units
636 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
637 const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
638 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
65aa45f2 639 return y; //fm^-3
5b3a5a5d 640}
641
642Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
643{
65aa45f2 644 //
645 // Wood Saxon Parameterisation
646 // as a function of z for fixed b
647 //
710a8d90 648 const Double_t kzz = x[0]; //fm
649 const Double_t kr0 = par[0]; //fm
650 const Double_t kd = par[1]; //fm
651 const Double_t kw = par[2]; //no units
652 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
653 const Double_t kbb = par[4]; //fm
654 const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
655 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
65aa45f2 656 return y; //fm^-3
5b3a5a5d 657}
658
f86dad79 659Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
5b3a5a5d 660{
65aa45f2 661 //
662 // Thickness function T_A
663 // as a function of b
664 //
710a8d90 665 const Double_t kb = x[0];
666 fgWSz->SetParameter(4, kb);
667 Double_t y = 2. * fgWSz->Integral(0., fgBMax);
65aa45f2 668 return y; //fm^-2
5b3a5a5d 669}
670
5b3a5a5d 671Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
672{
65aa45f2 673 //
674 // Kernel for overlap function: T_A(s)*T_A(s-b)
675 // as a function of r and phi
710a8d90 676 const Double_t kr1 = x[0];
677 const Double_t kphi = x[1];
678 const Double_t kb = par[0];
679 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
680 Double_t y = kr1 * fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
65aa45f2 681 return y; //fm^-3
682}
683
684Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
685{
686 //
687 // Overlap function
688 // T_{AB}=Int d2s T_A(s)*T_B(s-b)
689 // as a function of b
690 // (normalized to fA*fB)
691 //
710a8d90 692 const Double_t kb = x[0];
693 const Double_t ka = par[0];
694 fgWStarfi->SetParameter(0, kb);
65aa45f2 695
696 // root integration seems to fail
697 /*
698 Double_t al[2];
699 Double_t bl[2];
700 al[0] = 1e-6;
701 al[1] = fgBMax;
702 bl[0] = 0.;
703 bl[1] = TMath::Pi();
704 Double_t err;
705
706 Double_t y = 2. * 208. * 208. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
707 printf("WStaa: %.5e %.5e %.5e\n", b, y, err);
708 */
709
710 //
711 // MC Integration
712 //
713 Double_t y = 0;
a42548b0 714
715
65aa45f2 716 for (Int_t i = 0; i < fgkMCInts; i++)
717 {
a42548b0 718
710a8d90 719 const Double_t kphi = TMath::Pi() * gRandom->Rndm();
720 const Double_t kb1 = fgBMax * gRandom->Rndm();
721 y += fgWStarfi->Eval(kb1, kphi);
65aa45f2 722 }
723 y *= 2. * TMath::Pi() * fgBMax / fgkMCInts; //fm^-2
710a8d90 724 y *= ka * ka * 0.1; //mbarn^-1
65aa45f2 725 return y;
5b3a5a5d 726}
727
1bc228f5 728Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par)
729{
65aa45f2 730 //
731 // Kernel for number of participants
732 // as a function of r and phi
733 //
710a8d90 734 const Double_t kr1 = x[0];
735 const Double_t kphi = x[1];
736 const Double_t kb = par[0]; //fm
737 const Double_t ksig = par[1]; //mbarn
738 const Double_t ka = par[2]; //mass number
739 const Double_t kr2 = TMath::Sqrt(kr1*kr1 +kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
740 const Double_t kxsi = fgWSta->Eval(kr2) * ksig * 0.1; //no units
65aa45f2 741 /*
710a8d90 742 Double_t y=(1-TMath::Power((1-xsi),aa))
65aa45f2 743 */
710a8d90 744 Double_t a = ka;
745 Double_t sum = ka * kxsi;
65aa45f2 746 Double_t y = sum;
710a8d90 747 for (Int_t i = 1; i <= ka; i++)
1bc228f5 748 {
65aa45f2 749 a--;
710a8d90 750 sum *= (-kxsi) * a / Float_t(i+1);
65aa45f2 751 y += sum;
1bc228f5 752 }
710a8d90 753 y *= kr1 * fgWSta->Eval(kr1);
65aa45f2 754 return y; //fm^-1
1bc228f5 755}
756
65aa45f2 757Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* par)
758{
759 //
760 // Number of Participants as
761 // a function of b
762 //
710a8d90 763 const Double_t kb = x[0];
764 const Double_t ksig = par[0]; //mbarn
765 const Double_t ka = par[1]; //mass number
766 fgWKParticipants->SetParameter(0, kb);
767 fgWKParticipants->SetParameter(1, ksig);
768 fgWKParticipants->SetParameter(2, ka);
65aa45f2 769
770 //
771 // MC Integration
772 //
773 Double_t y = 0;
774 for (Int_t i = 0; i < fgkMCInts; i++)
775 {
710a8d90 776 const Double_t kphi = TMath::Pi() * gRandom->Rndm();
777 const Double_t kb1 = fgBMax * gRandom->Rndm();
778 y += fgWKParticipants->Eval(kb1, kphi);
65aa45f2 779 }
710a8d90 780 y *= 2. * ka * 2. * TMath::Pi() * fgBMax / fgkMCInts;
65aa45f2 781 return y; //no units
782}
783
784Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par)
785{
786 //
787 // Geometrical Cross-Section
788 // as a function of b
789 //
710a8d90 790 const Double_t kb = x[0]; //fm
791 const Double_t ksigNN = par[0]; //mbarn
792 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
793 Double_t y = 2. * TMath::Pi() * kb * (1. - TMath::Exp(- ksigNN * ktaa));
65aa45f2 794 return y; //fm
795}
796
797Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
798{
799 //
710a8d90 800 // Number of binary hard collisions
65aa45f2 801 // as a function of b
802 //
710a8d90 803 const Double_t kb = x[0]; //fm
804 const Double_t ksig = par[0]; //mbarn
805 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
806 Double_t y = 2. * TMath::Pi() * kb * ksig * ktaa;
65aa45f2 807 return y; //fm
808}
809
810Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
811{
812 //
813 // Number of hard processes per event
814 // as a function of b
710a8d90 815 const Double_t kb = x[0];
816 Double_t y = fgWSbinary->Eval(kb)/fgWSgeo->Eval(kb);
65aa45f2 817 return y; //no units
818}
819
820Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
821{
822 //
823 // Initial energy density
824 // as a function of the impact parameter
825 //
710a8d90 826 const Double_t kb = x[0];
827 const Double_t krA = par[0];
65aa45f2 828 //
829 // Attention: area of transverse reaction zone in hard-sphere approximation !
710a8d90 830 const Double_t krA2=krA*krA;
831 const Double_t kb2=kb*kb;
832 const Double_t ksaa = (TMath::Pi() - 2. * TMath::ASin(kb/ 2./ krA)) * krA2
833 - kb * TMath::Sqrt(krA2 - kb2/ 4.); //fm^2
834 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
835 Double_t y=ktaa/ksaa*10;
65aa45f2 836 return y; //fm^-4
837}
5b3a5a5d 838
f3a04204 839Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
840{
65aa45f2 841 //
842 // Almond shaped interaction region
843 // as a function of cartesian x,y.
844 //
710a8d90 845 const Double_t kb = par[0];
846 const Double_t kxx = x[0] + kb/2.;
847 const Double_t kyy = x[1];
848 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
849 const Double_t kphi = TMath::ATan2(kyy,kxx);
850 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
65aa45f2 851 //
852 // Interaction probability calculated as product of thicknesses
853 //
710a8d90 854 Double_t y = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
65aa45f2 855 return y; //fm^-4
f3a04204 856}
857
858Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
859{
65aa45f2 860 //
861 // Average interaction density over radius
862 // at which interaction takes place
863 // as a function of radius
864 //
710a8d90 865 const Double_t kr = x[0];
866 const Double_t kb = par[0];
867 fgWAlmond->SetParameter(0, kb);
65aa45f2 868 // Average over phi in small steps
710a8d90 869 const Double_t kdphi = 2. * TMath::Pi() / 100.;
65aa45f2 870 Double_t phi = 0.;
871 Double_t y = 0.;
872 for (Int_t i = 0; i < 100; i++) {
710a8d90 873 const Double_t kxx = kr * TMath::Cos(phi);
874 const Double_t kyy = kr * TMath::Sin(phi);
875 y += fgWAlmond->Eval(kxx,kyy);
876 phi += kdphi;
65aa45f2 877 } // phi loop
878 // Result multiplied by Jacobian (2 pi r)
710a8d90 879 y *= 2. * TMath::Pi() * kr / 100.;
65aa45f2 880 return y; //fm^-3
f3a04204 881}
882
883Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
884{
65aa45f2 885 //
886 // Path Length as a function of phi
887 // for interaction point fixed at (0,0)
888 // as a function of phi-direction
889 //
890 // Phi direction in Almond
710a8d90 891 const Double_t kphi0 = x[0];
892 const Double_t kb = par[0];
65aa45f2 893 // Path Length definition
710a8d90 894 const Int_t kiopt = Int_t(par[1]);
65aa45f2 895
896 // Step along radial direction phi
897 const Int_t kNp = 100; // Steps in r
898 const Double_t kDr = fgBMax/kNp;
899 Double_t r = 0.;
900 Double_t rw = 0.;
901 Double_t w = 0.;
902 for (Int_t i = 0; i < kNp; i++) {
903 //
904 // Transform into target frame
905 //
710a8d90 906 const Double_t kxx = r * TMath::Cos(kphi0) + kb / 2.;
907 const Double_t kyy = r * TMath::Sin(kphi0);
908 const Double_t kphi = TMath::ATan2(kyy, kxx);
909 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
65aa45f2 910 // Radius in projectile frame
710a8d90 911 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
912 const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
65aa45f2 913
710a8d90 914 rw += ky * r;
915 w += ky;
65aa45f2 916 r += kDr;
917 } // radial steps
918
919 Double_t y=0.;
710a8d90 920 if (!kiopt) // My length definition (is exact for hard disk)
65aa45f2 921 if(w) y= 2. * rw / w;
922 else {
710a8d90 923 const Double_t knorm=fgWSta->Eval(1e-4);
924 if(knorm) y = TMath::Sqrt(2. * rw * kDr / knorm / knorm);
65aa45f2 925 }
926 return y; //fm
927}
928
929Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
930{
931 //
932 // Path Length as a function of phi
933 // Interaction point from random distribution
934 // as a function of the phi-direction
710a8d90 935 const Double_t kphi0 = x[0];
936 const Double_t kb = par[0];
937 fgWAlmond->SetParameter(0, kb);
938 const Int_t kNpi = Int_t (par[1]); //Number of interactions
939 const Int_t kiopt = Int_t(par[2]); //Path Length definition
65aa45f2 940
941 //
942 // r-steps
943 //
944 const Int_t kNp = 100;
945 const Double_t kDr = fgBMax/Double_t(kNp);
946 Double_t l = 0.; // Path length
947 for (Int_t in = 0; in < kNpi; in ++) {
f3a04204 948 Double_t rw = 0.;
949 Double_t w = 0.;
65aa45f2 950 // Interaction point
951 Double_t x0, y0;
952 fgWAlmond->GetRandom2(x0, y0);
953 // Initial radius
710a8d90 954 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
955 const Int_t knps = Int_t ((fgBMax - kr0)/kDr) - 1;
f3a04204 956
65aa45f2 957 // Radial steps
958 Double_t r = 0.;
710a8d90 959 for (Int_t i = 0; (i < knps ); i++) {
65aa45f2 960 // Transform into target frame
710a8d90 961 const Double_t kxx = x0 + r * TMath::Cos(kphi0) + kb / 2.;
962 const Double_t kyy = y0 + r * TMath::Sin(kphi0);
963 const Double_t kphi = TMath::ATan2(kyy, kxx);
964 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
65aa45f2 965 // Radius in projectile frame
710a8d90 966 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
967 const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
65aa45f2 968
710a8d90 969 rw += ky * r;
970 w += ky;
65aa45f2 971 r += kDr;
972 } // steps
973 // Average over interactions
710a8d90 974 if (!kiopt) {
65aa45f2 975 if(w) l += (2. * rw / w);
8de7e046 976 } else {
710a8d90 977 const Double_t knorm=fgWSta->Eval(1e-4);
978 if(knorm) l+= 2. * rw * kDr / knorm / knorm;
8de7e046 979 }
65aa45f2 980 } // interactions
981 Double_t ret=0;
710a8d90 982 if (!kiopt)
65aa45f2 983 ret= l / kNpi;
984 else
985 ret=TMath::Sqrt( l / kNpi);
986 return ret; //fm
f3a04204 987}
988
710a8d90 989Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2) const
f3a04204 990{
65aa45f2 991 //
992 // Return the geometrical cross-section integrated from b1 to b2
993 //
994 return fgWSgeo->Integral(b1, b2)*10.; //mbarn
995}
f3a04204 996
710a8d90 997Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const
5b3a5a5d 998{
65aa45f2 999 //
1000 // Return the hard cross-section integrated from b1 to b2
1001 //
1002 return fgWSbinary->Integral(b1, b2)*10.; //mbarn
5b3a5a5d 1003}
1004
710a8d90 1005Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) const
1bc228f5 1006{
65aa45f2 1007 //
148c5ce5 1008 // Return fraction of hard cross-section integrated from b1 to b2
65aa45f2 1009 //
1010 return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
1bc228f5 1011}
1012
f762082f 1013Double_t AliFastGlauber::NHard(Double_t b1, Double_t b2) const
1014{
1015 //
1016 // Number of binary hard collisions
1017 // as a function of b (nucl/ex/0302016 eq. 19)
1018 //
1019 const Double_t kshard=HardCrossSection(b1,b2);
1020 const Double_t ksgeo=CrossSection(b1,b2);
1021 if(ksgeo>0)
1022 return kshard/ksgeo;
1023 else return -1;
1024}
1025
710a8d90 1026Double_t AliFastGlauber::Binaries(Double_t b) const
5b3a5a5d 1027{
65aa45f2 1028 //
710a8d90 1029 // Return number of binary hard collisions normalized to 1 at b=0
65aa45f2 1030 //
1031 if(b==0) b=1e-4;
1032 return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
5b3a5a5d 1033}
1034
a42548b0 1035Double_t AliFastGlauber::MeanOverlap(Double_t b1, Double_t b2)
1036{
1037//
1038// Calculate the mean overlap for impact parameter range b1 .. b2
1039//
1040 Double_t sum = 0.;
1041 Double_t sumc = 0.;
1042 Double_t b = b1;
1043
1044 while (b < b2-0.005) {
1045 Double_t nc = GetNumberOfCollisions(b);
1046 sum += 10. * fgWStaa->Eval(b) * fgWSgeo->Eval(b) * 0.01 / (1. - TMath::Exp(-nc));
1047 sumc += 10. * fgWSgeo->Eval(b) * 0.01;
1048 b += 0.01;
1049 }
1050 return (sum / CrossSection(b1, b2));
1051}
1052
1053
1054Double_t AliFastGlauber::MeanNumberOfCollisionsPerEvent(Double_t b1, Double_t b2)
1055{
1056//
1057// Calculate the mean number of collisions per event for impact parameter range b1 .. b2
1058//
1059 Double_t sum = 0.;
1060 Double_t sumc = 0.;
1061 Double_t b = b1;
1062
1063 while (b < b2-0.005) {
1064 Double_t nc = GetNumberOfCollisions(b);
1065 sum += nc / (1. - TMath::Exp(-nc)) * 10. * fgWSgeo->Eval(b) * 0.01;
1066 sumc += 10. * fgWSgeo->Eval(b) * 0.01;
1067 b += 0.01;
1068 }
1069 return (sum / CrossSection(b1, b2));
1070}
1071
1072
710a8d90 1073Double_t AliFastGlauber::GetNumberOfBinaries(Double_t b) const
5b3a5a5d 1074{
65aa45f2 1075 //
710a8d90 1076 // Return number of binary hard collisions at b
65aa45f2 1077 //
1078 if(b==0) b=1e-4;
1079 return fgWSN->Eval(b);
5b3a5a5d 1080}
1081
710a8d90 1082Double_t AliFastGlauber::Participants(Double_t b) const
5b3a5a5d 1083{
65aa45f2 1084 //
1085 // Return the number of participants normalized to 1 at b=0
1086 //
1087 if(b==0) b=1e-4;
1088 return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4));
5b3a5a5d 1089}
1090
710a8d90 1091Double_t AliFastGlauber::GetNumberOfParticipants(Double_t b) const
2a103154 1092{
65aa45f2 1093 //
1094 // Return the number of participants for impact parameter b
1095 //
1096 if(b==0) b=1e-4;
1097 return (fgWParticipants->Eval(b));
2a103154 1098}
1099
710a8d90 1100Double_t AliFastGlauber::GetNumberOfCollisions(Double_t b) const
1101{
1102 //
1103 // Return the number of collisions for impact parameter b
1104 //
1105 if(b==0) b=1e-4;
1106 return (fgWStaa->Eval(b)*fSigmaNN);
1107}
1108
148c5ce5 1109Double_t AliFastGlauber::GetNumberOfCollisionsPerEvent(Double_t b) const
1110{
1111 //
1112 // Return the number of collisions per event (at least one collision)
1113 // for impact parameter b
1114 //
1115 Double_t n = GetNumberOfCollisions(b);
1116 if (n > 0.) {
1117 return (n / (1. - TMath::Exp(- n)));
1118 } else {
1119 return (0.);
1120 }
1121}
1122
5b3a5a5d 1123void AliFastGlauber::SimulateTrigger(Int_t n)
1124{
65aa45f2 1125 //
1126 // Simulates Trigger
1127 //
1128 TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
1129 TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
1130 TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
1131 TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
5b3a5a5d 1132
65aa45f2 1133 mbtH->SetXTitle("b [fm]");
1134 hdtH->SetXTitle("b [fm]");
1135 mbmH->SetXTitle("Multiplicity");
1136 hdmH->SetXTitle("Multiplicity");
5b3a5a5d 1137
65aa45f2 1138 TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
1139 c0->Divide(2,1);
1140 TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
1141 c1->Divide(1,2);
5b3a5a5d 1142
65aa45f2 1143 //
1144 //
1145 Init(1);
1146 for (Int_t iev = 0; iev < n; iev++)
5b3a5a5d 1147 {
65aa45f2 1148 Float_t b, p, mult;
1149 GetRandom(b, p, mult);
1150 mbtH->Fill(b,1.);
1151 hdtH->Fill(b, p);
1152 mbmH->Fill(mult, 1.);
1153 hdmH->Fill(mult, p);
1154
1155 c0->cd(1);
1156 mbtH->Draw();
1157 c0->cd(2);
1158 hdtH->Draw();
1159 c0->Update();
1160
1161 c1->cd(1);
1162 mbmH->Draw();
1163 c1->cd(2);
1164 hdmH->Draw();
1165 c1->Update();
5b3a5a5d 1166 }
1167}
1168
1169void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
1170{
65aa45f2 1171 //
1172 // Gives back a random impact parameter, hard trigger probability and multiplicity
1173 //
1174 b = fgWSgeo->GetRandom();
710a8d90 1175 const Float_t kmu = fgWSN->Eval(b);
1176 p = 1.-TMath::Exp(-kmu);
1177 mult = 6000./fgWSN->Eval(1.) * kmu;
5b3a5a5d 1178}
1179
c2140715 1180void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
1181{
65aa45f2 1182 //
1183 // Gives back a random impact parameter bin, and hard trigger decission
1184 //
710a8d90 1185 const Float_t kb = fgWSgeo->GetRandom();
1186 const Float_t kmu = fgWSN->Eval(kb) * fSigmaHard;
1187 const Float_t kp = 1.-TMath::Exp(-kmu);
1188 if (kb < 5.) {
65aa45f2 1189 bin = 1;
710a8d90 1190 } else if (kb < 8.6) {
65aa45f2 1191 bin = 2;
710a8d90 1192 } else if (kb < 11.2) {
65aa45f2 1193 bin = 3;
710a8d90 1194 } else if (kb < 13.2) {
65aa45f2 1195 bin = 4;
710a8d90 1196 } else if (kb < 15.0) {
65aa45f2 1197 bin = 5;
1198 } else {
1199 bin = 6;
1200 }
1201 hard = kFALSE;
710a8d90 1202 const Float_t kr = gRandom->Rndm();
1203 if (kr < kp) hard = kTRUE;
c2140715 1204}
1205
65aa45f2 1206Double_t AliFastGlauber::GetRandomImpactParameter(Double_t bmin, Double_t bmax)
5b3a5a5d 1207{
65aa45f2 1208 //
1209 // Gives back a random impact parameter in the range bmin .. bmax
1210 //
1211 Float_t b = -1.;
1212 while(b < bmin || b > bmax)
1213 b = fgWSgeo->GetRandom();
1214 return b;
5b3a5a5d 1215}
1216
710a8d90 1217void AliFastGlauber::StoreFunctions() const
5b3a5a5d 1218{
65aa45f2 1219 //
1220 // Store in file functions
1221 //
1222 TFile* ff = new TFile(fName.Data(),"recreate");
1223 fgWStaa->Write("WStaa");
1224 fgWParticipants->Write("WParticipants");
1225 ff->Close();
1226 return;
5b3a5a5d 1227}
1228
65aa45f2 1229//=================== Added by A. Dainese 11/02/04 ===========================
1bc228f5 1230
710a8d90 1231void AliFastGlauber::StoreAlmonds() const
5b3a5a5d 1232{
65aa45f2 1233 //
1234 // Store in file
1235 // 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
1236 //
1237 Char_t almondName[100];
1238 TFile* ff = new TFile(fName.Data(),"update");
1239 for(Int_t k=0; k<40; k++) {
1240 sprintf(almondName,"WAlmondFixedB%d",k);
1241 Double_t b = 0.25+k*0.5;
1242 Info("StoreAlmonds"," b = %f\n",b);
1243 fgWAlmond->SetParameter(0,b);
1244 fgWAlmond->Write(almondName);
1245 }
1246 ff->Close();
1247 return;
5b3a5a5d 1248}
a2f2f511 1249
a2f2f511 1250void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
1251{
1252 //
65aa45f2 1253 // Set limits of centrality class as fractions
1254 // of the geomtrical cross section
a2f2f511 1255 //
1256 if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
65aa45f2 1257 Error("SetCentralityClass", "Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
a2f2f511 1258 return;
1259 }
1260
1261 Double_t bLow=0.,bUp=0.;
1262 Double_t xsecFr=0.;
710a8d90 1263 const Double_t knorm=fgWSgeo->Integral(0.,100.);
a2f2f511 1264 while(xsecFr<xsecFrLow) {
710a8d90 1265 xsecFr = fgWSgeo->Integral(0.,bLow)/knorm;
a2f2f511 1266 bLow += 0.1;
1267 }
1268 bUp = bLow;
1269 while(xsecFr<xsecFrUp) {
710a8d90 1270 xsecFr = fgWSgeo->Integral(0.,bUp)/knorm;
a2f2f511 1271 bUp += 0.1;
1272 }
1273
710a8d90 1274 Info("SetCentralityClass", "Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm",
a2f2f511 1275 xsecFrLow,xsecFrUp,bLow,bUp);
a2f2f511 1276 fgWSbinary->SetRange(bLow,bUp);
710a8d90 1277 fBmin=bLow;
1278 fBmax=bUp;
a2f2f511 1279 return;
1280}
1281
1282void AliFastGlauber::GetRandomBHard(Double_t& b)
1283{
1284 //
1285 // Get random impact parameter according to distribution of
1286 // hard (binary) cross-section, in the range defined by the centrality class
1287 //
1288 b = fgWSbinary->GetRandom();
1289 Int_t bin = 2*(Int_t)b;
1290 if( (b-(Int_t)b) > 0.5) bin++;
7f2f270b 1291 fgWAlmondCurrent = fgWAlmondFixedB[bin];
a2f2f511 1292 return;
1293}
1294
1295void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
1296{
1297 //
1298 // Get random position of parton production point according to
1299 // product of thickness functions
1300 //
65aa45f2 1301 fgWAlmondCurrent->GetRandom2(x,y);
a2f2f511 1302 return;
1303}
1304
1305void AliFastGlauber::GetRandomPhi(Double_t& phi)
1306{
1307 //
1308 // Get random parton azimuthal propagation direction
1309 //
1310 phi = 2.*TMath::Pi()*gRandom->Rndm();
1311 return;
1312}
1313
65aa45f2 1314Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Double_t phi0)
a2f2f511 1315{
1316 //
1317 // Calculate path length for a parton with production point (x0,y0)
1318 // and propagation direction (ux=cos(phi0),uy=sin(phi0))
1319 // in a collision with impact parameter b
1320 //
1321
1322 // number of steps in l
1323 const Int_t kNp = 100;
1324 const Double_t kDl = fgBMax/Double_t(kNp);
1325
a2f2f511 1326 if(fEllDef==1) {
1327 //
1328 // Definition 1:
1329 //
65aa45f2 1330 // ell = 2 * \int_0^\infty dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy) /
1331 // \int_0^\infty dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
a2f2f511 1332 //
a2f2f511 1333
a2f2f511 1334 // Initial radius
710a8d90 1335 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1336 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
65aa45f2 1337 Double_t l = 0.;
1338 Double_t integral1 = 0.;
1339 Double_t integral2 = 0.;
a2f2f511 1340 // Radial steps
710a8d90 1341 for (Int_t i = 0; i < knps; i++) {
a2f2f511 1342
1343 // Transform into target frame
710a8d90 1344 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1345 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1346 const Double_t kphi = TMath::ATan2(kyy, kxx);
1347 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
a2f2f511 1348 // Radius in projectile frame
710a8d90 1349 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1350 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
a2f2f511 1351
710a8d90 1352 integral1 += kprodTATB * l * kDl;
1353 integral2 += kprodTATB * kDl;
a2f2f511 1354 l += kDl;
1355 } // steps
1356
65aa45f2 1357 Double_t ell=0.;
1358 if(integral2)
1359 ell = (2. * integral1 / integral2);
a2f2f511 1360 return ell;
a2f2f511 1361 } else if(fEllDef==2) {
1362 //
1363 // Definition 2:
1364 //
1365 // ell = \int_0^\infty dl*
65aa45f2 1366 // \Theta((T_A*T_B)(x0+l*ux,y0+l*uy)-0.5*(T_A*T_B)(0,0))
a2f2f511 1367 //
1368
a2f2f511 1369 // Initial radius
710a8d90 1370 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1371 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1372 const Double_t kprodTATBHalfMax = 0.5*fgWAlmondCurrent->Eval(0.,0.);
a2f2f511 1373 // Radial steps
65aa45f2 1374 Double_t l = 0.;
1375 Double_t integral = 0.;
710a8d90 1376 for (Int_t i = 0; i < knps; i++) {
a2f2f511 1377 // Transform into target frame
710a8d90 1378 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1379 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1380 const Double_t kphi = TMath::ATan2(kyy, kxx);
1381 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
a2f2f511 1382 // Radius in projectile frame
710a8d90 1383 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1384 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1385 if(kprodTATB>kprodTATBHalfMax) integral += kDl;
a2f2f511 1386 l += kDl;
1387 } // steps
65aa45f2 1388 Double_t ell = integral;
a2f2f511 1389 return ell;
a2f2f511 1390 } else {
65aa45f2 1391 Error("CalculateLength","Wrong length definition setting: %d !\n",fEllDef);
a2f2f511 1392 return -1.;
1393 }
1394}
1395
83f67d08 1396void AliFastGlauber::GetLengthAndPhi(Double_t& ell,Double_t& phi,Double_t b)
a2f2f511 1397{
1398 //
1399 // Return length from random b, x0, y0, phi0
83f67d08 1400 // Return also phi0
a2f2f511 1401 //
1402 Double_t x0,y0,phi0;
1403 if(b<0.) GetRandomBHard(b);
1404 GetRandomXY(x0,y0);
1405 GetRandomPhi(phi0);
83f67d08 1406 phi = phi0;
a2f2f511 1407 ell = CalculateLength(b,x0,y0,phi0);
a2f2f511 1408 return;
1409}
1410
83f67d08 1411void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
a2f2f511 1412{
1413 //
83f67d08 1414 // Return length from random b, x0, y0, phi0
a2f2f511 1415 //
83f67d08 1416 Double_t phi;
1417 GetLengthAndPhi(ell,phi,b);
1418 return;
1419}
1420
1421void AliFastGlauber::GetLengthsBackToBackAndPhi(Double_t& ell1,Double_t& ell2,Double_t &phi,Double_t b)
1422{
1423 //
1424 // Return 2 lengths back to back from random b, x0, y0, phi0
1425 // Return also phi0
1426 //
a2f2f511 1427 Double_t x0,y0,phi0;
1428 if(b<0.) GetRandomBHard(b);
1429 GetRandomXY(x0,y0);
1430 GetRandomPhi(phi0);
710a8d90 1431 const Double_t kphi0plusPi = phi0+TMath::Pi();
83f67d08 1432 phi = phi0;
a2f2f511 1433 ell1 = CalculateLength(b,x0,y0,phi0);
710a8d90 1434 ell2 = CalculateLength(b,x0,y0,kphi0plusPi);
a2f2f511 1435 return;
1436}
1437
83f67d08 1438void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
1439 Double_t b)
1440{
1441 //
1442 // Return 2 lengths back to back from random b, x0, y0, phi0
1443 //
1444 Double_t phi;
1445 GetLengthsBackToBackAndPhi(ell1,ell2,phi,b);
1446 return;
1447}
1448
1bc228f5 1449void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Double_t b)
a2f2f511 1450{
1451 //
1452 // Returns lenghts for n partons with azimuthal angles phi[n]
1453 // from random b, x0, y0
1454 //
1bc228f5 1455 Double_t x0, y0;
1456 if(b < 0.) GetRandomBHard(b);
a2f2f511 1457 GetRandomXY(x0,y0);
1bc228f5 1458 for(Int_t i = 0; i< n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
a2f2f511 1459 return;
1460}
1461
1462void AliFastGlauber::PlotBDistr(Int_t n)
65aa45f2 1463{
a2f2f511 1464 //
1465 // Plot distribution of n impact parameters
1466 //
1467 Double_t b;
1468 TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax);
1469 hB->SetXTitle("b [fm]");
1470 hB->SetYTitle("dN/db [a.u.]");
1471 hB->SetFillColor(3);
a2f2f511 1472 for(Int_t i=0; i<n; i++) {
1473 GetRandomBHard(b);
1474 hB->Fill(b);
1475 }
a2f2f511 1476 TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
1477 cB->cd();
1478 hB->Draw();
a2f2f511 1479 return;
1480}
1481
d3d4a92f 1482void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,const char *fname)
a2f2f511 1483{
1484 //
1485 // Plot length distribution
1486 //
1487 Double_t ell;
710a8d90 1488 TH1F *hEll = new TH1F("hEll","Length distribution",64,-0.5,15);
a2f2f511 1489 hEll->SetXTitle("Transverse path length, L [fm]");
1490 hEll->SetYTitle("Probability");
1491 hEll->SetFillColor(2);
a2f2f511 1492 for(Int_t i=0; i<n; i++) {
1493 GetLength(ell);
1494 hEll->Fill(ell);
1495 }
1496 hEll->Scale(1/(Double_t)n);
a2f2f511 1497 TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
1498 cL->cd();
1499 hEll->Draw();
1500
1501 if(save) {
1502 TFile *f = new TFile(fname,"recreate");
1503 hEll->Write();
1504 f->Close();
1505 }
1506 return;
1507}
1508
d3d4a92f 1509void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,const char *fname)
a2f2f511 1510{
1511 //
1512 // Plot lengths back-to-back distributions
1513 //
1514 Double_t ell1,ell2;
1515 TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15);
1516 hElls->SetXTitle("Transverse path length, L [fm]");
1517 hElls->SetYTitle("Transverse path length, L [fm]");
a2f2f511 1518 for(Int_t i=0; i<n; i++) {
1519 GetLengthsBackToBack(ell1,ell2);
1520 hElls->Fill(ell1,ell2);
1521 }
1522 hElls->Scale(1/(Double_t)n);
a2f2f511 1523 TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
1524 gStyle->SetPalette(1,0);
1525 cLs->cd();
1526 hElls->Draw("col,Z");
a2f2f511 1527 if(save) {
1528 TFile *f = new TFile(fname,"recreate");
1529 hElls->Write();
1530 f->Close();
1531 }
1532 return;
1533}
1534
710a8d90 1535void AliFastGlauber::PlotAlmonds() const
1536{
1537 //
1538 // Plot almonds for some impact parameters
1539 //
1540 TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
1541 gStyle->SetPalette(1,0);
1542 c->Divide(2,2);
1543 c->cd(1);
7f2f270b 1544 fgWAlmondFixedB[0]->Draw("cont1");
710a8d90 1545 c->cd(2);
7f2f270b 1546 fgWAlmondFixedB[10]->Draw("cont1");
710a8d90 1547 c->cd(3);
7f2f270b 1548 fgWAlmondFixedB[20]->Draw("cont1");
710a8d90 1549 c->cd(4);
7f2f270b 1550 fgWAlmondFixedB[30]->Draw("cont1");
710a8d90 1551 return;
1552}
1553
65aa45f2 1554//=================== Added by A. Dainese 05/03/04 ===========================
1555
1556void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1,
1557 Double_t b,Double_t x0,Double_t y0,
710a8d90 1558 Double_t phi0,Double_t ellCut) const
a2f2f511 1559{
65aa45f2 1560 //
1561 // Calculate integrals:
e9663638 1562 // integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
1563 // integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
a2f2f511 1564 //
65aa45f2 1565 // for a parton with production point (x0,y0)
1566 // and propagation direction (ux=cos(phi0),uy=sin(phi0))
1567 // in a collision with impact parameter b
1568 //
1569
1570 // number of steps in l
1571 const Int_t kNp = 100;
1572 const Double_t kDl = fgBMax/Double_t(kNp);
1573
1574 // Initial radius
710a8d90 1575 const Double_t kr0 = TMath::Sqrt(x0 * x0 + y0 * y0);
1576 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
65aa45f2 1577
1578 // Radial steps
1579 Double_t l = 0.;
1580 integral0 = 0.;
1581 integral1 = 0.;
1582 Int_t i = 0;
710a8d90 1583 while((i < knps) && (l < ellCut)) {
65aa45f2 1584 // Transform into target frame
710a8d90 1585 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1586 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1587 const Double_t kphi = TMath::ATan2(kyy, kxx);
1588 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
65aa45f2 1589 // Radius in projectile frame
710a8d90 1590 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1591 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1592 integral0 += kprodTATB * kDl;
1593 integral1 += kprodTATB * l * kDl;
65aa45f2 1594 l += kDl;
1595 i++;
1596 } // steps
1597 return;
1598}
1599
83f67d08 1600void AliFastGlauber::GetI0I1AndPhi(Double_t& integral0,Double_t& integral1,
1601 Double_t& phi,
1602 Double_t ellCut,Double_t b)
65aa45f2 1603{
a2f2f511 1604 //
65aa45f2 1605 // Return I0 and I1 from random b, x0, y0, phi0
83f67d08 1606 // Return also phi
65aa45f2 1607 //
1608 Double_t x0,y0,phi0;
1609 if(b<0.) GetRandomBHard(b);
1610 GetRandomXY(x0,y0);
1611 GetRandomPhi(phi0);
83f67d08 1612 phi = phi0;
65aa45f2 1613 CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut);
1614 return;
1615}
a2f2f511 1616
83f67d08 1617void AliFastGlauber::GetI0I1(Double_t& integral0,Double_t& integral1,
1618 Double_t ellCut,Double_t b)
1619{
1620 //
1621 // Return I0 and I1 from random b, x0, y0, phi0
1622 //
1623 Double_t phi;
1624 GetI0I1AndPhi(integral0,integral1,phi,ellCut,b);
1625 return;
1626}
1627
1628void AliFastGlauber::GetI0I1BackToBackAndPhi(Double_t& integral01,Double_t& integral11,
1629 Double_t& integral02,Double_t& integral12,
1630 Double_t& phi,
1631 Double_t ellCut,Double_t b)
65aa45f2 1632{
1633 //
1634 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
83f67d08 1635 // Return also phi0
65aa45f2 1636 //
1637 Double_t x0,y0,phi0;
1638 if(b<0.) GetRandomBHard(b);
1639 GetRandomXY(x0,y0);
1640 GetRandomPhi(phi0);
83f67d08 1641 phi = phi0;
710a8d90 1642 const Double_t kphi0plusPi = phi0+TMath::Pi();
c54404bf 1643 CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
1644 CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
1645 return;
1646}
1647
1648void AliFastGlauber::GetI0I1BackToBackAndPhiAndXY(Double_t& integral01,Double_t& integral11,
1649 Double_t& integral02,Double_t& integral12,
1650 Double_t& phi,Double_t &x,Double_t &y,
1651 Double_t ellCut,Double_t b)
1652{
1653 //
1654 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1655 // Return also phi0
1656 //
1657 Double_t x0,y0,phi0;
1658 if(b<0.) GetRandomBHard(b);
1659 GetRandomXY(x0,y0);
1660 GetRandomPhi(phi0);
1661 phi = phi0; x=x0; y=y0;
1662 const Double_t kphi0plusPi = phi0+TMath::Pi();
65aa45f2 1663 CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
710a8d90 1664 CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
65aa45f2 1665 return;
1666}
a2f2f511 1667
83f67d08 1668void AliFastGlauber::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11,
1669 Double_t& integral02,Double_t& integral12,
1670 Double_t ellCut,Double_t b)
1671{
1672 //
1673 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1674 //
1675 Double_t phi;
1676 GetI0I1BackToBackAndPhi(integral01,integral11,integral02,integral12,
1677 phi,ellCut,b);
1678 return;
1679}
1680
65aa45f2 1681void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
1682 Double_t* integral0,Double_t* integral1,
1683 Double_t ellCut,Double_t b)
1684{
1685 //
1686 // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
1687 // from random b, x0, y0
1688 //
1689 Double_t x0,y0;
1690 if(b<0.) GetRandomBHard(b);
1691 GetRandomXY(x0,y0);
1692 for(Int_t i=0; i<n; i++)
1693 CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1694 return;
1695}
1696
2e3b5c95 1697void AliFastGlauber::GetI0I1ForPythiaAndXY(Int_t n,Double_t* phi,
1698 Double_t* integral0,Double_t* integral1,
1699 Double_t &x,Double_t& y,
1700 Double_t ellCut,Double_t b)
1701{
1702 //
1703 // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
1704 // from random b, x0, y0 and return x0,y0
1705 //
1706 Double_t x0,y0;
1707 if(b<0.) GetRandomBHard(b);
1708 GetRandomXY(x0,y0);
1709 for(Int_t i=0; i<n; i++)
1710 CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1711 x=x0;
1712 y=y0;
1713 return;
1714}
1715
65aa45f2 1716void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
d3d4a92f 1717 Bool_t save,const char *fname)
65aa45f2 1718{
1719 //
bbf8513d 1720 // Plot I0-I1 distribution
65aa45f2 1721 //
1722 Double_t i0,i1;
bbf8513d 1723 TH2F *hI0I1s = new TH2F("hI0I1s","I_{0} versus I_{1}",1000,0,0.001,1000,0,0.01);
1724 hI0I1s->SetXTitle("I_{0} [fm^{-3}]");
1725 hI0I1s->SetYTitle("I_{1} [fm^{-2}]");
1726
65aa45f2 1727 TH1F *hI0 = new TH1F("hI0","I_{0} = #hat{q}L / k",
bbf8513d 1728 1000,0,0.001);
65aa45f2 1729 hI0->SetXTitle("I_{0} [fm^{-3}]");
1730 hI0->SetYTitle("Probability");
1731 hI0->SetFillColor(3);
1732 TH1F *hI1 = new TH1F("hI1","I_{1} = #omega_{c} / k",
bbf8513d 1733 1000,0,0.01);
65aa45f2 1734 hI1->SetXTitle("I_{1} [fm^{-2}]");
1735 hI1->SetYTitle("Probability");
1736 hI1->SetFillColor(4);
1737 TH1F *h2 = new TH1F("h2","2 I_{1}^{2}/I_{0} = R / k",
1738 100,0,0.02);
1739 h2->SetXTitle("2 I_{1}^{2}/I_{0} [fm^{-1}]");
1740 h2->SetYTitle("Probability");
1741 h2->SetFillColor(2);
1742 TH1F *h3 = new TH1F("h3","2 I_{1}/I_{0} = L",
1743 100,0,15);
1744 h3->SetXTitle("2 I_{1}/I_{0} [fm]");
1745 h3->SetYTitle("Probability");
1746 h3->SetFillColor(5);
1747 TH1F *h4 = new TH1F("h4","I_{0}^{2}/(2 I_{1}) = #hat{q} / k",
1748 100,0,0.00015);
1749 h4->SetXTitle("I_{0}^{2}/(2 I_{1}) [fm^{-4}]");
1750 h4->SetYTitle("Probability");
1751 h4->SetFillColor(7);
1752
1753 for(Int_t i=0; i<n; i++) {
1754 GetI0I1(i0,i1,ellCut);
bbf8513d 1755 hI0I1s->Fill(i0,i1);
65aa45f2 1756 hI0->Fill(i0);
1757 hI1->Fill(i1);
1758 h2->Fill(2.*i1*i1/i0);
1759 h3->Fill(2.*i1/i0);
1760 h4->Fill(i0*i0/2./i1);
a2f2f511 1761 }
65aa45f2 1762 hI0->Scale(1/(Double_t)n);
1763 hI1->Scale(1/(Double_t)n);
bbf8513d 1764 h2->Scale(1/(Double_t)n);
1765 h3->Scale(1/(Double_t)n);
1766 h4->Scale(1/(Double_t)n);
1767 hI0I1s->Scale(1/(Double_t)n);
65aa45f2 1768
1769 TCanvas *cI0I1 = new TCanvas("cI0I1","I0 and I1",0,0,900,700);
1770 cI0I1->Divide(3,2);
1771 cI0I1->cd(1);
1772 hI0->Draw();
1773 cI0I1->cd(2);
1774 hI1->Draw();
1775 cI0I1->cd(3);
1776 h2->Draw();
1777 cI0I1->cd(4);
1778 h3->Draw();
1779 cI0I1->cd(5);
1780 h4->Draw();
bbf8513d 1781 cI0I1->cd(6);
1782 gStyle->SetPalette(1,0);
1783 hI0I1s->Draw("col,Z");
a2f2f511 1784
65aa45f2 1785 if(save) {
1786 TFile *f = new TFile(fname,"recreate");
bbf8513d 1787 hI0I1s->Write();
65aa45f2 1788 hI0->Write();
1789 hI1->Write();
1790 h2->Write();
1791 h3->Write();
1792 h4->Write();
1793 f->Close();
1794 }
a2f2f511 1795 return;
1796}
1797
65aa45f2 1798void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
d3d4a92f 1799 Bool_t save,const char *fname)
a2f2f511 1800{
1801 //
bbf8513d 1802 // Plot I0-I1 back-to-back distributions
a2f2f511 1803 //
65aa45f2 1804 Double_t i01,i11,i02,i12;
1805 TH2F *hI0s = new TH2F("hI0s","I_{0}'s back-to-back",100,0,100,100,0,100);
1806 hI0s->SetXTitle("I_{0} [fm^{-3}]");
1807 hI0s->SetYTitle("I_{0} [fm^{-3}]");
1808 TH2F *hI1s = new TH2F("hI1s","I_{1}'s back-to-back",100,0,100,100,0,100);
1809 hI1s->SetXTitle("I_{1} [fm^{-2}]");
1810 hI1s->SetYTitle("I_{1} [fm^{-2}]");
a2f2f511 1811
65aa45f2 1812 for(Int_t i=0; i<n; i++) {
1813 GetI0I1BackToBack(i01,i11,i02,i12,ellCut);
1814 hI0s->Fill(i01,i02);
1815 hI1s->Fill(i11,i12);
1816 }
1817 hI0s->Scale(1/(Double_t)n);
1818 hI1s->Scale(1/(Double_t)n);
1819
1820 TCanvas *cI0I1s = new TCanvas("cI0I1s","I0 and I1 back-to-back distributions",0,0,800,400);
a2f2f511 1821 gStyle->SetPalette(1,0);
65aa45f2 1822 cI0I1s->Divide(2,1);
1823 cI0I1s->cd(1);
1824 hI0s->Draw("col,Z");
1825 cI0I1s->cd(2);
1826 hI1s->Draw("col,Z");
a2f2f511 1827
65aa45f2 1828 if(save) {
1829 TFile *f = new TFile(fname,"recreate");
1830 hI0s->Write();
1831 hI1s->Write();
1832 f->Close();
1833 }
a2f2f511 1834 return;
1835}
710a8d90 1836
a42548b0 1837AliFastGlauber& AliFastGlauber::operator=(const AliFastGlauber& rhs)
1838{
1839// Assignment operator
1840 rhs.Copy(*this);
1841 return *this;
1842}
1843
1844void AliFastGlauber::Copy(TObject&) const
1845{
1846 //
1847 // Copy
1848 //
1849 Fatal("Copy","Not implemented!\n");
1850}
1851