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