+
+ const Double_t kphi = TMath::Pi() * gRandom->Rndm();
+ const Double_t kb1 = fgBMax * gRandom->Rndm();
+ y += fgWStarfi->Eval(kb1, kphi);
+ }
+ y *= 2. * TMath::Pi() * fgBMax / fgkMCInts; //fm^-2
+ y *= ka * ka * 0.1; //mbarn^-1
+ return y;
+}
+
+Double_t AliFastGlauber::WKParticipants(const Double_t* x, const Double_t* par)
+{
+ //
+ // Kernel for number of participants
+ // as a function of r and phi
+ //
+ const Double_t kr1 = x[0];
+ const Double_t kphi = x[1];
+ const Double_t kb = par[0]; //fm
+ const Double_t ksig = par[1]; //mbarn
+ const Double_t ka = par[2]; //mass number
+ const Double_t kr2 = TMath::Sqrt(kr1*kr1 +kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
+ const Double_t kxsi = fgWSta->Eval(kr2) * ksig * 0.1; //no units
+ /*
+ Double_t y=(1-TMath::Power((1-xsi),aa))
+ */
+ Double_t a = ka;
+ Double_t sum = ka * kxsi;
+ Double_t y = sum;
+ for (Int_t i = 1; i <= ka; i++)
+ {
+ a--;
+ sum *= (-kxsi) * a / Float_t(i+1);
+ y += sum;
+ }
+ y *= kr1 * fgWSta->Eval(kr1);
+ return y; //fm^-1
+}
+
+Double_t AliFastGlauber::WParticipants(const Double_t* x, const Double_t* par)
+{
+ //
+ // Number of Participants as
+ // a function of b
+ //
+ const Double_t kb = x[0];
+ const Double_t ksig = par[0]; //mbarn
+ const Double_t ka = par[1]; //mass number
+ fgWKParticipants->SetParameter(0, kb);
+ fgWKParticipants->SetParameter(1, ksig);
+ fgWKParticipants->SetParameter(2, ka);
+
+ //
+ // MC Integration
+ //
+ Double_t y = 0;
+ for (Int_t i = 0; i < fgkMCInts; i++)
+ {
+ const Double_t kphi = TMath::Pi() * gRandom->Rndm();
+ const Double_t kb1 = fgBMax * gRandom->Rndm();
+ y += fgWKParticipants->Eval(kb1, kphi);
+ }
+ y *= 2. * ka * 2. * TMath::Pi() * fgBMax / fgkMCInts;
+ return y; //no units
+}
+
+Double_t AliFastGlauber::WSgeo(const Double_t* x, const Double_t* par)
+{
+ //
+ // Geometrical Cross-Section
+ // as a function of b
+ //
+ const Double_t kb = x[0]; //fm
+ const Double_t ksigNN = par[0]; //mbarn
+ const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
+ Double_t y = 2. * TMath::Pi() * kb * (1. - TMath::Exp(- ksigNN * ktaa));
+ return y; //fm
+}
+
+Double_t AliFastGlauber::WSbinary(const Double_t* x, const Double_t* par)
+{
+ //
+ // Number of binary hard collisions
+ // as a function of b
+ //
+ const Double_t kb = x[0]; //fm
+ const Double_t ksig = par[0]; //mbarn
+ const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
+ Double_t y = 2. * TMath::Pi() * kb * ksig * ktaa;
+ return y; //fm
+}
+
+Double_t AliFastGlauber::WSN(const Double_t* x, const Double_t* /*par*/)
+{
+ //
+ // Number of hard processes per event
+ // as a function of b
+ const Double_t kb = x[0];
+ Double_t y = fgWSbinary->Eval(kb)/fgWSgeo->Eval(kb);
+ return y; //no units
+}
+
+Double_t AliFastGlauber::WEnergyDensity(const Double_t* x, const Double_t* par)
+{
+ //
+ // Initial energy density
+ // as a function of the impact parameter
+ //
+ const Double_t kb = x[0];
+ const Double_t krA = par[0];
+ //
+ // Attention: area of transverse reaction zone in hard-sphere approximation !
+ const Double_t krA2=krA*krA;
+ const Double_t kb2=kb*kb;
+ const Double_t ksaa = (TMath::Pi() - 2. * TMath::ASin(kb/ 2./ krA)) * krA2
+ - kb * TMath::Sqrt(krA2 - kb2/ 4.); //fm^2
+ const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
+ Double_t y=ktaa/ksaa*10;
+ return y; //fm^-4
+}
+
+Double_t AliFastGlauber::WAlmond(const Double_t* x, const Double_t* par)
+{
+ //
+ // Almond shaped interaction region
+ // as a function of cartesian x,y.
+ //
+ const Double_t kb = par[0];
+ const Double_t kxx = x[0] + kb/2.;
+ const Double_t kyy = x[1];
+ const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
+ const Double_t kphi = TMath::ATan2(kyy,kxx);
+ const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
+ //
+ // Interaction probability calculated as product of thicknesses
+ //
+ Double_t y = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
+ return y; //fm^-4
+}
+
+Double_t AliFastGlauber::WIntRadius(const Double_t* x, const Double_t* par)
+{
+ //
+ // Average interaction density over radius
+ // at which interaction takes place
+ // as a function of radius
+ //
+ const Double_t kr = x[0];
+ const Double_t kb = par[0];
+ fgWAlmond->SetParameter(0, kb);
+ // Average over phi in small steps
+ const Double_t kdphi = 2. * TMath::Pi() / 100.;
+ Double_t phi = 0.;
+ Double_t y = 0.;
+ for (Int_t i = 0; i < 100; i++) {
+ const Double_t kxx = kr * TMath::Cos(phi);
+ const Double_t kyy = kr * TMath::Sin(phi);
+ y += fgWAlmond->Eval(kxx,kyy);
+ phi += kdphi;
+ } // phi loop
+ // Result multiplied by Jacobian (2 pi r)
+ y *= 2. * TMath::Pi() * kr / 100.;
+ return y; //fm^-3
+}
+
+Double_t AliFastGlauber::WPathLength0(const Double_t* x, const Double_t* par)
+{
+ //
+ // Path Length as a function of phi
+ // for interaction point fixed at (0,0)
+ // as a function of phi-direction
+ //
+ // Phi direction in Almond
+ const Double_t kphi0 = x[0];
+ const Double_t kb = par[0];
+ // Path Length definition
+ const Int_t kiopt = Int_t(par[1]);
+
+ // Step along radial direction phi
+ const Int_t kNp = 100; // Steps in r
+ const Double_t kDr = fgBMax/kNp;
+ Double_t r = 0.;
+ Double_t rw = 0.;
+ Double_t w = 0.;
+ for (Int_t i = 0; i < kNp; i++) {
+ //
+ // Transform into target frame
+ //
+ const Double_t kxx = r * TMath::Cos(kphi0) + kb / 2.;
+ const Double_t kyy = r * TMath::Sin(kphi0);
+ const Double_t kphi = TMath::ATan2(kyy, kxx);
+ const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
+ // Radius in projectile frame
+ const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
+ const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
+
+ rw += ky * r;
+ w += ky;
+ r += kDr;
+ } // radial steps
+
+ Double_t y=0.;
+ if (!kiopt) { // My length definition (is exact for hard disk)
+ if(w) y= 2. * rw / w;
+ } else {
+ const Double_t knorm=fgWSta->Eval(1e-4);
+ if(knorm) y = TMath::Sqrt(2. * rw * kDr / knorm / knorm);
+ }
+ return y; //fm
+}
+
+Double_t AliFastGlauber::WPathLength(const Double_t* x, const Double_t* par)
+{
+ //
+ // Path Length as a function of phi
+ // Interaction point from random distribution
+ // as a function of the phi-direction
+ const Double_t kphi0 = x[0];
+ const Double_t kb = par[0];
+ fgWAlmond->SetParameter(0, kb);
+ const Int_t kNpi = Int_t (par[1]); //Number of interactions
+ const Int_t kiopt = Int_t(par[2]); //Path Length definition
+
+ //
+ // r-steps
+ //
+ const Int_t kNp = 100;
+ const Double_t kDr = fgBMax/Double_t(kNp);
+ Double_t l = 0.; // Path length
+ for (Int_t in = 0; in < kNpi; in ++) {
+ Double_t rw = 0.;
+ Double_t w = 0.;
+ // Interaction point
+ Double_t x0, y0;
+ fgWAlmond->GetRandom2(x0, y0);
+ // Initial radius
+ const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
+ const Int_t knps = Int_t ((fgBMax - kr0)/kDr) - 1;
+
+ // Radial steps
+ Double_t r = 0.;
+ for (Int_t i = 0; (i < knps ); i++) {
+ // Transform into target frame
+ const Double_t kxx = x0 + r * TMath::Cos(kphi0) + kb / 2.;
+ const Double_t kyy = y0 + r * TMath::Sin(kphi0);
+ const Double_t kphi = TMath::ATan2(kyy, kxx);
+ const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
+ // Radius in projectile frame
+ const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
+ const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
+
+ rw += ky * r;
+ w += ky;
+ r += kDr;
+ } // steps
+ // Average over interactions
+ if (!kiopt) {
+ if(w) l += (2. * rw / w);
+ } else {
+ const Double_t knorm=fgWSta->Eval(1e-4);
+ if(knorm) l+= 2. * rw * kDr / knorm / knorm;