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