PathLength functions added.

diff --git a/FASTSIM/AliFastGlauber.cxx b/FASTSIM/AliFastGlauber.cxx
index cffcdabedf0a22b6650a0a644faa9227b8e110a6..c533fec364a71755f86541a5ea06e250a6261981 100644 (file)
--- a/FASTSIM/AliFastGlauber.cxx
+++ b/FASTSIM/AliFastGlauber.cxx
@@ -27,16 +27,20 @@
 
 ClassImp(AliFastGlauber)
 
-TF1*    AliFastGlauber::fWSb      = NULL;     
-TF2*    AliFastGlauber::fWSbz     = NULL;    
-TF1*    AliFastGlauber::fWSz      = NULL;     
-TF1*    AliFastGlauber::fWSta     = NULL;    
-TF2*    AliFastGlauber::fWStarfi  = NULL; 
-TF1*    AliFastGlauber::fWStaa    = NULL;   
-TF1*    AliFastGlauber::fWSgeo    = NULL;   
-TF1*    AliFastGlauber::fWSbinary = NULL;   
-TF1*    AliFastGlauber::fWSN      = NULL;   
-Float_t AliFastGlauber::fbMax     = 0.;
+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::fWIntRadius   = NULL;   
+Float_t AliFastGlauber::fbMax         = 0.;
 
 AliFastGlauber::AliFastGlauber()
 {
@@ -85,6 +89,26 @@ void AliFastGlauber::Init(Int_t mode)
     fWStarfi->SetNpx(200);     
     fWStarfi->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);    
+//
+//  Path Length as a function of Phi
+//    
+    fWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 1);
+    fWPathLength0->SetParameter(0, 0.);     
+    fWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 2);
+//  Impact Parameter
+    fWPathLength->SetParameter(0, 0.);    
+//  Number of interactions used for average
+    fWPathLength->SetParameter(1, 1000.);    
+
+    fWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fbMax, 1);
+    fWIntRadius->SetParameter(0, 0.);     
+//
 //  Overlap
 //
     if (! mode) {
@@ -173,16 +197,69 @@ void AliFastGlauber::DrawN()
     fWSN->Draw();
 }
 
-void AliFastGlauber::DrawKernel()
+void AliFastGlauber::DrawKernel(Double_t b)
 {
 //
 //  Draw Kernel
 //
     TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
     c6->cd();
+    fWStarfi->SetParameter(0, b);
     fWStarfi->Draw();
 }
 
+void AliFastGlauber::DrawAlmond(Double_t b)
+{
+//
+//  Draw Kernel
+//
+    TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
+    c7->cd();
+    fWAlmond->SetParameter(0, b);
+    fWAlmond->Draw();
+}
+
+void AliFastGlauber::DrawPathLength0(Double_t b)
+{
+//
+//  Draw Kernel
+//
+    TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
+    c8->cd();
+    fWPathLength0->SetParameter(0, b);
+    fWPathLength0->SetMinimum(0.); 
+    fWPathLength0->SetMaximum(10.); 
+    fWPathLength0->Draw();
+}
+
+void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni)
+{
+//
+//  Draw Kernel
+//
+    TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
+    c9->cd();
+    fWAlmond->SetParameter(0, b);
+
+    fWPathLength->SetParameter(0, b);
+    fWPathLength->SetParameter(1, Double_t (ni));
+    fWPathLength->SetMinimum(0.); 
+    fWPathLength->SetMaximum(10.); 
+    fWPathLength->Draw();
+}
+
+void AliFastGlauber::DrawIntRadius(Double_t b)
+{
+//
+//  Draw Kernel
+//
+    TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
+    c10->cd();
+    fWIntRadius->SetParameter(0, b);
+    fWIntRadius->SetMinimum(0.);
+    fWIntRadius->Draw();
+}
+
 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
 {
 //
@@ -261,6 +338,149 @@ Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
 }
 
 
+Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
+{
+//
+//  Almond shaped interaction region
+//
+    Double_t b    = par[0];
+    Double_t xx   = x[0] + b/2.;
+    Double_t yy   = x[1];
+    Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
+    Double_t phi  = TMath::ATan2(yy,xx);
+    
+    Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
+//
+//  Interaction probability calculated as product of thicknesses
+//
+    Double_t y    = fWSta->Eval(r1) * fWSta->Eval(r2);
+    return y;
+}
+
+Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
+{
+//
+//  Average radius at which interaction takes place
+//
+//  Radius in the Almond
+    Double_t r    = x[0];
+//  Impact parameter
+    Double_t b    = par[0];
+    fWAlmond->SetParameter(0, b);
+//  Steps in phi
+    Double_t dphi = 2. * TMath::Pi() / 100.;
+//  Average over phi    
+    Double_t phi  = 0.;
+    Double_t y    = 0.;
+
+    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);
+       phi += dphi;
+    } // phi loop
+// Result multiplied by Jacobian (2 pi r)     
+    return (2. * TMath::Pi() * y * r / 100.);
+}
+
+Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
+{
+//
+//  Path Length as a function of phi for interaction point fixed at (0,0)
+//
+//
+//  Steps in r 
+    const Int_t    np  = 100;
+    const Double_t dr  = fbMax/Double_t(np);
+//  Impact parameter    
+    Double_t b      = par[0];
+//  Phi direction in Almond
+    Double_t phi0   = x[0];
+    Double_t r  = 0.;
+    Double_t rw = 0.;
+    Double_t w  = 0.;
+//  Step along radial direction phi   
+    for (Int_t i = 0; i < np; i++) {
+//
+//  Transform into target frame
+//
+       Double_t xx   = r * TMath::Cos(phi0) + b / 2.;
+       Double_t yy   = r * TMath::Sin(phi0);
+       Double_t phi  = TMath::ATan2(yy, xx);
+       
+       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);
+
+       rw += y * r;
+       w  += y;
+       r  += dr;
+    } // radial steps
+//
+//  My length definition (is exact for hard sphere)
+    return (2. * rw / w);
+}
+
+Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
+{
+//
+//  Path Length as a function of phi 
+//  Interaction point from random distribution
+//
+//
+//  r-steps
+// 
+    const Int_t    np   = 100;
+    const Double_t dr  = fbMax/Double_t(np);
+//  Number of interactions
+    const Int_t    npi  = Int_t (par[1]);
+
+//
+//  Impact parameter    
+    Double_t b      = par[0];
+//  Phi direction
+    Double_t phi0   = x[0];
+
+    printf("phi0 %f \n", phi0);
+    
+//  Path length 
+    Double_t l = 0.;
+    
+    for (Int_t in = 0; in < npi; in ++) {
+       Double_t rw = 0.;
+       Double_t w  = 0.;
+       
+       // Interaction point
+       Double_t x0, y0;
+       fWAlmond->GetRandom2(x0, y0);
+// Initial radius
+       Double_t r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
+       Int_t    nps = Int_t ((fbMax - r0)/dr) - 1;
+       
+       Double_t r  = 0.;
+// Radial steps
+       for (Int_t i = 0; (i < nps ); i++) {
+           
+// Transform into target frame
+           Double_t xx   = x0 + r * TMath::Cos(phi0) + b / 2.;
+           Double_t yy   = y0 + r * TMath::Sin(phi0);
+           Double_t phi  = TMath::ATan2(yy, xx);
+           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);
+           
+           rw += y * r;
+           w  += y;
+           r  += dr;
+       } // steps
+// Average over interactions
+       l += 2. * rw / w;
+    } // interactions
+    return (l / Double_t(npi));
+}
+
 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
 {
 //

diff --git a/FASTSIM/AliFastGlauber.h b/FASTSIM/AliFastGlauber.h
index 3670fd18bd353163783967abad683873d17a1366..df1e4d46a3cbc8641cf18ef637c14572e8c240e7 100644 (file)
--- a/FASTSIM/AliFastGlauber.h
+++ b/FASTSIM/AliFastGlauber.h
@@ -18,15 +18,20 @@ class AliFastGlauber : public TObject {
     void SetMaxImpact(Float_t bmax = 20.) {fbMax = bmax;};
     void SetHardCrossSection(Float_t xs = 6.6) {fSigmaHard = xs;}
     
-    static Double_t WSb      (Double_t *xx, Double_t *par);
-    static Double_t WSbz     (Double_t *xx, Double_t *par);
-    static Double_t WSz      (Double_t *xx, Double_t *par);
-    static Double_t WSta     (Double_t *xx, Double_t *par);
-    static Double_t WStarfi  (Double_t *xx, Double_t *par);
-    static Double_t WStaa    (Double_t *xx, Double_t *par);
-    static Double_t WSgeo    (Double_t *xx, Double_t *par);
-    static Double_t WSbinary (Double_t *xx, Double_t *par);
-    static Double_t WSN      (Double_t *xx, Double_t *par);
+    static Double_t WSb          (Double_t *xx, Double_t *par);
+    static Double_t WSbz         (Double_t *xx, Double_t *par);
+    static Double_t WSz          (Double_t *xx, Double_t *par);
+    static Double_t WSta         (Double_t *xx, Double_t *par);
+    static Double_t WStarfi      (Double_t *xx, Double_t *par);
+    static Double_t WStaa        (Double_t *xx, Double_t *par);
+    static Double_t WSgeo        (Double_t *xx, Double_t *par);
+    static Double_t WSbinary     (Double_t *xx, Double_t *par);
+    static Double_t WSN          (Double_t *xx, Double_t *par);
+    static Double_t WAlmond      (Double_t *xx, Double_t *par);
+    static Double_t WPathLength0 (Double_t *xx, Double_t *par);
+    static Double_t WPathLength  (Double_t *xx, Double_t *par);
+    static Double_t WIntRadius   (Double_t *xx, Double_t *par);
+    
     void Init(Int_t mode = 0);
     void DrawWSb();
     void DrawThickness();
@@ -34,7 +39,12 @@ class AliFastGlauber : public TObject {
     void DrawGeo();
     void DrawBinary();
     void DrawN();    
-    void DrawKernel();
+    void DrawKernel(Double_t b);
+    void DrawAlmond(Double_t b);
+    void DrawPathLength0(Double_t b);
+    void DrawPathLength(Double_t b, Int_t ni = 1000);
+    void DrawIntRadius(Double_t b);
+    
     Double_t CrossSection(Double_t b1, Double_t b2);
     Double_t FractionOfHardCrossSection(Double_t b1, Double_t b2);
     Double_t Binaries(Double_t b);
@@ -45,15 +55,19 @@ class AliFastGlauber : public TObject {
     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 TF1*    fWSgeo;    // dSigma/db geometric
-    static TF1*    fWSbinary; // dSigma/db binary
-    static TF1*    fWSN;      // dN/db binary
+    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
     
     Float_t fWSr0;      // Wood-Saxon Parameter r0
     Float_t fWSd;       // Wood-Saxon Parameter d