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