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