**************************************************************************/
/* $Id$ */
+//
+// Utility class to make simple Glauber type calculations for collision geometries:
+// Impact parameter, production points, reaction plane dependence
+// The SimulateTrigger method can be used for simple MB and hard-process
+// (binary scaling) trigger studies.
+// Some basic quantities can be visualized directly.
+// The default set-up for PbPb collisions can be read from a file calling Init(1).
+//
+//
+// Author: andreas.morsch@cern.ch
// from AliRoot
#include "AliFastGlauber.h"
ClassImp(AliFastGlauber)
-TF1* AliFastGlauber::fWSb = NULL;
-TF2* AliFastGlauber::fWSbz = NULL;
-TF1* AliFastGlauber::fWSz = NULL;
-TF1* AliFastGlauber::fWSta = NULL;
-TF2* AliFastGlauber::fWStarfi = NULL;
-TF2* AliFastGlauber::fWAlmond = NULL;
-TF1* AliFastGlauber::fWStaa = NULL;
-TF1* AliFastGlauber::fWSgeo = NULL;
-TF1* AliFastGlauber::fWSbinary = NULL;
-TF1* AliFastGlauber::fWSN = NULL;
-TF1* AliFastGlauber::fWPathLength0 = NULL;
-TF1* AliFastGlauber::fWPathLength = NULL;
-TF1* AliFastGlauber::fWEnergyDensity = NULL;
-TF1* AliFastGlauber::fWIntRadius = NULL;
-Float_t AliFastGlauber::fbMax = 0.;
+TF1* AliFastGlauber::fgWSb = NULL;
+TF2* AliFastGlauber::fgWSbz = NULL;
+TF1* AliFastGlauber::fgWSz = NULL;
+TF1* AliFastGlauber::fgWSta = NULL;
+TF2* AliFastGlauber::fgWStarfi = NULL;
+TF2* AliFastGlauber::fgWAlmond = NULL;
+TF1* AliFastGlauber::fgWStaa = NULL;
+TF1* AliFastGlauber::fgWSgeo = NULL;
+TF1* AliFastGlauber::fgWSbinary = NULL;
+TF1* AliFastGlauber::fgWSN = NULL;
+TF1* AliFastGlauber::fgWPathLength0 = NULL;
+TF1* AliFastGlauber::fgWPathLength = NULL;
+TF1* AliFastGlauber::fgWEnergyDensity = NULL;
+TF1* AliFastGlauber::fgWIntRadius = NULL;
+Float_t AliFastGlauber::fgBMax = 0.;
AliFastGlauber::AliFastGlauber()
{
//
// Wood-Saxon
//
- fWSb = new TF1("WSb", WSb, 0, fbMax, 4);
- fWSb->SetParameter(0, fWSr0);
- fWSb->SetParameter(1, fWSd);
- fWSb->SetParameter(2, fWSw);
- fWSb->SetParameter(3, fWSn);
+ fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
+ fgWSb->SetParameter(0, fWSr0);
+ fgWSb->SetParameter(1, fWSd);
+ fgWSb->SetParameter(2, fWSw);
+ fgWSb->SetParameter(3, fWSn);
- fWSbz = new TF2("WSbz", WSbz, 0, fbMax, 4);
- fWSbz->SetParameter(0, fWSr0);
- fWSbz->SetParameter(1, fWSd);
- fWSbz->SetParameter(2, fWSw);
- fWSbz->SetParameter(3, fWSn);
+ fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
+ fgWSbz->SetParameter(0, fWSr0);
+ fgWSbz->SetParameter(1, fWSd);
+ fgWSbz->SetParameter(2, fWSw);
+ fgWSbz->SetParameter(3, fWSn);
- fWSz = new TF1("WSz", WSz, 0, fbMax, 5);
- fWSz->SetParameter(0, fWSr0);
- fWSz->SetParameter(1, fWSd);
- fWSz->SetParameter(2, fWSw);
- fWSz->SetParameter(3, fWSn);
+ fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
+ fgWSz->SetParameter(0, fWSr0);
+ fgWSz->SetParameter(1, fWSd);
+ fgWSz->SetParameter(2, fWSw);
+ fgWSz->SetParameter(3, fWSn);
//
// Thickness
//
- fWSta = new TF1("WSta", WSta, 0., fbMax, 0);
+ fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
//
// Overlap Kernel
//
- fWStarfi = new TF2("WStarfi", WStarfi, 0., fbMax, 0., TMath::Pi(), 1);
- fWStarfi->SetParameter(0, 0.);
- fWStarfi->SetNpx(200);
- fWStarfi->SetNpy(20);
+ fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
+ fgWStarfi->SetParameter(0, 0.);
+ fgWStarfi->SetNpx(200);
+ fgWStarfi->SetNpy(20);
//
// Almond shaped interaction region
//
- fWAlmond = new TF2("WAlmond", WAlmond, -fbMax, fbMax, -fbMax, fbMax, 1);
- fWAlmond->SetParameter(0, 0.);
- fWAlmond->SetNpx(200);
- fWAlmond->SetNpy(200);
+ fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
+ fgWAlmond->SetParameter(0, 0.);
+ fgWAlmond->SetNpx(200);
+ fgWAlmond->SetNpy(200);
//
// Path Length as a function of Phi
//
- fWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
- fWPathLength0->SetParameter(0, 0.);
+ fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
+ fgWPathLength0->SetParameter(0, 0.);
// Pathlength definition
- fWPathLength0->SetParameter(1, 0.);
+ fgWPathLength0->SetParameter(1, 0.);
- fWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
+ fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
// Impact Parameter
- fWPathLength->SetParameter(0, 0.);
+ fgWPathLength->SetParameter(0, 0.);
// Number of interactions used for average
- fWPathLength->SetParameter(1, 1000.);
+ fgWPathLength->SetParameter(1, 1000.);
// Pathlength definition
- fWPathLength->SetParameter(2, 0);
+ fgWPathLength->SetParameter(2, 0);
- fWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fbMax, 1);
- fWIntRadius->SetParameter(0, 0.);
+ fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
+ fgWIntRadius->SetParameter(0, 0.);
//
// Overlap
//
if (! mode) {
- fWStaa = new TF1("WStaa", WStaa, 0., fbMax, 0);
- fWStaa->SetNpx(100);
+ fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 0);
+ fgWStaa->SetNpx(100);
} else {
TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
- fWStaa = (TF1*) f->Get("WStaa");
+ fgWStaa = (TF1*) f->Get("WStaa");
}
//
- fWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
- fWEnergyDensity->SetParameter(0, fWSr0 + 1.);
+ fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
+ fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
//
// Geometrical Cross-Section
//
- fWSgeo = new TF1("WSgeo", WSgeo, 0., fbMax, 0);
- fWSgeo->SetNpx(100);
+ fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 0);
+ fgWSgeo->SetNpx(100);
//
// Hard cross section (~ binary collisions)
//
- fWSbinary = new TF1("WSbinary", WSbinary, 0., fbMax, 1);
- fWSbinary->SetParameter(0, fSigmaHard); // mb
- fWSbinary->SetNpx(100);
+ fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
+ fgWSbinary->SetParameter(0, fSigmaHard); // mb
+ fgWSbinary->SetNpx(100);
//
// Hard collisions per event
//
- fWSN = new TF1("WSN", WSN, 0., fbMax, 1);
- fWSN->SetNpx(100);
+ fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
+ fgWSN->SetNpx(100);
}
void AliFastGlauber::DrawWSb()
//
TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
c1->cd();
- fWSb->Draw();
+ fgWSb->Draw();
}
void AliFastGlauber::DrawOverlap()
//
TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
c2->cd();
- fWStaa->Draw();
+ fgWStaa->Draw();
}
void AliFastGlauber::DrawThickness()
//
TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
c3->cd();
- fWSta->Draw();
+ fgWSta->Draw();
}
void AliFastGlauber::DrawGeo()
//
TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
c3->cd();
- fWSgeo->Draw();
+ fgWSgeo->Draw();
}
void AliFastGlauber::DrawBinary()
//
TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
c4->cd();
- fWSbinary->Draw();
+ fgWSbinary->Draw();
}
void AliFastGlauber::DrawN()
//
TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
c5->cd();
- fWSN->Draw();
+ fgWSN->Draw();
}
void AliFastGlauber::DrawKernel(Double_t b)
//
TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
c6->cd();
- fWStarfi->SetParameter(0, b);
- fWStarfi->Draw();
+ fgWStarfi->SetParameter(0, b);
+ fgWStarfi->Draw();
}
void AliFastGlauber::DrawAlmond(Double_t b)
//
TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
c7->cd();
- fWAlmond->SetParameter(0, b);
- fWAlmond->Draw();
+ fgWAlmond->SetParameter(0, b);
+ fgWAlmond->Draw();
}
void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
//
TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
c8->cd();
- fWPathLength0->SetParameter(0, b);
- fWPathLength0->SetParameter(1, Double_t(iopt));
- fWPathLength0->SetMinimum(0.);
- fWPathLength0->SetMaximum(10.);
- fWPathLength0->Draw();
+ fgWPathLength0->SetParameter(0, b);
+ fgWPathLength0->SetParameter(1, Double_t(iopt));
+ fgWPathLength0->SetMinimum(0.);
+ fgWPathLength0->SetMaximum(10.);
+ fgWPathLength0->Draw();
}
void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
//
TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
c9->cd();
- fWAlmond->SetParameter(0, b);
+ fgWAlmond->SetParameter(0, b);
- fWPathLength->SetParameter(0, b);
- fWPathLength->SetParameter(1, Double_t (ni));
- fWPathLength->SetParameter(2, Double_t (iopt));
- fWPathLength->SetMinimum(0.);
- fWPathLength->SetMaximum(10.);
- fWPathLength->Draw();
+ fgWPathLength->SetParameter(0, b);
+ fgWPathLength->SetParameter(1, Double_t (ni));
+ fgWPathLength->SetParameter(2, Double_t (iopt));
+ fgWPathLength->SetMinimum(0.);
+ fgWPathLength->SetMaximum(10.);
+ fgWPathLength->Draw();
}
void AliFastGlauber::DrawIntRadius(Double_t b)
//
TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
c10->cd();
- fWIntRadius->SetParameter(0, b);
- fWIntRadius->SetMinimum(0.);
- fWIntRadius->Draw();
+ fgWIntRadius->SetParameter(0, b);
+ fgWIntRadius->SetMinimum(0.);
+ fgWIntRadius->Draw();
}
void AliFastGlauber::DrawEnergyDensity()
//
TCanvas *c11 = new TCanvas("c11","Energy Density",400, 10, 600, 700);
c11->cd();
- fWEnergyDensity->SetMinimum(0.);
- fWEnergyDensity->Draw();
+ fgWEnergyDensity->SetMinimum(0.);
+ fgWEnergyDensity->Draw();
}
Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
// Thickness function
//
Double_t b = x[0];
- fWSz->SetParameter(4, b);
- Double_t y = 2. * fWSz->Integral(0., fbMax);
+ fgWSz->SetParameter(4, b);
+ Double_t y = 2. * fgWSz->Integral(0., fgBMax);
return y;
}
Double_t r1 = x[0];
Double_t phi = x[1];
Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- Double_t y = r1 * fWSta->Eval(r1) * fWSta->Eval(r2);
+ Double_t y = r1 * fgWSta->Eval(r1) * fgWSta->Eval(r2);
return y;
}
//
// Interaction probability calculated as product of thicknesses
//
- Double_t y = fWSta->Eval(r1) * fWSta->Eval(r2);
+ Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
return y;
}
Double_t r = x[0];
// Impact parameter
Double_t b = par[0];
- fWAlmond->SetParameter(0, b);
+ fgWAlmond->SetParameter(0, b);
// Steps in phi
Double_t dphi = 2. * TMath::Pi() / 100.;
// Average over phi
for (Int_t i = 0; i < 100; i++) {
Double_t xx = r * TMath::Cos(phi);
Double_t yy = r * TMath::Sin(phi);
- y += fWAlmond->Eval(xx,yy);
+ y += fgWAlmond->Eval(xx,yy);
phi += dphi;
} // phi loop
// Result multiplied by Jacobian (2 pi r)
//
//
// Steps in r
- const Int_t np = 100;
- const Double_t dr = fbMax/Double_t(np);
+ const Int_t kNp = 100;
+ const Double_t kDr = fgBMax/Double_t(kNp);
// Impact parameter
Double_t b = par[0];
// Path Length definition
Double_t rw = 0.;
Double_t w = 0.;
// Step along radial direction phi
- for (Int_t i = 0; i < np; i++) {
+ for (Int_t i = 0; i < kNp; i++) {
//
// Transform into target frame
//
Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
// Radius in projectile frame
Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- Double_t y = fWSta->Eval(r1) * fWSta->Eval(r2);
+ Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
rw += y * r;
w += y;
- r += dr;
+ r += kDr;
} // radial steps
//
// My length definition (is exact for hard disk)
if (!iopt) {
return (2. * rw / w);
} else {
- return TMath::Sqrt(2. * rw * dr / fWSta->Eval(0.01) / fWSta->Eval(0.01));
+ return TMath::Sqrt(2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01));
}
}
//
// r-steps
//
- const Int_t np = 100;
- const Double_t dr = fbMax/Double_t(np);
+ const Int_t kNp = 100;
+ const Double_t kDr = fgBMax/Double_t(kNp);
// Number of interactions
- const Int_t npi = Int_t (par[1]);
+ const Int_t kNpi = Int_t (par[1]);
//
// Impact parameter
// Path length
Double_t l = 0.;
- for (Int_t in = 0; in < npi; in ++) {
+ for (Int_t in = 0; in < kNpi; in ++) {
Double_t rw = 0.;
Double_t w = 0.;
// Interaction point
Double_t x0, y0;
- fWAlmond->GetRandom2(x0, y0);
+ fgWAlmond->GetRandom2(x0, y0);
// Initial radius
Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
- Int_t nps = Int_t ((fbMax - r0)/dr) - 1;
+ Int_t nps = Int_t ((fgBMax - r0)/kDr) - 1;
Double_t r = 0.;
// Radial steps
Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
// Radius in projectile frame
Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- Double_t y = fWSta->Eval(r1) * fWSta->Eval(r2);
+ Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
rw += y * r;
w += y;
- r += dr;
+ r += kDr;
} // steps
// Average over interactions
if (!iopt) {
l += (2. * rw / w);
} else {
- l+= 2. * rw * dr / fWSta->Eval(0.01) / fWSta->Eval(0.01);
+ l+= 2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01);
}
} // interactions
if (!iopt)
- return (l / Double_t(npi));
+ return (l / Double_t(kNpi));
else
- return (TMath::Sqrt(l / Double_t(npi)));
+ return (TMath::Sqrt(l / Double_t(kNpi)));
}
Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
// Overlap function
//
Double_t b = x[0];
- fWStarfi->SetParameter(0, b);
+ fgWStarfi->SetParameter(0, b);
/*
Double_t al[2];
Double_t bl[2];
bl[1] = TMath::Pi();
Double_t err;
- Double_t y = 2. * fWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
+ Double_t y = 2. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
printf("WStaa: %f %f %f\n", b, y, err);
*/
//
for (Int_t i = 0; i < 100000; i++)
{
Double_t phi = TMath::Pi() * gRandom->Rndm();
- Double_t b1 = fbMax * gRandom->Rndm();
- y += fWStarfi->Eval(b1, phi);
+ Double_t b1 = fgBMax * gRandom->Rndm();
+ y += fgWStarfi->Eval(b1, phi);
}
- y *= 2. * 0.1 * 208. * 208. * TMath::Pi() * fbMax / 100000.;
+ y *= 2. * 0.1 * 208. * 208. * TMath::Pi() * fgBMax / 100000.;
return y;
}
// Geometrical Cross-Section
//
Double_t b = x[0];
- Double_t taa = fWStaa->Eval(b);
- const Double_t sigma = 55.6; // mbarn
+ Double_t taa = fgWStaa->Eval(b);
+ const Double_t kSigma = 55.6; // mbarn
- Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigma * taa)); // fm
+ Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- kSigma * taa)); // fm
return y;
}
//
Double_t b = x[0];
Double_t sigma = par[0];
- Double_t taa = fWStaa->Eval(b);
+ Double_t taa = fgWStaa->Eval(b);
Double_t y = 2. * TMath::Pi() * b * sigma * taa; // fm
return y;
// Number of hard processes per event
//
Double_t b = x[0];
- Double_t y = fWSbinary->Eval(b)/fWSgeo->Eval(b);
+ Double_t y = fgWSbinary->Eval(b)/fgWSgeo->Eval(b);
return y;
}
// Attention: area of transverse reaction zone in hard-sphere approximation !
Double_t saa = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA * rA
- b * TMath::Sqrt(rA * rA - b * b/ 4.);
- Double_t taa = fWStaa->Eval(b);
+ Double_t taa = fgWStaa->Eval(b);
return (taa/saa);
}
//
// Gives back a random impact parameter, hard trigger probability and multiplicity
//
- b = fWSgeo->GetRandom();
- Float_t mu = fWSN->Eval(b);
+ b = fgWSgeo->GetRandom();
+ Float_t mu = fgWSN->Eval(b);
p = 1.-TMath::Exp(-mu);
- mult = 6000./fWSN->Eval(1.) * mu;
+ mult = 6000./fgWSN->Eval(1.) * mu;
}
void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
//
// Gives back a random impact parameter bin, and hard trigger decission
//
- Float_t b = fWSgeo->GetRandom();
- Float_t mu = fWSN->Eval(b) * fSigmaHard;
+ Float_t b = fgWSgeo->GetRandom();
+ Float_t mu = fgWSN->Eval(b) * fSigmaHard;
Float_t p = 1.-TMath::Exp(-mu);
if (b < 5.) {
bin = 1;
Float_t b = -1.;
while(b < bmin || b > bmax)
- b = fWSgeo->GetRandom();
+ b = fgWSgeo->GetRandom();
return b;
}
// Return cross-section integrated from b1 to b2
//
- return fWSgeo->Integral(b1, b2)/100.;
+ return fgWSgeo->Integral(b1, b2)/100.;
}
Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
// Return raction of hard cross-section integrated from b1 to b2
//
- return fWSbinary->Integral(b1, b2)/fWSbinary->Integral(0., 100.);
+ return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
}
// Return number of binary collisions normalized to 1 at b=0
//
- return fWSN->Eval(b)/fWSN->Eval(0.001);
+ return fgWSN->Eval(b)/fgWSN->Eval(0.001);
}
* See cxx source for full Copyright notice */
/* $Id$ */
+//
+// Utility class to make simple Glauber type calculations for collision geometries:
+// Impact parameter, production points, reaction plane dependence
+//
+// Author: andreas.morsch@cern.ch
#include <TObject.h>
class TF1;
virtual ~AliFastGlauber(){;}
void SetWoodSaxonParameters(Double_t r0, Double_t d, Double_t w, Double_t n)
{fWSr0 = r0; fWSd = d; fWSw = w; fWSn = n;}
- void SetMaxImpact(Float_t bmax = 20.) {fbMax = bmax;};
+ void SetMaxImpact(Float_t bmax = 20.) {fgBMax = bmax;};
void SetHardCrossSection(Float_t xs = 6.6) {fSigmaHard = xs;}
static Double_t WSb (Double_t *xx, Double_t *par);
Double_t CrossSection(Double_t b1, Double_t b2);
Double_t FractionOfHardCrossSection(Double_t b1, Double_t b2);
Double_t Binaries(Double_t b);
- TF2* Kernel() {return fWStarfi;}
- TF1* Overlap() {return fWStaa;}
+ TF2* Kernel() {return fgWStarfi;}
+ TF1* Overlap() {return fgWStaa;}
void SimulateTrigger(Int_t n);
void GetRandom(Float_t& b, Float_t& p, Float_t& mult);
void GetRandom(Int_t& bin, Bool_t& hard);
Float_t GetRandomImpactParameter(Float_t bmin, Float_t bmax);
protected:
- static TF1* fWSb; // Wood-Saxon Function (b)
- static TF2* fWSbz; // Wood-Saxon Function (b, z)
- static TF1* fWSz; // Wood-Saxon Function (b = b0, z)
- static TF1* fWSta; // Thickness Function
- static TF2* fWStarfi; // Kernel for Overlap Function
- static TF1* fWStaa; // Overlap Function
- static TF2* fWAlmond; // Interaction Almond
- static TF1* fWPathLength0; // Path Length as a function of phi
- static TF1* fWPathLength; // Path Length as a function of phi
- static TF1* fWIntRadius; // Interaction Radius
- static TF1* fWSgeo; // dSigma/db geometric
- static TF1* fWSbinary; // dSigma/db binary
- static TF1* fWSN; // dN/db binary
- static TF1* fWEnergyDensity; // Energy density as a function of impact parameter
+ static TF1* fgWSb; // Wood-Saxon Function (b)
+ static TF2* fgWSbz; // Wood-Saxon Function (b, z)
+ static TF1* fgWSz; // Wood-Saxon Function (b = b0, z)
+ static TF1* fgWSta; // Thickness Function
+ static TF2* fgWStarfi; // Kernel for Overlap Function
+ static TF1* fgWStaa; // Overlap Function
+ static TF2* fgWAlmond; // Interaction Almond
+ static TF1* fgWPathLength0; // Path Length as a function of phi
+ static TF1* fgWPathLength; // Path Length as a function of phi
+ static TF1* fgWIntRadius; // Interaction Radius
+ static TF1* fgWSgeo; // dSigma/db geometric
+ static TF1* fgWSbinary; // dSigma/db binary
+ static TF1* fgWSN; // dN/db binary
+ static TF1* fgWEnergyDensity; // Energy density as a function of impact parameter
Float_t fWSr0; // Wood-Saxon Parameter r0
Float_t fWSd; // Wood-Saxon Parameter d
Float_t fWSw; // Wood-Saxon Parameter w
Float_t fWSn; // Wood-Saxon Parameter n
Float_t fSigmaHard; // Hard Cross Section
- static Float_t fbMax; // Maximum Impact Parameter
+ static Float_t fgBMax; // Maximum Impact Parameter
ClassDef(AliFastGlauber,1) // Event geometry simulation in the Glauber Model
};
/*
$Log$
+Revision 1.5 2003/08/13 17:37:29 hristov
+Bug fix (Alpha)
+
Revision 1.4 2003/08/05 16:14:20 morsch
Some problems with too big fluctuations corrected. (A. de Falco)
*/
+// Implementation of AliFastResponse for the Muon Spectrometer resolution.
+// The response depends on the charge of the muon and
+// the background level.
+// The class uses the instance of an object of type AliMUONFastTracking to
+// obtain the smearing parameters.
+// Author: andreas.morsch@cern.ch
+
#include "AliFastMuonTrackingRes.h"
#include "AliMUONFastTracking.h"
#include <TRandom.h>
AliFastMuonTrackingRes::AliFastMuonTrackingRes() :
AliFastResponse("Resolution", "Muon Tracking Resolution")
{
+// Deafault constructor
SetBackground();
}
void AliFastMuonTrackingRes::Init()
{
+// Initialisation
fFastTracking = AliMUONFastTracking::Instance();
fFastTracking->Init(fBackground);
}
void AliFastMuonTrackingRes::Evaluate(Float_t p, Float_t theta , Float_t phi,
Float_t& pS, Float_t& thetaS, Float_t& phiS)
{
+//
+// Evaluate Gaussian smearing from given kinematics
+//
Double_t meanp = fFastTracking->MeanP (p, theta, phi, Int_t(fCharge));
Double_t sigmap = fFastTracking->SigmaP (p, theta, phi, Int_t(fCharge));
-#ifndef ALIFASTMUONTRACKINGRES
-#define ALIFASTMUONTRACKINGRES
+#ifndef ALIFASTMUONTRACKINGRES_H
+#define ALIFASTMUONTRACKINGRES_H
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
/* $Id$ */
+// Implementation of AliFastResponse for the Muon Spectrometer resolution.
+// Author: andreas.morsch@cern.ch
+//
#include "AliFastResponse.h"
class AliMUONFastTracking;
class AliFastMuonTrackingRes : public AliFastResponse {
public:
AliFastMuonTrackingRes();
- ~AliFastMuonTrackingRes(){;}
+ virtual ~AliFastMuonTrackingRes(){;}
void SetBackground(Float_t bg = 1.) {fBackground = bg;}
void SetCharge(Float_t charge = 1.) {fCharge = charge;}
virtual void Init();
/*
$Log$
+Revision 1.6 2003/08/12 15:16:25 morsch
+Saver initialisation of fFitp array. (Lenaic COUEDEL)
+
Revision 1.5 2003/08/05 16:14:20 morsch
Some problems with too big fluctuations corrected. (A. de Falco)
*/
+//-------------------------------------------------------------------------
+// Class AliMUONFastTracking
+//
+// Manager for the fast simulation of tracking in the muon spectrometer
+// This class reads the lookup tables containing the parameterization
+// of the deltap, deltatheta, deltaphi for different background levels
+// and provides the related smeared parameters.
+// Used by AliFastMuonTrackingEff, AliFastMuonTrackingAcc,
+// AliFastMuonTrackingRes.
+//-------------------------------------------------------------------------
+
#include "AliMUONFastTracking.h"
#include "AliMUONFastTrackingEntry.h"
-#include <TMatrixD.h>
#include <TSpline.h>
#include <TFile.h>
-#include <TH1.h>
#include <TH3.h>
#include <TF1.h>
#include <TRandom.h>
ClassImp(AliMUONFastTracking)
+
AliMUONFastTracking* AliMUONFastTracking::fgMUONFastTracking=NULL;
static Double_t FitP(Double_t *x, Double_t *par){
return value;
}
+AliMUONFastTracking::AliMUONFastTracking(const AliMUONFastTracking & ft):TObject()
+{
+// Copy constructor
+ ft.Copy(*this);
+}
+
+
AliMUONFastTracking* AliMUONFastTracking::Instance()
{
// Set random number generator
AliMUONFastTracking::AliMUONFastTracking()
{
+//
+// constructor
+//
for (Int_t i = 0; i<20;i++) {
for (Int_t j = 0; j<20; j++) {
for (Int_t k = 0; k<20; k++) {
char filename [100];
if (fClusterFinder==kOld) sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root");
- else sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root");
+ else sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT-AZ.root");
TFile *file = new TFile(filename);
ReadLUT(file);
void AliMUONFastTracking::ReadLUT(TFile* file)
{
+ //
+ // read the lookup tables from file
+ //
TH3F *heff[5][3], *hacc[5][3], *hmeanp, *hsigmap, *hsigma1p, *hchi2p;
TH3F *hnormg2, *hmeang2, *hsigmag2, *hmeantheta, *hsigmatheta, *hchi2theta;
TH3F *hmeanphi, *hsigmaphi, *hchi2phi;
printf ("Reading parameters from LUT file %s...\n",file->GetName());
- const Float_t bkg[4] = {0, 0.5, 1, 2};
+ const Float_t kBkg[4] = {0, 0.5, 1, 2};
for (Int_t ibkg=0; ibkg<4; ibkg++) {
- sprintf (tag,"BKG%g",bkg[ibkg]);
+ sprintf (tag,"BKG%g",kBkg[ibkg]);
file->cd(tag);
for (Int_t isplp = 0; isplp<kSplitP; isplp++) {
for (Int_t ispltheta = 0; ispltheta<kSplitTheta; ispltheta++) {
void AliMUONFastTracking::GetBinning(Int_t &nbinp, Float_t &pmin, Float_t &pmax,
Int_t &nbintheta, Float_t &thetamin,
Float_t &thetamax,
- Int_t &nbinphi, Float_t &phimin, Float_t &phimax)
+ Int_t &nbinphi, Float_t &phimin, Float_t &phimax) const
{
+ //
+ // gets the binning for the discrete parametrizations in the lookup table
+ //
nbinp = fNbinp;
pmin = fPmin;
pmax = fPmax;
void AliMUONFastTracking::GetIpIthetaIphi(Float_t p, Float_t theta, Float_t phi,
Int_t charge, Int_t &ip, Int_t &itheta,
- Int_t &iphi)
+ Int_t &iphi) const
{
+ //
+ // gets the id of the cells in the LUT for a given (p,theta,phi, charge)
+ //
if (charge < 0) phi = -phi;
ip = Int_t (( p - fPmin ) / fDeltaP);
itheta = Int_t (( theta - fThetamin ) / fDeltaTheta);
}
void AliMUONFastTracking::GetSplit(Int_t ip, Int_t itheta,
- Int_t &nSplitP, Int_t &nSplitTheta) {
+ Int_t &nSplitP, Int_t &nSplitTheta) const
+{
+ //
+ // the first cell is splitted in more bins for theta and momentum
+ // parameterizations. Get the number of divisions for the splitted bins
+ //
if (ip==0) nSplitP = 5;
else nSplitP = 2;
if (itheta==0) nSplitTheta = 3;
Float_t AliMUONFastTracking::Efficiency(Float_t p, Float_t theta,
Float_t phi, Int_t charge){
+ //
+ // gets the tracking efficiency
+ //
Int_t ip=0, itheta=0, iphi=0;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
Int_t nSplitP, nSplitTheta;
Float_t AliMUONFastTracking::Acceptance(Float_t p, Float_t theta,
Float_t phi, Int_t charge){
+ //
+ // gets the geometrical acceptance
+ //
if (theta<fThetamin || theta>fThetamax) return 0;
Int_t ip=0, itheta=0, iphi=0;
}
Float_t AliMUONFastTracking::MeanP(Float_t p, Float_t theta,
- Float_t phi, Int_t charge)
+ Float_t phi, Int_t charge) const
{
+ //
+ // gets the mean value of the prec-pgen distribution
+ //
Int_t ip=0, itheta=0, iphi=0;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
return fCurrentEntry[ip][itheta][iphi]->fMeanp;
}
Float_t AliMUONFastTracking::SigmaP(Float_t p, Float_t theta,
- Float_t phi, Int_t charge)
+ Float_t phi, Int_t charge) const
{
+ //
+ // gets the width of the prec-pgen distribution
+ //
Int_t ip=0, itheta=0, iphi=0;
Int_t index;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
}
Float_t AliMUONFastTracking::Sigma1P(Float_t p, Float_t theta,
- Float_t phi, Int_t charge)
+ Float_t phi, Int_t charge) const
{
+ //
+ // gets the width correction of the prec-pgen distribution (see FitP)
+ //
Int_t ip=0, itheta=0, iphi=0;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
if (p>fPmax) {
}
Float_t AliMUONFastTracking::NormG2(Float_t p, Float_t theta,
- Float_t phi, Int_t charge)
+ Float_t phi, Int_t charge) const
{
+ //
+ // gets the relative normalization of the background
+ // (gaussian) component in the prec-pgen distribution
+ //
Int_t ip=0, itheta=0, iphi=0;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
if (p>fPmax) {
}
Float_t AliMUONFastTracking::MeanG2(Float_t p, Float_t theta,
- Float_t phi, Int_t charge)
+ Float_t phi, Int_t charge) const
{
+ //
+ // gets the mean value of the background
+ // (gaussian) component in the prec-pgen distribution
+ //
Int_t ip=0, itheta=0, iphi=0;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
if (p>fPmax) {
}
Float_t AliMUONFastTracking::SigmaG2(Float_t p, Float_t theta,
- Float_t phi, Int_t charge)
+ Float_t phi, Int_t charge) const
{
+ //
+ // gets the width of the background
+ // (gaussian) component in the prec-pgen distribution
+ //
Int_t ip=0, itheta=0, iphi=0;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
if (p>fPmax) {
Float_t AliMUONFastTracking::MeanTheta(Float_t p, Float_t theta,
- Float_t phi, Int_t charge)
+ Float_t phi, Int_t charge) const
{
+ //
+ // gets the mean value of the thetarec-thetagen distribution
+ //
Int_t ip=0, itheta=0, iphi=0;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
return fCurrentEntry[ip][itheta][iphi]->fMeantheta;
}
-Float_t AliMUONFastTracking::SigmaTheta(Float_t p, Float_t theta,
- Float_t phi, Int_t charge){
+Float_t AliMUONFastTracking::SigmaTheta(Float_t p, Float_t theta,
+ Float_t phi, Int_t charge) const
+{
+ //
+ // gets the width of the thetarec-thetagen distribution
+ //
Int_t ip=0, itheta=0, iphi=0;
Int_t index;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
Float_t AliMUONFastTracking::MeanPhi(Float_t p, Float_t theta,
- Float_t phi, Int_t charge){
+ Float_t phi, Int_t charge) const
+{
+ //
+ // gets the mean value of the phirec-phigen distribution
+ //
Int_t ip=0, itheta=0, iphi=0;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
return fCurrentEntry[ip][itheta][iphi]->fMeanphi;
Float_t AliMUONFastTracking::SigmaPhi(Float_t p, Float_t theta,
Float_t phi, Int_t charge){
+ //
+ // gets the width of the phirec-phigen distribution
+ //
Int_t ip=0, itheta=0, iphi=0;
Int_t index;
GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
}
void AliMUONFastTracking::SetSpline(){
+ //
+ // sets the spline functions for a smooth behaviour of the parameters
+ // when going from one cell to another
+ //
printf ("Setting spline functions...");
char splname[40];
Double_t x[20][3];
}
void AliMUONFastTracking::SetBackground(Float_t bkg){
+ //
// linear interpolation of the parameters in the LUT between 2 values where
// the background has been actually calculated
-
+ //
if (bkg>2) printf ("WARNING: unsafe extrapolation!\n");
fBkg = bkg;
- Float_t BKG[4] = {0, 0.5, 1, 2}; // bkg values for which LUT is calculated
+ Float_t bkgLevel[4] = {0, 0.5, 1, 2}; // bkg values for which LUT is calculated
Int_t ibkg;
- for (ibkg=0; ibkg<4; ibkg++) if ( bkg < BKG[ibkg]) break;
+ for (ibkg=0; ibkg<4; ibkg++) if ( bkg < bkgLevel[ibkg]) break;
if (ibkg == 4) ibkg--;
if (ibkg == 0) ibkg++;
- Float_t x0 = BKG[ibkg-1];
- Float_t x1 = BKG[ibkg];
+ Float_t x0 = bkgLevel[ibkg-1];
+ Float_t x1 = bkgLevel[ibkg];
Float_t x = (bkg - x0) / (x1 - x0);
Float_t y0, y1;
}
TF1* AliMUONFastTracking::GetFitP(Int_t ip,Int_t itheta,Int_t iphi) {
+ // gets the correct prec-pgen distribution for a given LUT cell
if (!fFitp[ip][itheta][iphi]) {
fFitp[ip][itheta][iphi] = new TF1("fit1",FitP,-20.,20.,6);
fFitp[ip][itheta][iphi]->SetNpx(500);
return fFitp[ip][itheta][iphi];
}
- // to guarantee a safe extrapolation for sigmag2 to 0<bkg<0.5, let's fit
- // with a straight line sigmag2 vs bkg for bkg=0.5, 1 and 2, and put the
- // sigma2(BKG=0) as the extrapolation of this fit
+AliMUONFastTracking& AliMUONFastTracking::operator=(const AliMUONFastTracking& rhs)
+{
+// Assignment operator
+ rhs.Copy(*this);
+ return *this;
+}
+
+void AliMUONFastTracking::Copy(AliMUONFastTracking&) const
+{
+ //
+ // Copy
+ //
+ Fatal("Copy","Not implemented!\n");
+}
-#ifndef ALIMUONFASTTRACKING
-#define ALIMUONFASTTRACKING
+#ifndef ALIMUONFASTTRACKING_H
+#define ALIMUONFASTTRACKING_H
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
/* $Id$ */
+
+//-------------------------------------------------------------------------
+// Class AliMUONFastTracking
+//
+// Manager for the fast simulation of tracking in the muon spectrometer
+// This class reads the lookup tables containing the parameterization
+// of the deltap, deltatheta, deltaphi for different background levels
+// and provides the related smeared parameters
+//-------------------------------------------------------------------------
+
class TF1;
class TSpline3;
class TFile;
void ReadLUT(TFile *file);
void GetBinning(Int_t &nbinp, Float_t &pmin, Float_t &pmax,
Int_t &nbintheta, Float_t &thetamin, Float_t &thetamax,
- Int_t &nbinphi, Float_t &phimin, Float_t &phimax);
+ Int_t &nbinphi, Float_t &phimin, Float_t &phimax) const;
void GetIpIthetaIphi(Float_t p, Float_t theta, Float_t phi, Int_t charge,
- Int_t &ip, Int_t &itheta, Int_t &iphi);
- void GetSplit(Int_t ip, Int_t itheta, Int_t &nSplitP, Int_t &nSplitTheta);
+ Int_t &ip, Int_t &itheta, Int_t &iphi) const;
+ void GetSplit(Int_t ip, Int_t itheta, Int_t &nSplitP, Int_t &nSplitTheta) const;
Float_t Efficiency(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t Acceptance(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t MeanP(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t SigmaP(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t Sigma1P(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t NormG2(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t MeanG2(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t SigmaG2(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t MeanTheta(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t SigmaTheta(Float_t p, Float_t theta, Float_t phi, Int_t charge);
- Float_t MeanPhi(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+ Float_t Acceptance(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+ Float_t MeanP(Float_t p, Float_t theta, Float_t phi, Int_t charge) const;
+ Float_t SigmaP(Float_t p, Float_t theta, Float_t phi, Int_t charge) const;
+ Float_t Sigma1P(Float_t p, Float_t theta, Float_t phi, Int_t charge) const;
+ Float_t NormG2(Float_t p, Float_t theta, Float_t phi, Int_t charge) const;
+ Float_t MeanG2(Float_t p, Float_t theta, Float_t phi, Int_t charge) const;
+ Float_t SigmaG2(Float_t p, Float_t theta, Float_t phi, Int_t charge) const;
+ Float_t MeanTheta(Float_t p, Float_t theta, Float_t phi, Int_t charge) const;
+ Float_t SigmaTheta(Float_t p, Float_t theta, Float_t phi, Int_t charge)const;
+ Float_t MeanPhi(Float_t p, Float_t theta, Float_t phi, Int_t charge) const;
Float_t SigmaPhi(Float_t p, Float_t theta, Float_t phi, Int_t charge);
void SetSpline();
- Float_t GetBackground() {return fBkg;}
+ Float_t GetBackground() const {return fBkg;}
void SetLUTClusterFinder(LUTClusterType clusterFinder) { fClusterFinder = clusterFinder;}
void SetBackground(Float_t bkg);
void UseSpline (Int_t splineSwitch=1) {fSpline = splineSwitch;}
TF1* GetFitP(Int_t ip, Int_t itheta, Int_t iphi);
- private:
- AliMUONFastTracking();
- AliMUONFastTracking(Float_t bkg){;}
protected:
- Int_t fNentries;
- Int_t fNbinp;
- Float_t fPmin;
- Float_t fPmax;
- Float_t fDeltaP;
- Int_t fNbintheta;
- Float_t fThetamin;
- Float_t fThetamax;
- Float_t fDeltaTheta;
- Int_t fNbinphi;
- Float_t fPhimin;
- Float_t fPhimax;
- Float_t fDeltaPhi;
- Int_t fPrintLevel;
- Float_t fBkg;
+ Int_t fNbinp; // n. of momentum bins in the lookup table
+ Float_t fPmin; // min. value of momentum parameterized in LUT
+ Float_t fPmax; // max. value of momentum parameterized in LUT
+ Float_t fDeltaP; // momentum bin width
+ Int_t fNbintheta; // n. of theta bins in the lookup table
+ Float_t fThetamin; // min. value of theta parameterized in LUT
+ Float_t fThetamax; // max. value of theta parameterized in LUT
+ Float_t fDeltaTheta; // theta bin width
+ Int_t fNbinphi; // n. of phi bins in the lookup table
+ Float_t fPhimin; // min. value of phi parameterized in LUT
+ Float_t fPhimax; // min. value of phi parameterized in LUT
+ Float_t fDeltaPhi; // phi bin width
+ Int_t fPrintLevel; // level of information printed for debugging
+ Float_t fBkg; // soft background level
TF1 *fFitp[20][20][20]; // func for psmear-pgen distr
AliMUONFastTrackingEntry *fEntry[20][20][20][4]; // array of LUT parameters
AliMUONFastTrackingEntry *fCurrentEntry[20][20][20]; // array of LUT parameters
- public:
- TSpline3 *fSplineEff[200][3]; // spline funcs for efficiency
- TSpline3 *fSplineAcc[200][3]; // spline funcs for acceptance
- TSpline3 *fSplineSigmap[200][3]; //
- TSpline3 *fSplineSigma1p[200][3]; //!
- TSpline3 *fSplineSigmatheta[200][3]; //!
- TSpline3 *fSplineSigmaphi[200][3]; //!
- protected:
- Int_t fSpline;
- LUTClusterType fClusterFinder;
+ TSpline3 *fSplineEff[200][3]; // spline funcs for efficiency
+ TSpline3 *fSplineAcc[200][3]; // spline funcs for acceptance
+ TSpline3 *fSplineSigmap[200][3]; // spl.funcs for dp distribution width
+ TSpline3 *fSplineSigma1p[200][3]; // spl.funcs for dp distr. width correction (see function FitP)
+ TSpline3 *fSplineSigmatheta[200][3]; // spl.funcs for dtheta distr. width
+ TSpline3 *fSplineSigmaphi[200][3]; // spl.funcs for dphi distr. width
+ Int_t fSpline; // switches on/off the use of spline
+ LUTClusterType fClusterFinder; // type of cluster finder (old/new)
static AliMUONFastTracking* fgMUONFastTracking; //!Pointer to single instance
- ClassDef(AliMUONFastTracking,1) // Fast MUON Tracking Data Handler
+ ClassDef(AliMUONFastTracking,1) // Fast MUON Tracking Data Handler
+ private:
+ AliMUONFastTracking();
+ AliMUONFastTracking(Float_t bkg){;}
+ AliMUONFastTracking(const AliMUONFastTracking &ft);
+ void Copy(AliMUONFastTracking &) const;
+ AliMUONFastTracking& operator=(const AliMUONFastTracking & rhs);
};
#endif