1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 // Utility class to make simple Glauber type calculations
19 // for SYMMETRIC collision geometries (AA):
20 // Impact parameter, production points, reaction plane dependence
22 // The SimulateTrigger method can be used for simple MB and hard-process
23 // (binary scaling) trigger studies.
25 // Some basic quantities can be visualized directly.
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.
30 // ***** If you change settings dont forget to call init afterwards, *****
31 // ***** in order to update the formulas with the new parameters. *****
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
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
51 #include "AliFastGlauber.h"
64 ClassImp(AliFastGlauber)
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;
90 AliFastGlauber::AliFastGlauber():
103 // Default Constructor
106 SetLengthDefinition();
109 fI0I1[0] = fI0I1[1] = 0;
112 AliFastGlauber::AliFastGlauber(const AliFastGlauber & gl)
129 fI0I1[0] = fI0I1[1] = 0;
132 AliFastGlauber* AliFastGlauber::Instance()
134 // Set random number generator
138 fgGlauber = new AliFastGlauber();
143 AliFastGlauber::~AliFastGlauber()
146 for(Int_t k=0; k<40; k++) delete fgWAlmondFixedB[k];
149 void AliFastGlauber::SetAuAuRhic()
151 //Set all parameters for RHIC
152 SetWoodSaxonParametersAu();
153 SetHardCrossSection();
154 SetNNCrossSection(42);
156 SetFileName("$(ALICE_ROOT)/FASTSIM/data/glauberAuAu.root");
159 void AliFastGlauber::SetPbPbLHC()
161 //Set all parameters for LHC
162 SetWoodSaxonParametersPb();
163 SetHardCrossSection();
169 void AliFastGlauber::Init(Int_t mode)
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.!
185 fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
186 fgWSb->SetParameter(0, fWSr0);
187 fgWSb->SetParameter(1, fWSd);
188 fgWSb->SetParameter(2, fWSw);
189 fgWSb->SetParameter(3, fWSn);
191 fgRWSb = new TF1("RWSb", RWSb, 0, fgBMax, 4);
192 fgRWSb->SetParameter(0, fWSr0);
193 fgRWSb->SetParameter(1, fWSd);
194 fgRWSb->SetParameter(2, fWSw);
195 fgRWSb->SetParameter(3, fWSn);
197 fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 0, fgBMax, 4);
198 fgWSbz->SetParameter(0, fWSr0);
199 fgWSbz->SetParameter(1, fWSd);
200 fgWSbz->SetParameter(2, fWSw);
201 fgWSbz->SetParameter(3, fWSn);
203 fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
204 fgWSz->SetParameter(0, fWSr0);
205 fgWSz->SetParameter(1, fWSd);
206 fgWSz->SetParameter(2, fWSw);
207 fgWSz->SetParameter(3, fWSn);
212 fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
217 fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
218 fgWStarfi->SetParameter(0, 0.);
219 fgWStarfi->SetNpx(200);
220 fgWStarfi->SetNpy(20);
223 // Participants Kernel
225 fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 3);
226 fgWKParticipants->SetParameter(0, 0.);
227 fgWKParticipants->SetParameter(1, fSigmaNN);
228 fgWKParticipants->SetParameter(2, fA);
229 fgWKParticipants->SetNpx(200);
230 fgWKParticipants->SetNpy(20);
233 // Overlap and Participants
236 fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 1);
237 fgWStaa->SetNpx(100);
238 fgWStaa->SetParameter(0,fA);
239 fgWStaa->SetNpx(100);
240 fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 2);
241 fgWParticipants->SetParameter(0, fSigmaNN);
242 fgWParticipants->SetParameter(1, fA);
243 fgWParticipants->SetNpx(100);
245 Info("Init","Reading overlap function from file %s",fName.Data());
246 TFile* f = new TFile(fName.Data());
248 Fatal("Init", "Could not open file %s",fName.Data());
250 fgWStaa = (TF1*) f->Get("WStaa");
251 fgWParticipants = (TF1*) f->Get("WParticipants");
258 fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
259 fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
262 // Geometrical Cross-Section
264 fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 1);
265 fgWSgeo->SetParameter(0,fSigmaNN); //mbarn
266 fgWSgeo->SetNpx(100);
269 // Hard cross section (binary collisions)
271 fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
272 fgWSbinary->SetParameter(0, fSigmaHard); //mbarn
273 fgWSbinary->SetNpx(100);
276 // Hard collisions per event
278 fgWSN = new TF1("WSN", WSN, 0.01, fgBMax, 1);
282 // Almond shaped interaction region
284 fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
285 fgWAlmond->SetParameter(0, 0.);
286 fgWAlmond->SetNpx(200);
287 fgWAlmond->SetNpy(200);
290 Info("Init","Reading interaction almonds from file: %s",fName.Data());
291 Char_t almondName[100];
292 TFile* ff = new TFile(fName.Data());
293 for(Int_t k=0; k<40; k++) {
294 snprintf(almondName,100, "WAlmondFixedB%d",k);
295 fgWAlmondCurrent = (TF2*)ff->Get(almondName);
296 fgWAlmondFixedB[k] = fgWAlmondCurrent;
301 fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
302 fgWIntRadius->SetParameter(0, 0.);
305 // Path Length as a function of Phi
307 fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
308 fgWPathLength0->SetParameter(0, 0.);
309 fgWPathLength0->SetParameter(1, 0.); //Pathlength definition
311 fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
312 fgWPathLength->SetParameter(0, 0.); //Impact Parameter
313 fgWPathLength->SetParameter(1, 1000.); //Number of interactions used for average
314 fgWPathLength->SetParameter(2, 0); //Pathlength definition
317 void AliFastGlauber::Reset() const
320 // Reset dynamic allocated formulas
321 // in case init is called twice
323 if(fgWSb) delete fgWSb;
324 if(fgRWSb) delete fgRWSb;
325 if(fgWSbz) delete fgWSbz;
326 if(fgWSz) delete fgWSz;
327 if(fgWSta) delete fgWSta;
328 if(fgWStarfi) delete fgWStarfi;
329 if(fgWAlmond) delete fgWAlmond;
330 if(fgWStaa) delete fgWStaa;
331 if(fgWSgeo) delete fgWSgeo;
332 if(fgWSbinary) delete fgWSbinary;
333 if(fgWSN) delete fgWSN;
334 if(fgWPathLength0) delete fgWPathLength0;
335 if(fgWPathLength) delete fgWPathLength;
336 if(fgWEnergyDensity) delete fgWEnergyDensity;
337 if(fgWIntRadius) delete fgWIntRadius;
338 if(fgWKParticipants) delete fgWKParticipants;
339 if(fgWParticipants) delete fgWParticipants;
342 void AliFastGlauber::DrawWSb() const
345 // Draw Wood-Saxon Nuclear Density Function
347 TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
349 Double_t max=fgWSb->GetMaximum(0.,fgBMax)*1.01;
350 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);
352 h2f->GetXaxis()->SetTitle("r [fm]");
353 h2f->GetYaxis()->SetNoExponent(kTRUE);
354 h2f->GetYaxis()->SetTitle("#rho [fm^{-3}]");
357 TLegend *l1a = new TLegend(0.45,0.6,.90,0.8);
358 l1a->SetFillStyle(0);
359 l1a->SetBorderSize(0);
361 snprintf(label,100, "r_{0} = %.2f fm",fWSr0);
362 l1a->AddEntry(fgWSb,label,"");
363 snprintf(label,100, "d = %.2f fm",fWSd);
364 l1a->AddEntry(fgWSb,label,"");
365 snprintf(label,100, "n = %.2e fm^{-3}",fWSn);
366 l1a->AddEntry(fgWSb,label,"");
367 snprintf(label,100, "#omega = %.2f",fWSw);
368 l1a->AddEntry(fgWSb,label,"");
373 void AliFastGlauber::DrawOverlap() const
376 // Draw Overlap Function
378 TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
380 Double_t max=fgWStaa->GetMaximum(0,fgBMax)*1.01;
381 TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0, max);
383 h2f->GetXaxis()->SetTitle("b [fm]");
384 h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]");
386 fgWStaa->Draw("same");
389 void AliFastGlauber::DrawParticipants() const
392 // Draw Number of Participants Npart
394 TCanvas *c3 = new TCanvas("c3","Participants",400,10,600,700);
396 Double_t max=fgWParticipants->GetMaximum(0,fgBMax)*1.01;
397 TH2F *h2f=new TH2F("h2fpart","Number of Participants",2,0,fgBMax,2,0,max);
399 h2f->GetXaxis()->SetTitle("b [fm]");
400 h2f->GetYaxis()->SetTitle("N_{part}");
402 fgWParticipants->Draw("same");
403 TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
404 l1a->SetFillStyle(0);
405 l1a->SetBorderSize(0);
407 snprintf(label,100, "#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
408 l1a->AddEntry(fgWParticipants,label,"");
413 void AliFastGlauber::DrawThickness() const
416 // Draw Thickness Function
418 TCanvas *c4 = new TCanvas("c4","Thickness",400,10,600,700);
420 Double_t max=fgWSta->GetMaximum(0,fgBMax)*1.01;
421 TH2F *h2f=new TH2F("h2fta","Thickness function: T_{A} [fm^{-2}]",2,0,fgBMax,2,0,max);
423 h2f->GetXaxis()->SetTitle("b [fm]");
424 h2f->GetYaxis()->SetTitle("T_{A} [fm^{-2}]");
426 fgWSta->Draw("same");
429 void AliFastGlauber::DrawGeo() const
432 // Draw Geometrical Cross-Section
434 TCanvas *c5 = new TCanvas("c5","Geometrical Cross-Section",400,10,600,700);
436 Double_t max=fgWSgeo->GetMaximum(0,fgBMax)*1.01;
437 TH2F *h2f=new TH2F("h2fgeo","Differential Geometrical Cross-Section: d#sigma^{geo}_{AB}/db [fm]",2,0,fgBMax,2,0,max);
439 h2f->GetXaxis()->SetTitle("b [fm]");
440 h2f->GetYaxis()->SetTitle("d#sigma^{geo}_{AB}/db [fm]");
442 fgWSgeo->Draw("same");
443 TLegend *l1a = new TLegend(0.10,0.8,.40,0.9);
444 l1a->SetFillStyle(0);
445 l1a->SetBorderSize(0);
447 snprintf(label,100, "#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
448 l1a->AddEntry(fgWSgeo,label,"");
453 void AliFastGlauber::DrawBinary() const
456 // Draw Binary Cross-Section
458 TCanvas *c6 = new TCanvas("c6","Binary Cross-Section",400,10,600,700);
460 Double_t max=fgWSbinary->GetMaximum(0,fgBMax)*1.01;
461 TH2F *h2f=new TH2F("h2fbinary","Differential Binary Cross-Section: #sigma^{hard}_{NN} dT_{AB}/db [fm]",2,0,fgBMax,2,0,max);
463 h2f->GetXaxis()->SetTitle("b [fm]");
464 h2f->GetYaxis()->SetTitle("d#sigma^{hard}_{AB}/db [fm]");
466 fgWSbinary->Draw("same");
467 TLegend *l1a = new TLegend(0.50,0.8,.90,0.9);
468 l1a->SetFillStyle(0);
469 l1a->SetBorderSize(0);
471 snprintf(label,100, "#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
472 l1a->AddEntry(fgWSb,label,"");
477 void AliFastGlauber::DrawN() const
480 // Draw Binaries per event (Ncoll)
482 TCanvas *c7 = new TCanvas("c7","Binaries per event",400,10,600,700);
484 Double_t max=fgWSN->GetMaximum(0.01,fgBMax)*1.01;
485 TH2F *h2f=new TH2F("h2fhardcols","Number of hard collisions: T_{AB} #sigma^{hard}_{NN}/#sigma_{AB}^{geo}",2,0,fgBMax,2,0,max);
487 h2f->GetXaxis()->SetTitle("b [fm]");
488 h2f->GetYaxis()->SetTitle("N_{coll}");
491 TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
492 l1a->SetFillStyle(0);
493 l1a->SetBorderSize(0);
495 snprintf(label,100, "#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
496 l1a->AddEntry(fgWSN,label,"");
497 snprintf(label,100, "#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
498 l1a->AddEntry(fgWSN,label,"");
503 void AliFastGlauber::DrawKernel(Double_t b) const
508 TCanvas *c8 = new TCanvas("c8","Kernel",400,10,600,700);
510 fgWStarfi->SetParameter(0, b);
511 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());
513 h2f->GetXaxis()->SetTitle("r [fm]");
514 h2f->GetYaxis()->SetTitle("#phi [rad]");
516 fgWStarfi->Draw("same");
517 TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
518 l1a->SetFillStyle(0);
519 l1a->SetBorderSize(0);
521 snprintf(label, 100, "b = %.1f fm",b);
522 l1a->AddEntry(fgWStarfi,label,"");
527 void AliFastGlauber::DrawAlmond(Double_t b) const
530 // Draw Interaction Almond
532 TCanvas *c9 = new TCanvas("c9","Almond",400,10,600,700);
534 fgWAlmond->SetParameter(0, b);
535 TH2F *h2f=new TH2F("h2falmond","Interaction Almond [fm^{-4}]",2,-fgBMax, fgBMax, 2, -fgBMax, fgBMax);
537 h2f->GetXaxis()->SetTitle("x [fm]");
538 h2f->GetYaxis()->SetTitle("y [fm]");
540 gStyle->SetPalette(1);
541 fgWAlmond->Draw("colzsame");
542 TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
543 l1a->SetFillStyle(0);
544 l1a->SetBorderSize(0);
546 snprintf(label, 100, "b = %.1f fm",b);
547 l1a->AddEntry(fgWAlmond,label,"");
552 void AliFastGlauber::DrawEnergyDensity() const
555 // Draw energy density
557 TCanvas *c10 = new TCanvas("c10","Energy Density",400, 10, 600, 700);
559 fgWEnergyDensity->SetMinimum(0.);
560 Double_t max=fgWEnergyDensity->GetMaximum(0,fgWEnergyDensity->GetParameter(0))*1.01;
561 TH2F *h2f=new TH2F("h2fenergydens","Energy density",2,0,fgBMax,2,0,max);
563 h2f->GetXaxis()->SetTitle("b [fm]");
564 h2f->GetYaxis()->SetTitle("fm^{-4}");
566 fgWEnergyDensity->Draw("same");
570 void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt) const
575 TCanvas *c11 = new TCanvas("c11","Path Length",400,10,600,700);
577 fgWPathLength0->SetParameter(0, b);
578 fgWPathLength0->SetParameter(1, Double_t(iopt));
579 fgWPathLength0->SetMinimum(0.);
580 fgWPathLength0->SetMaximum(10.);
581 TH2F *h2f=new TH2F("h2fpathlength0","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
583 h2f->GetXaxis()->SetTitle("#phi [rad]");
584 h2f->GetYaxis()->SetTitle("l [fm]");
586 fgWPathLength0->Draw("same");
589 void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt) const
594 TCanvas *c12 = new TCanvas("c12","Path Length",400,10,600,700);
596 fgWAlmond->SetParameter(0, b);
597 fgWPathLength->SetParameter(0, b);
598 fgWPathLength->SetParameter(1, Double_t (ni));
599 fgWPathLength->SetParameter(2, Double_t (iopt));
600 fgWPathLength->SetMinimum(0.);
601 fgWPathLength->SetMaximum(10.);
602 TH2F *h2f=new TH2F("h2fpathlength","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
604 h2f->GetXaxis()->SetTitle("#phi [rad]");
605 h2f->GetYaxis()->SetTitle("l [fm]");
607 fgWPathLength->Draw("same");
610 void AliFastGlauber::DrawIntRadius(Double_t b) const
613 // Draw Interaction Radius
615 TCanvas *c13 = new TCanvas("c13","Interaction Radius",400,10,600,700);
617 fgWIntRadius->SetParameter(0, b);
618 fgWIntRadius->SetMinimum(0);
619 Double_t max=fgWIntRadius->GetMaximum(0,fgBMax)*1.01;
620 TH2F *h2f=new TH2F("h2fintradius","Interaction Density",2,0.,fgBMax,2,0,max);
622 h2f->GetXaxis()->SetTitle("r [fm]");
623 h2f->GetYaxis()->SetTitle("[fm^{-3}]");
625 fgWIntRadius->Draw("same");
628 Double_t AliFastGlauber::WSb(const Double_t* x, const Double_t* par)
631 // Woods-Saxon Parameterisation
632 // as a function of radius (xx)
634 const Double_t kxx = x[0]; //fm
635 const Double_t kr0 = par[0]; //fm
636 const Double_t kd = par[1]; //fm
637 const Double_t kw = par[2]; //no units
638 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
639 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
643 Double_t AliFastGlauber::RWSb(const Double_t* x, const Double_t* par)
646 // Woods-Saxon Parameterisation
647 // as a function of radius (xx)
649 const Double_t kxx = x[0]; //fm
650 const Double_t kr0 = par[0]; //fm
651 const Double_t kd = par[1]; //fm
652 const Double_t kw = par[2]; //no units
653 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
654 Double_t y = kxx * kxx * kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
659 Double_t AliFastGlauber::WSbz(const Double_t* x, const Double_t* par)
662 // Wood Saxon Parameterisation
663 // as a function of z and b
665 const Double_t kbb = x[0]; //fm
666 const Double_t kzz = x[1]; //fm
667 const Double_t kr0 = par[0]; //fm
668 const Double_t kd = par[1]; //fm
669 const Double_t kw = par[2]; //no units
670 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
671 const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
672 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
676 Double_t AliFastGlauber::WSz(const Double_t* x, const Double_t* par)
679 // Wood Saxon Parameterisation
680 // as a function of z for fixed b
682 const Double_t kzz = x[0]; //fm
683 const Double_t kr0 = par[0]; //fm
684 const Double_t kd = par[1]; //fm
685 const Double_t kw = par[2]; //no units
686 const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
687 const Double_t kbb = par[4]; //fm
688 const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
689 Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
693 Double_t AliFastGlauber::WSta(const Double_t* x, const Double_t* /*par*/)
696 // Thickness function T_A
697 // as a function of b
699 const Double_t kb = x[0];
700 fgWSz->SetParameter(4, kb);
701 Double_t y = 2. * fgWSz->Integral(0., fgBMax);
705 Double_t AliFastGlauber::WStarfi(const Double_t* x, const Double_t* par)
708 // Kernel for overlap function: T_A(s)*T_A(s-b)
709 // as a function of r and phi
710 const Double_t kr1 = x[0];
711 const Double_t kphi = x[1];
712 const Double_t kb = par[0];
713 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
714 Double_t y = kr1 * fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
718 Double_t AliFastGlauber::WStaa(const Double_t* x, const Double_t* par)
722 // T_{AB}=Int d2s T_A(s)*T_B(s-b)
723 // as a function of b
724 // (normalized to fA*fB)
726 const Double_t kb = x[0];
727 const Double_t ka = par[0];
728 fgWStarfi->SetParameter(0, kb);
730 // root integration seems to fail
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);
750 for (Int_t i = 0; i < fgkMCInts; i++)
753 const Double_t kphi = TMath::Pi() * gRandom->Rndm();
754 const Double_t kb1 = fgBMax * gRandom->Rndm();
755 y += fgWStarfi->Eval(kb1, kphi);
757 y *= 2. * TMath::Pi() * fgBMax / fgkMCInts; //fm^-2
758 y *= ka * ka * 0.1; //mbarn^-1
762 Double_t AliFastGlauber::WKParticipants(const Double_t* x, const Double_t* par)
765 // Kernel for number of participants
766 // as a function of r and phi
768 const Double_t kr1 = x[0];
769 const Double_t kphi = x[1];
770 const Double_t kb = par[0]; //fm
771 const Double_t ksig = par[1]; //mbarn
772 const Double_t ka = par[2]; //mass number
773 const Double_t kr2 = TMath::Sqrt(kr1*kr1 +kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
774 const Double_t kxsi = fgWSta->Eval(kr2) * ksig * 0.1; //no units
776 Double_t y=(1-TMath::Power((1-xsi),aa))
779 Double_t sum = ka * kxsi;
781 for (Int_t i = 1; i <= ka; i++)
784 sum *= (-kxsi) * a / Float_t(i+1);
787 y *= kr1 * fgWSta->Eval(kr1);
791 Double_t AliFastGlauber::WParticipants(const Double_t* x, const Double_t* par)
794 // Number of Participants as
797 const Double_t kb = x[0];
798 const Double_t ksig = par[0]; //mbarn
799 const Double_t ka = par[1]; //mass number
800 fgWKParticipants->SetParameter(0, kb);
801 fgWKParticipants->SetParameter(1, ksig);
802 fgWKParticipants->SetParameter(2, ka);
808 for (Int_t i = 0; i < fgkMCInts; i++)
810 const Double_t kphi = TMath::Pi() * gRandom->Rndm();
811 const Double_t kb1 = fgBMax * gRandom->Rndm();
812 y += fgWKParticipants->Eval(kb1, kphi);
814 y *= 2. * ka * 2. * TMath::Pi() * fgBMax / fgkMCInts;
818 Double_t AliFastGlauber::WSgeo(const Double_t* x, const Double_t* par)
821 // Geometrical Cross-Section
822 // as a function of b
824 const Double_t kb = x[0]; //fm
825 const Double_t ksigNN = par[0]; //mbarn
826 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
827 Double_t y = 2. * TMath::Pi() * kb * (1. - TMath::Exp(- ksigNN * ktaa));
831 Double_t AliFastGlauber::WSbinary(const Double_t* x, const Double_t* par)
834 // Number of binary hard collisions
835 // as a function of b
837 const Double_t kb = x[0]; //fm
838 const Double_t ksig = par[0]; //mbarn
839 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
840 Double_t y = 2. * TMath::Pi() * kb * ksig * ktaa;
844 Double_t AliFastGlauber::WSN(const Double_t* x, const Double_t* /*par*/)
847 // Number of hard processes per event
848 // as a function of b
849 const Double_t kb = x[0];
850 Double_t y = fgWSbinary->Eval(kb)/fgWSgeo->Eval(kb);
854 Double_t AliFastGlauber::WEnergyDensity(const Double_t* x, const Double_t* par)
857 // Initial energy density
858 // as a function of the impact parameter
860 const Double_t kb = x[0];
861 const Double_t krA = par[0];
863 // Attention: area of transverse reaction zone in hard-sphere approximation !
864 const Double_t krA2=krA*krA;
865 const Double_t kb2=kb*kb;
866 const Double_t ksaa = (TMath::Pi() - 2. * TMath::ASin(kb/ 2./ krA)) * krA2
867 - kb * TMath::Sqrt(krA2 - kb2/ 4.); //fm^2
868 const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
869 Double_t y=ktaa/ksaa*10;
873 Double_t AliFastGlauber::WAlmond(const Double_t* x, const Double_t* par)
876 // Almond shaped interaction region
877 // as a function of cartesian x,y.
879 const Double_t kb = par[0];
880 const Double_t kxx = x[0] + kb/2.;
881 const Double_t kyy = x[1];
882 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
883 const Double_t kphi = TMath::ATan2(kyy,kxx);
884 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
886 // Interaction probability calculated as product of thicknesses
888 Double_t y = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
892 Double_t AliFastGlauber::WIntRadius(const Double_t* x, const Double_t* par)
895 // Average interaction density over radius
896 // at which interaction takes place
897 // as a function of radius
899 const Double_t kr = x[0];
900 const Double_t kb = par[0];
901 fgWAlmond->SetParameter(0, kb);
902 // Average over phi in small steps
903 const Double_t kdphi = 2. * TMath::Pi() / 100.;
906 for (Int_t i = 0; i < 100; i++) {
907 const Double_t kxx = kr * TMath::Cos(phi);
908 const Double_t kyy = kr * TMath::Sin(phi);
909 y += fgWAlmond->Eval(kxx,kyy);
912 // Result multiplied by Jacobian (2 pi r)
913 y *= 2. * TMath::Pi() * kr / 100.;
917 Double_t AliFastGlauber::WPathLength0(const Double_t* x, const Double_t* par)
920 // Path Length as a function of phi
921 // for interaction point fixed at (0,0)
922 // as a function of phi-direction
924 // Phi direction in Almond
925 const Double_t kphi0 = x[0];
926 const Double_t kb = par[0];
927 // Path Length definition
928 const Int_t kiopt = Int_t(par[1]);
930 // Step along radial direction phi
931 const Int_t kNp = 100; // Steps in r
932 const Double_t kDr = fgBMax/kNp;
936 for (Int_t i = 0; i < kNp; i++) {
938 // Transform into target frame
940 const Double_t kxx = r * TMath::Cos(kphi0) + kb / 2.;
941 const Double_t kyy = r * TMath::Sin(kphi0);
942 const Double_t kphi = TMath::ATan2(kyy, kxx);
943 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
944 // Radius in projectile frame
945 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
946 const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
954 if (!kiopt) { // My length definition (is exact for hard disk)
955 if(w) y= 2. * rw / w;
957 const Double_t knorm=fgWSta->Eval(1e-4);
958 if(knorm) y = TMath::Sqrt(2. * rw * kDr / knorm / knorm);
963 Double_t AliFastGlauber::WPathLength(const Double_t* x, const Double_t* par)
966 // Path Length as a function of phi
967 // Interaction point from random distribution
968 // as a function of the phi-direction
969 const Double_t kphi0 = x[0];
970 const Double_t kb = par[0];
971 fgWAlmond->SetParameter(0, kb);
972 const Int_t kNpi = Int_t (par[1]); //Number of interactions
973 const Int_t kiopt = Int_t(par[2]); //Path Length definition
978 const Int_t kNp = 100;
979 const Double_t kDr = fgBMax/Double_t(kNp);
980 Double_t l = 0.; // Path length
981 for (Int_t in = 0; in < kNpi; in ++) {
986 fgWAlmond->GetRandom2(x0, y0);
988 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
989 const Int_t knps = Int_t ((fgBMax - kr0)/kDr) - 1;
993 for (Int_t i = 0; (i < knps ); i++) {
994 // Transform into target frame
995 const Double_t kxx = x0 + r * TMath::Cos(kphi0) + kb / 2.;
996 const Double_t kyy = y0 + r * TMath::Sin(kphi0);
997 const Double_t kphi = TMath::ATan2(kyy, kxx);
998 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
999 // Radius in projectile frame
1000 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
1001 const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1007 // Average over interactions
1009 if(w) l += (2. * rw / w);
1011 const Double_t knorm=fgWSta->Eval(1e-4);
1012 if(knorm) l+= 2. * rw * kDr / knorm / knorm;
1019 ret=TMath::Sqrt( l / kNpi);
1023 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2) const
1026 // Return the geometrical cross-section integrated from b1 to b2
1028 return fgWSgeo->Integral(b1, b2)*10.; //mbarn
1031 Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const
1034 // Return the hard cross-section integrated from b1 to b2
1036 return fgWSbinary->Integral(b1, b2)*10.; //mbarn
1039 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) const
1042 // Return fraction of hard cross-section integrated from b1 to b2
1044 return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
1047 Double_t AliFastGlauber::NHard(const Double_t b1, const Double_t b2) const
1050 // Number of binary hard collisions
1051 // as a function of b (nucl/ex/0302016 eq. 19)
1053 const Double_t kshard=HardCrossSection(b1,b2);
1054 const Double_t ksgeo=CrossSection(b1,b2);
1056 return kshard/ksgeo;
1060 Double_t AliFastGlauber::Binaries(Double_t b) const
1063 // Return number of binary hard collisions normalized to 1 at b=0
1065 if(b < 1.e-4) b = 1e-4;
1066 return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
1069 Double_t AliFastGlauber::MeanOverlap(Double_t b1, Double_t b2)
1072 // Calculate the mean overlap for impact parameter range b1 .. b2
1078 while (b < b2-0.005) {
1079 Double_t nc = GetNumberOfCollisions(b);
1080 sum += 10. * fgWStaa->Eval(b) * fgWSgeo->Eval(b) * 0.01 / (1. - TMath::Exp(-nc));
1081 sumc += 10. * fgWSgeo->Eval(b) * 0.01;
1084 return (sum / CrossSection(b1, b2));
1088 Double_t AliFastGlauber::MeanNumberOfCollisionsPerEvent(Double_t b1, Double_t b2)
1091 // Calculate the mean number of collisions per event for impact parameter range b1 .. b2
1097 while (b < b2-0.005) {
1098 Double_t nc = GetNumberOfCollisions(b);
1099 sum += nc / (1. - TMath::Exp(-nc)) * 10. * fgWSgeo->Eval(b) * 0.01;
1100 sumc += 10. * fgWSgeo->Eval(b) * 0.01;
1103 return (sum / CrossSection(b1, b2));
1107 Double_t AliFastGlauber::GetNumberOfBinaries(Double_t b) const
1110 // Return number of binary hard collisions at b
1113 return fgWSN->Eval(b);
1116 Double_t AliFastGlauber::Participants(Double_t b) const
1119 // Return the number of participants normalized to 1 at b=0
1122 return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4));
1125 Double_t AliFastGlauber::GetNumberOfParticipants(Double_t b) const
1128 // Return the number of participants for impact parameter b
1131 return (fgWParticipants->Eval(b));
1134 Double_t AliFastGlauber::GetNumberOfCollisions(Double_t b) const
1137 // Return the number of collisions for impact parameter b
1140 return (fgWStaa->Eval(b)*fSigmaNN);
1143 Double_t AliFastGlauber::GetNumberOfCollisionsPerEvent(Double_t b) const
1146 // Return the number of collisions per event (at least one collision)
1147 // for impact parameter b
1149 Double_t n = GetNumberOfCollisions(b);
1151 return (n / (1. - TMath::Exp(- n)));
1157 void AliFastGlauber::SimulateTrigger(Int_t n)
1160 // Simulates Trigger
1162 TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
1163 TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
1164 TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
1165 TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
1167 mbtH->SetXTitle("b [fm]");
1168 hdtH->SetXTitle("b [fm]");
1169 mbmH->SetXTitle("Multiplicity");
1170 hdmH->SetXTitle("Multiplicity");
1172 TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
1174 TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
1180 for (Int_t iev = 0; iev < n; iev++)
1183 GetRandom(b, p, mult);
1186 mbmH->Fill(mult, 1.);
1187 hdmH->Fill(mult, p);
1203 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
1206 // Gives back a random impact parameter, hard trigger probability and multiplicity
1208 b = fgWSgeo->GetRandom();
1209 const Float_t kmu = fgWSN->Eval(b);
1210 p = 1.-TMath::Exp(-kmu);
1211 mult = 6000./fgWSN->Eval(1.) * kmu;
1214 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
1217 // Gives back a random impact parameter bin, and hard trigger decission
1219 const Float_t kb = fgWSgeo->GetRandom();
1220 const Float_t kmu = fgWSN->Eval(kb) * fSigmaHard;
1221 const Float_t kp = 1.-TMath::Exp(-kmu);
1224 } else if (kb < 8.6) {
1226 } else if (kb < 11.2) {
1228 } else if (kb < 13.2) {
1230 } else if (kb < 15.0) {
1236 const Float_t kr = gRandom->Rndm();
1237 if (kr < kp) hard = kTRUE;
1240 Double_t AliFastGlauber::GetRandomImpactParameter(Double_t bmin, Double_t bmax)
1243 // Gives back a random impact parameter in the range bmin .. bmax
1246 while(b < bmin || b > bmax)
1247 b = fgWSgeo->GetRandom();
1251 void AliFastGlauber::StoreFunctions() const
1254 // Store in file functions
1256 TFile* ff = new TFile(fName.Data(),"recreate");
1257 fgWStaa->Write("WStaa");
1258 fgWParticipants->Write("WParticipants");
1263 //=================== Added by A. Dainese 11/02/04 ===========================
1265 void AliFastGlauber::StoreAlmonds() const
1269 // 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
1271 Char_t almondName[100];
1272 TFile* ff = new TFile(fName.Data(),"update");
1273 for(Int_t k=0; k<40; k++) {
1274 snprintf(almondName, 100, "WAlmondFixedB%d",k);
1275 Double_t b = 0.25+k*0.5;
1276 Info("StoreAlmonds"," b = %f\n",b);
1277 fgWAlmond->SetParameter(0,b);
1278 fgWAlmond->Write(almondName);
1284 void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
1287 // Set limits of centrality class as fractions
1288 // of the geomtrical cross section
1290 if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
1291 Error("SetCentralityClass", "Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
1295 Double_t bLow=0.,bUp=0.;
1297 const Double_t knorm=fgWSgeo->Integral(0.,100.);
1298 while(xsecFr<xsecFrLow) {
1299 xsecFr = fgWSgeo->Integral(0.,bLow)/knorm;
1303 while(xsecFr<xsecFrUp) {
1304 xsecFr = fgWSgeo->Integral(0.,bUp)/knorm;
1308 Info("SetCentralityClass", "Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm",
1309 xsecFrLow,xsecFrUp,bLow,bUp);
1310 fgWSbinary->SetRange(bLow,bUp);
1316 void AliFastGlauber::GetRandomBHard(Double_t& b)
1319 // Get random impact parameter according to distribution of
1320 // hard (binary) cross-section, in the range defined by the centrality class
1322 b = fgWSbinary->GetRandom();
1323 Int_t bin = 2*(Int_t)b;
1324 if( (b-(Int_t)b) > 0.5) bin++;
1325 fgWAlmondCurrent = fgWAlmondFixedB[bin];
1329 void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
1332 // Get random position of parton production point according to
1333 // product of thickness functions
1335 fgWAlmondCurrent->GetRandom2(x,y);
1339 void AliFastGlauber::GetRandomPhi(Double_t& phi)
1342 // Get random parton azimuthal propagation direction
1344 phi = 2.*TMath::Pi()*gRandom->Rndm();
1348 Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Double_t phi0)
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
1356 // number of steps in l
1357 const Int_t kNp = 100;
1358 const Double_t kDl = fgBMax/Double_t(kNp);
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)
1369 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1370 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1372 Double_t integral1 = 0.;
1373 Double_t integral2 = 0.;
1375 for (Int_t i = 0; i < knps; i++) {
1377 // Transform into target frame
1378 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1379 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1380 const Double_t kphi = TMath::ATan2(kyy, kxx);
1381 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
1382 // Radius in projectile frame
1383 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1384 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1386 integral1 += kprodTATB * l * kDl;
1387 integral2 += kprodTATB * kDl;
1393 ell = (2. * integral1 / integral2);
1395 } else if(fEllDef==2) {
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))
1404 const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1405 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1406 const Double_t kprodTATBHalfMax = 0.5*fgWAlmondCurrent->Eval(0.,0.);
1409 Double_t integral = 0.;
1410 for (Int_t i = 0; i < knps; i++) {
1411 // Transform into target frame
1412 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1413 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1414 const Double_t kphi = TMath::ATan2(kyy, kxx);
1415 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
1416 // Radius in projectile frame
1417 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1418 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1419 if(kprodTATB>kprodTATBHalfMax) integral += kDl;
1422 Double_t ell = integral;
1425 Error("CalculateLength","Wrong length definition setting: %d !\n",fEllDef);
1430 void AliFastGlauber::GetLengthAndPhi(Double_t& ell,Double_t& phi,Double_t b)
1433 // Return length from random b, x0, y0, phi0
1436 Double_t x0,y0,phi0;
1437 if(b<0.) GetRandomBHard(b);
1441 ell = CalculateLength(b,x0,y0,phi0);
1445 void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
1448 // Return length from random b, x0, y0, phi0
1451 GetLengthAndPhi(ell,phi,b);
1455 void AliFastGlauber::GetLengthsBackToBackAndPhi(Double_t& ell1,Double_t& ell2,Double_t &phi,Double_t b)
1458 // Return 2 lengths back to back from random b, x0, y0, phi0
1461 Double_t x0,y0,phi0;
1462 if(b<0.) GetRandomBHard(b);
1465 const Double_t kphi0plusPi = phi0+TMath::Pi();
1467 ell1 = CalculateLength(b,x0,y0,phi0);
1468 ell2 = CalculateLength(b,x0,y0,kphi0plusPi);
1472 void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
1476 // Return 2 lengths back to back from random b, x0, y0, phi0
1479 GetLengthsBackToBackAndPhi(ell1,ell2,phi,b);
1483 void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* const phi,Double_t* ell, Double_t b)
1486 // Returns lenghts for n partons with azimuthal angles phi[n]
1487 // from random b, x0, y0
1490 if(b < 0.) GetRandomBHard(b);
1492 for(Int_t i = 0; i< n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
1496 void AliFastGlauber::PlotBDistr(Int_t n)
1499 // Plot distribution of n impact parameters
1502 TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax);
1503 hB->SetXTitle("b [fm]");
1504 hB->SetYTitle("dN/db [a.u.]");
1505 hB->SetFillColor(3);
1506 for(Int_t i=0; i<n; i++) {
1510 TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
1516 void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,const char *fname)
1519 // Plot length distribution
1522 TH1F *hEll = new TH1F("hEll","Length distribution",64,-0.5,15);
1523 hEll->SetXTitle("Transverse path length, L [fm]");
1524 hEll->SetYTitle("Probability");
1525 hEll->SetFillColor(2);
1526 for(Int_t i=0; i<n; i++) {
1530 hEll->Scale(1/(Double_t)n);
1531 TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
1536 TFile *f = new TFile(fname,"recreate");
1543 void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,const char *fname)
1546 // Plot lengths back-to-back distributions
1549 TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15);
1550 hElls->SetXTitle("Transverse path length, L [fm]");
1551 hElls->SetYTitle("Transverse path length, L [fm]");
1552 for(Int_t i=0; i<n; i++) {
1553 GetLengthsBackToBack(ell1,ell2);
1554 hElls->Fill(ell1,ell2);
1556 hElls->Scale(1/(Double_t)n);
1557 TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
1558 gStyle->SetPalette(1,0);
1560 hElls->Draw("col,Z");
1562 TFile *f = new TFile(fname,"recreate");
1569 void AliFastGlauber::PlotAlmonds() const
1572 // Plot almonds for some impact parameters
1574 TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
1575 gStyle->SetPalette(1,0);
1578 fgWAlmondFixedB[0]->Draw("cont1");
1580 fgWAlmondFixedB[10]->Draw("cont1");
1582 fgWAlmondFixedB[20]->Draw("cont1");
1584 fgWAlmondFixedB[30]->Draw("cont1");
1588 //=================== Added by A. Dainese 05/03/04 ===========================
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
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)
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
1604 // number of steps in l
1605 const Int_t kNp = 100;
1606 const Double_t kDl = fgBMax/Double_t(kNp);
1609 const Double_t kr0 = TMath::Sqrt(x0 * x0 + y0 * y0);
1610 const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1617 while((i < knps) && (l < ellCut)) {
1618 // Transform into target frame
1619 const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
1620 const Double_t kyy = y0 + l * TMath::Sin(phi0);
1621 const Double_t kphi = TMath::ATan2(kyy, kxx);
1622 const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
1623 // Radius in projectile frame
1624 const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1625 const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1626 integral0 += kprodTATB * kDl;
1627 integral1 += kprodTATB * l * kDl;
1634 void AliFastGlauber::GetI0I1AndPhi(Double_t& integral0,Double_t& integral1,
1636 Double_t ellCut,Double_t b)
1639 // Return I0 and I1 from random b, x0, y0, phi0
1642 Double_t x0,y0,phi0;
1643 if(b<0.) GetRandomBHard(b);
1647 CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut);
1651 void AliFastGlauber::GetI0I1(Double_t& integral0,Double_t& integral1,
1652 Double_t ellCut,Double_t b)
1655 // Return I0 and I1 from random b, x0, y0, phi0
1658 GetI0I1AndPhi(integral0,integral1,phi,ellCut,b);
1662 void AliFastGlauber::GetI0I1BackToBackAndPhi(Double_t& integral01,Double_t& integral11,
1663 Double_t& integral02,Double_t& integral12,
1665 Double_t ellCut,Double_t b)
1668 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1671 Double_t x0,y0,phi0;
1672 if(b<0.) GetRandomBHard(b);
1676 const Double_t kphi0plusPi = phi0+TMath::Pi();
1677 CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
1678 CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
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)
1688 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1691 Double_t x0,y0,phi0;
1692 if(b<0.) GetRandomBHard(b);
1695 phi = phi0; x=x0; y=y0;
1696 const Double_t kphi0plusPi = phi0+TMath::Pi();
1697 CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
1698 CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
1702 void AliFastGlauber::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11,
1703 Double_t& integral02,Double_t& integral12,
1704 Double_t ellCut,Double_t b)
1707 // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1710 GetI0I1BackToBackAndPhi(integral01,integral11,integral02,integral12,
1715 void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
1716 Double_t* integral0,Double_t* integral1,
1717 Double_t ellCut,Double_t b)
1720 // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
1721 // from random b, x0, y0
1724 if(b<0.) GetRandomBHard(b);
1726 for(Int_t i=0; i<n; i++)
1727 CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
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)
1737 // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
1738 // from random b, x0, y0 and return x0,y0
1741 if(b<0.) GetRandomBHard(b);
1743 for(Int_t i=0; i<n; i++)
1744 CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1750 void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
1751 Bool_t save,const char *fname)
1754 // Plot I0-I1 distribution
1757 TH2F *hI0I1s = new TH2F("hI0I1s","I_{0} versus I_{1}",1000,0,0.001,1000,0,0.01);
1758 hI0I1s->SetXTitle("I_{0} [fm^{-3}]");
1759 hI0I1s->SetYTitle("I_{1} [fm^{-2}]");
1761 TH1F *hI0 = new TH1F("hI0","I_{0} = #hat{q}L / k",
1763 hI0->SetXTitle("I_{0} [fm^{-3}]");
1764 hI0->SetYTitle("Probability");
1765 hI0->SetFillColor(3);
1766 TH1F *hI1 = new TH1F("hI1","I_{1} = #omega_{c} / k",
1768 hI1->SetXTitle("I_{1} [fm^{-2}]");
1769 hI1->SetYTitle("Probability");
1770 hI1->SetFillColor(4);
1771 TH1F *h2 = new TH1F("h2","2 I_{1}^{2}/I_{0} = R / k",
1773 h2->SetXTitle("2 I_{1}^{2}/I_{0} [fm^{-1}]");
1774 h2->SetYTitle("Probability");
1775 h2->SetFillColor(2);
1776 TH1F *h3 = new TH1F("h3","2 I_{1}/I_{0} = L",
1778 h3->SetXTitle("2 I_{1}/I_{0} [fm]");
1779 h3->SetYTitle("Probability");
1780 h3->SetFillColor(5);
1781 TH1F *h4 = new TH1F("h4","I_{0}^{2}/(2 I_{1}) = #hat{q} / k",
1783 h4->SetXTitle("I_{0}^{2}/(2 I_{1}) [fm^{-4}]");
1784 h4->SetYTitle("Probability");
1785 h4->SetFillColor(7);
1787 for(Int_t i=0; i<n; i++) {
1788 GetI0I1(i0,i1,ellCut);
1789 hI0I1s->Fill(i0,i1);
1792 h2->Fill(2.*i1*i1/i0);
1794 h4->Fill(i0*i0/2./i1);
1796 hI0->Scale(1/(Double_t)n);
1797 hI1->Scale(1/(Double_t)n);
1798 h2->Scale(1/(Double_t)n);
1799 h3->Scale(1/(Double_t)n);
1800 h4->Scale(1/(Double_t)n);
1801 hI0I1s->Scale(1/(Double_t)n);
1803 TCanvas *cI0I1 = new TCanvas("cI0I1","I0 and I1",0,0,900,700);
1816 gStyle->SetPalette(1,0);
1817 hI0I1s->Draw("col,Z");
1820 TFile *f = new TFile(fname,"recreate");
1832 void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
1833 Bool_t save,const char *fname)
1836 // Plot I0-I1 back-to-back distributions
1838 Double_t i01,i11,i02,i12;
1839 TH2F *hI0s = new TH2F("hI0s","I_{0}'s back-to-back",100,0,100,100,0,100);
1840 hI0s->SetXTitle("I_{0} [fm^{-3}]");
1841 hI0s->SetYTitle("I_{0} [fm^{-3}]");
1842 TH2F *hI1s = new TH2F("hI1s","I_{1}'s back-to-back",100,0,100,100,0,100);
1843 hI1s->SetXTitle("I_{1} [fm^{-2}]");
1844 hI1s->SetYTitle("I_{1} [fm^{-2}]");
1846 for(Int_t i=0; i<n; i++) {
1847 GetI0I1BackToBack(i01,i11,i02,i12,ellCut);
1848 hI0s->Fill(i01,i02);
1849 hI1s->Fill(i11,i12);
1851 hI0s->Scale(1/(Double_t)n);
1852 hI1s->Scale(1/(Double_t)n);
1854 TCanvas *cI0I1s = new TCanvas("cI0I1s","I0 and I1 back-to-back distributions",0,0,800,400);
1855 gStyle->SetPalette(1,0);
1856 cI0I1s->Divide(2,1);
1858 hI0s->Draw("col,Z");
1860 hI1s->Draw("col,Z");
1863 TFile *f = new TFile(fname,"recreate");
1871 AliFastGlauber& AliFastGlauber::operator=(const AliFastGlauber& rhs)
1873 // Assignment operator
1878 void AliFastGlauber::Copy(TObject&) const
1883 Fatal("Copy","Not implemented!\n");