Re-organised scripts
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Jul 2011 08:32:12 +0000 (08:32 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Jul 2011 08:32:12 +0000 (08:32 +0000)
PWG2/FORWARD/analysis2/tests/TestAcc.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/tests/TestELossDist.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/tests/TestFitELoss.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/tests/TestMakeELossFits.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/tests/TestMarkers.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/tests/TestPoisson.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/tests/TestRunMakeELossFit.C [new file with mode: 0644]

diff --git a/PWG2/FORWARD/analysis2/tests/TestAcc.C b/PWG2/FORWARD/analysis2/tests/TestAcc.C
new file mode 100644 (file)
index 0000000..b012f93
--- /dev/null
@@ -0,0 +1,244 @@
+#include <TMath.h>
+#include <TGraph.h>
+#include <TMultiGraph.h>
+#include <TMarker.h>
+#include <TLine.h>
+#include <TArc.h>
+#include <TCanvas.h>
+#include <TH2F.h>
+#include <THStack.h>
+
+//_____________________________________________________________________
+void AcceptanceCorrection(Char_t r, UShort_t t, Float_t& oldm, Float_t& newm)
+{
+  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+  
+  //AliFMDRing fmdring(ring);
+  //fmdring.Init();
+  // Float_t   rad       = pars->GetMaxR(r)-pars->GetMinR(r);
+  // Float_t   slen      = pars->GetStripLength(r,t);
+  // Float_t   sblen     = pars->GetBaseStripLength(r,t);
+  const Double_t ic1[] = { 4.9895, 15.3560 };
+  const Double_t ic2[] = { 1.8007, 17.2000 };
+  const Double_t oc1[] = { 4.2231, 26.6638 };
+  const Double_t oc2[] = { 1.8357, 27.9500 };
+
+  Double_t  minR      = (r == 'I' ?  4.5213 : 15.4);
+  Double_t  maxR      = (r == 'I' ? 17.2    : 28.0);
+  Double_t  rad       = maxR - minR;
+  
+  Int_t     nStrips   = (r == 'I' ? 512 : 256);
+  Int_t     nSec      = (r == 'I' ?  20 :  40);
+  Float_t   segment   = rad / nStrips;
+  Float_t   basearc   = 2 * TMath::Pi() / (0.5 * nSec);
+  Float_t   radius    = minR + t * segment;
+  Float_t   baselen   = 0.5 * basearc * radius;
+
+  const Double_t* c1  = (r == 'I' ? ic1 : oc1);
+  const Double_t* c2  = (r == 'I' ? ic2 : oc2);
+
+  // If the radius of the strip is smaller than the radius corresponding 
+  // to the first corner we have a full strip length 
+  Float_t   cr        = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
+  if (radius <= cr) newm = 1;
+  else {
+    // Next, we should find the end-point of the strip - that is, 
+    // the coordinates where the line from c1 to c2 intersects a circle 
+    // with radius given by the strip. 
+    // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
+    Float_t D   = c1[0] * c2[1] - c1[1] * c2[0];
+    Float_t dx  = c2[0] - c1[0];
+    Float_t dy  = c2[1] - c1[1];
+    Float_t dr  = TMath::Sqrt(dx*dx+dy*dy);
+    Float_t det = radius * radius * dr * dr - D*D;
+    
+    if      (det <  0) newm = 1; // No intersection
+    else if (det == 0) newm = 1; // Exactly tangent
+    else {
+      Float_t x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
+      Float_t y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
+      Float_t th  = TMath::ATan2(x, y);
+
+      newm = th / (basearc/2);
+    }
+  }
+
+  Float_t   slope     = (c1[1] - c2[1]) / (c1[0] - c2[0]);
+  Float_t   constant  = (c2[1] * c1[0] - c2[0] * c1[1]) / (c1[0]-c2[0]);
+  
+  Float_t   d         = (TMath::Power(TMath::Abs(radius*slope),2) + 
+                        TMath::Power(radius,2) - TMath::Power(constant,2));
+  
+  // If below corners return 1
+  Float_t arclen = baselen;
+  if (d > 0) {
+    Float_t   x         = ((-TMath::Sqrt(d) - slope * constant) / 
+                          (1+TMath::Power(slope, 2)));
+    Float_t   y         = slope*x + constant;
+    
+    // If x is larger than corner x or y less than corner y, we have full
+    // length strip
+    if(x < c1[0] && y> c1[1]) {
+      //One sector since theta is by definition half-hybrid
+      Float_t   theta     = TMath::ATan2(x,y);
+      arclen              = radius * theta;
+    }
+  }
+  // Calculate the area of a strip with no cut
+  Float_t   basearea1 = 0.5 * baselen * TMath::Power(radius,2);
+  Float_t   basearea2 = 0.5 * baselen * TMath::Power((radius-segment),2);
+  Float_t   basearea  = basearea1 - basearea2;
+  
+  // Calculate the area of a strip with cut
+  Float_t   area1     = 0.5 * arclen * TMath::Power(radius,2);
+  Float_t   area2     = 0.5 * arclen * TMath::Power((radius-segment),2);
+  Float_t   area      = area1 - area2;
+  
+  // Acceptance is ratio 
+  oldm = area/basearea;
+}
+
+void DrawSolution(Char_t r, UShort_t dt=16, UShort_t offT=128)
+{
+  TCanvas* c = new TCanvas(Form("c%c", r), r == 'I' ? 
+                          "Inner" : "Outer", 800, 500);
+  
+  const Double_t ic1[] = { 4.9895, 15.3560 };
+  const Double_t ic2[] = { 1.8007, 17.2000 };
+  const Double_t oc1[] = { 4.2231, 26.6638 };
+  const Double_t oc2[] = { 1.8357, 27.9500 };
+
+  const Double_t* c1   = (r == 'I' ? ic1 : oc1);
+  const Double_t* c2   = (r == 'I' ? ic2 : oc2);
+  Double_t  minR       = (r == 'I' ?  4.5213 : 15.4);
+  Double_t  maxR       = (r == 'I' ? 17.2    : 28.0);
+  Double_t  rad        = maxR - minR;
+  Float_t   cr         = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
+  Double_t  theta      = (r == 'I' ? 18 : 9);
+  TH2*   frame = new TH2F("frame", "Frame", 100, -10, 10, 100, 0, 30);
+  frame->Draw();
+
+  TLine* line = new TLine(c1[0], c1[1], c2[0], c2[1]);
+  line->SetLineColor(kBlue+1);
+  line->Draw();
+
+  UShort_t nT = (r == 'I' ? 512 : 256);
+  for (Int_t t = offT; t < nT; t += dt) {
+    Double_t radius = minR + t * rad / nT;
+
+    TArc*    circle = new TArc(0, 0, radius, 90-theta, 90+theta);
+    circle->SetFillColor(0);
+    circle->SetFillStyle(0);
+    circle->Draw();
+
+    // Now find the intersection 
+    if (radius <= cr) continue;
+
+    // Next, we should find the end-point of the strip - that is, 
+    // the coordinates where the line from c1 to c2 intersects a circle 
+    // with radius given by the strip. 
+    // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
+    Float_t D   = c1[0] * c2[1] - c1[1] * c2[0];
+    Float_t dx  = c2[0] - c1[0];
+    Float_t dy  = c2[1] - c1[1];
+    Float_t dr  = TMath::Sqrt(dx*dx+dy*dy);
+    Float_t det = radius * radius * dr * dr - D*D;
+    
+    if (det <  0) continue;
+    if (det == 0) continue;
+
+
+    Float_t x   = (+D * dy - dx * TMath::Sqrt(det)) / dr / dr;
+    Float_t y   = (-D * dx - dy * TMath::Sqrt(det)) / dr / dr;
+    
+    TMarker* m = new TMarker(x, y, 20);
+    m->SetMarkerColor(kRed+1);
+    m->Draw();
+
+    x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
+    y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
+
+    // Float_t th = TMath::ATan2(x,y);
+    // Printf("theta of strip %3d: %f degrees", 180. / TMath::Pi() * th);
+
+    m = new TMarker(x, y, 20);
+    m->SetMarkerColor(kGreen+1);
+    m->Draw();
+  }
+  c->cd();
+}
+    
+  
+void TestAcc()
+{
+  TCanvas* c =  new TCanvas("c", "C");
+  TGraph*      innerOld = new TGraph(512);
+  TGraph*      outerOld = new TGraph(256);
+  TGraph*      innerNew = new TGraph(512);
+  TGraph*      outerNew = new TGraph(256);
+  innerOld->SetName("innerOld");
+  innerOld->SetName("Inner type, old");
+  innerOld->SetMarkerStyle(21);
+  innerOld->SetMarkerColor(kRed+1);
+  outerOld->SetName("outerOld");
+  outerOld->SetName("Outer type, old");
+  outerOld->SetMarkerStyle(21);
+  outerOld->SetMarkerColor(kBlue+1);
+  innerNew->SetName("innerNew");
+  innerNew->SetName("Inner type, new");
+  innerNew->SetMarkerStyle(20);
+  innerNew->SetMarkerColor(kGreen+1);
+  outerNew->SetName("outerNew");
+  outerNew->SetName("Outer type, new");
+  outerNew->SetMarkerStyle(20);
+  outerNew->SetMarkerColor(kMagenta+1);
+
+  TMultiGraph* all      = new TMultiGraph("all", "Acceptances");
+  all->Add(innerOld);
+  all->Add(outerOld);
+  all->Add(innerNew);
+  all->Add(outerNew);
+
+  TH2* innerCor = new TH2F("innerCor", "Inner correlation", 100,0,1,100,0,1);
+  TH2* outerCor = new TH2F("outerCor", "Outer correlation", 100,0,1,100,0,1);
+  innerCor->SetMarkerStyle(20);
+  outerCor->SetMarkerStyle(20);
+  innerCor->SetMarkerColor(kRed+1);
+  outerCor->SetMarkerColor(kBlue+1);
+  // innerCor->SetMarkerSize(1.2);
+  // outerCor->SetMarkerSize(1.2);
+  THStack* stack = new THStack("corr", "Correlations");
+  stack->Add(innerCor);
+  stack->Add(outerCor);
+
+  for (Int_t i = 0; i < 512; i++) { 
+    Float_t oldAcc, newAcc;
+    
+    AcceptanceCorrection('I', i, oldAcc, newAcc);
+    innerOld->SetPoint(i, i, oldAcc);
+    innerNew->SetPoint(i, i, newAcc);
+    // Printf("Inner strip # %3d:  old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
+
+    innerCor->Fill(oldAcc, newAcc);
+    if (i >= 256) continue;
+    
+    AcceptanceCorrection('O', i, oldAcc, newAcc);
+    outerOld->SetPoint(i, i, oldAcc);
+    outerNew->SetPoint(i, i, newAcc);
+    // Printf("Outer strip # %3d:  old=%8.6f, new=%8.6f", i, oldAcc, newAcc);
+    outerCor->Fill(oldAcc, newAcc);
+  }
+
+  all->Draw("ap");
+
+  DrawSolution('I',4,256);
+  DrawSolution('O',4,128);
+
+  TCanvas* c2 = new TCanvas("cc", "cc");
+  stack->Draw("nostack p");
+  TLine* l = new TLine(0,0,1,1);
+  l->SetLineStyle(2);
+  l->Draw();
+
+  c2->cd();
+}
diff --git a/PWG2/FORWARD/analysis2/tests/TestELossDist.C b/PWG2/FORWARD/analysis2/tests/TestELossDist.C
new file mode 100644 (file)
index 0000000..d48c538
--- /dev/null
@@ -0,0 +1,1426 @@
+#include <iomanip>
+#ifndef __CINT__
+#include <iostream>
+#include <TH1.h>
+#include <TH2.h>
+#include <TRandom.h>
+#include <TCanvas.h>
+#include <TArrayD.h>
+#include <TStyle.h>
+// #include <TFitResult.h>
+#include <TGraphErrors.h>
+#include <TF1.h>
+#include <TMath.h>
+#include <THStack.h>
+#include <TList.h>
+#include <TLatex.h>
+#include <TProfile.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#else
+class TH1;
+class TH2;
+class TArrayD;
+class TRandom;
+class TCanvas;
+class TF1;
+#endif
+
+//___________________________________________________________________
+static Double_t landauGaus1(Double_t* xp, Double_t* pp);
+static Double_t landauGausN(Double_t* xp, Double_t* pp);
+static Double_t landauGausI(Double_t* xp, Double_t* pp);
+
+//====================================================================
+struct Function
+{
+  /**
+   * Enumeration of parameters 
+   * 
+   */
+  enum { 
+    /// Constant 
+    kC, 
+    /// MPV of landau
+    kDelta, 
+    /// Width of landau
+    kXi, 
+    /// Sigma of gaussian
+    kSigma, 
+    /// Number of particles
+    kN, 
+    /// Weights 
+    kA 
+  };
+  /// MP shift
+  static const Double_t fgkMPShift;
+  /// 1 / sqrt(2 * pi)
+  static const Double_t fgkInvRoot2Pi;
+  /// Number of sigmas to integrate over
+  static const Double_t fgkConvNSigma;
+  /// Number of steps in integration 
+  static const Double_t fgkConvNSteps;
+  /** 
+   * Evaluate shifted landay 
+   * @f[                                                                        
+   *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)                      
+   * @f]                                                                        
+   *                                                                            
+   * where @f$ f_{L}@f$ is the ROOT implementation of the Landau                
+   * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for              
+   * @f$\Delta=0,\xi=1@f$.                                                      
+   *                                                                            
+   * @param x      Where to evaluate @f$ f'_{L}@f$                              
+   * @param delta  Most probable value                                          
+   * @param xi     The 'width' of the distribution                              
+   *                                                                            
+   * @return @f$ f'_{L}(x;\Delta,\xi) @f$                                       
+   */
+  static Double_t  Landau(Double_t x, Double_t delta, Double_t xi)
+  {
+    return TMath::Landau(x, delta - xi * fgkMPShift, xi);
+  }
+  /** 
+   * Calculate the value of a Landau convolved with a Gaussian                  
+   *                                                                            
+   * @f[                                                                        
+   * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}                   
+   *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)                 
+   *    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma^2}}                            
+   * @f]                                                                        
+   *                                                                            
+   * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
+   * energy loss, @f$ \xi@f$ the width of the Landau, and @f$\sigma@f$
+   * is the variance of the Gaussian.
+   *                                                                            
+   * Note that this function uses the constants fgkConvNSteps and          
+   * fgkConvNSigma                                                        
+   *                                                                            
+   * References: 
+   *  - <a href="dx.doi.org/10.1016/0168-583X(84)90472-5">
+   *   Nucl.Instrum.Meth.B1:16</a>
+   *  - <a href="dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
+   *  - <a href="root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">
+   *   ROOT implementation</a> 
+   *                                                                            
+   * @param x         where to evaluate @f$ f@f$                                
+   * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$           
+   * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$  
+   * @param sigma     @f$ \sigma@f$ of Gaussian
+   *                                                                            
+   * @return @f$ f@f$ evaluated at @f$ x@f$.                                    
+   */
+  static Double_t  LandauGaus(Double_t x,  Double_t delta, 
+                             Double_t xi, Double_t sigma) 
+  {
+    Double_t deltap = delta - xi * fgkMPShift;
+    Double_t xlow   = x - fgkConvNSigma * sigma;
+    Double_t xhigh  = x + fgkConvNSigma * sigma;
+    Double_t step   = (xhigh - xlow) / fgkConvNSteps;
+    Double_t sum    = 0;
+
+    for (Int_t i = 0; i < fgkConvNSteps / 2; i++) { 
+      Double_t x1 = xlow  + (i - .5) * step;
+      Double_t x2 = xhigh - (i - .5) * step;
+
+      Double_t s1 = 
+       TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma);
+      Double_t s2 = 
+       TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma);
+      sum += s1 + s2;
+    }
+    Double_t ret = step * sum * fgkInvRoot2Pi / sigma;
+    return ret;
+  }
+  /** 
+   * Evaluate                                                                   
+   * @f[                                                                        
+   *    f_i(x;\Delta,\xi,\sigma) = f(x;\Delta_i,\xi_i,\sigma_i)               
+   * @f]                                                                        
+   * corresponding to @f$ i@f$ particles i.e., with the substitutions           
+   * @f{eqnarray*}{                                                             
+   *    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))\\          
+   *    \xi       \rightarrow \xi_i       &=& i \xi\\                           
+   *    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma\\
+   * @f}                                                                        
+   *                                                                            
+   * @param x        Where to evaluate                                          
+   * @param delta    @f$ \Delta@f$                                              
+   * @param xi       @f$ \xi@f$                                                 
+   * @param sigma    @f$ \sigma@f$                                              
+   * @param i        @f$ i @f$                                                  
+   *                                                                            
+   * @return @f$ f_i @f$ evaluated        * 
+   */
+  static Double_t ILandauGaus(Int_t i, Double_t x, Double_t delta, 
+                             Double_t xi, Double_t sigma)
+  {
+    if (i <= 1) return LandauGaus(x, delta, xi, sigma);
+    Double_t di      = i;
+    Double_t deltai  = i * (delta + xi * TMath::Log(di));
+    Double_t xii     = i * xi;
+    Double_t sigmai  = TMath::Sqrt(di) * sigma;
+    
+    if (sigmai < 1e-10) return Landau(x, deltai, xii);
+    Double_t ret =  LandauGaus(x, deltai, xii, sigmai);;
+    // Info("ILandauGaus", "Fi(%f;%d,%f,%f,%f)->%f",
+    //      x, i, deltai, xii, sigmai, ret);
+    return ret;
+  }
+  /** 
+   * Numerically evaluate 
+   * @f[ 
+   *    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
+   * @f] 
+   * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping 
+   * of the parameters is given by 
+   *
+   * - 0: @f$\Delta@f$ 
+   * - 1: @f$\xi@f$ 
+   * - 2: @f$\sigma@f$ 
+   *
+   * This is the partial derivative with respect to the parameter of
+   * the response function corresponding to @f$ i@f$ particles i.e.,
+   * with the substitutions
+   * @f[ 
+   *    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))\\
+   *    \xi       \rightarrow \xi_i       = i \xi\\
+   *    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma\\
+   * @f] 
+   * 
+   * @param x        Where to evaluate 
+   * @param ipar     Parameter number 
+   * @param dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
+   * @param delta    @f$ \Delta@f$ 
+   * @param xi       @f$ \xi@f$ 
+   * @param sigma    @f$ \sigma@f$ 
+   * @param i        @f$ i@f$
+   * 
+   * @return @f$ f_i@f$ evaluated
+   */  
+  static Double_t IdLandauGausdPar(Int_t    i,     Double_t x, 
+                                  UShort_t ipar,  Double_t dPar, 
+                                  Double_t delta, Double_t xi, 
+                                  Double_t sigma)
+  {
+    if (dPar == 0) return 0;
+    Double_t dp      = dPar;
+    Double_t d2      = dPar / 2;
+    Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
+    Double_t xiI     =  i * xi;
+    Double_t si      =  TMath::Sqrt(Double_t(i));
+    Double_t sigmaI  =  si*sigma;
+    Double_t y1      = 0;
+    Double_t y2      = 0;
+    Double_t y3      = 0;
+    Double_t y4      = 0;
+    switch (ipar) {
+    case 0: 
+      y1 = ILandauGaus(i, x, deltaI+i*dp, xiI, sigmaI);
+      y2 = ILandauGaus(i, x, deltaI+i*d2, xiI, sigmaI);
+      y3 = ILandauGaus(i, x, deltaI-i*d2, xiI, sigmaI);
+      y4 = ILandauGaus(i, x, deltaI-i*dp, xiI, sigmaI);
+      break;
+    case 1: 
+      y1 = ILandauGaus(i, x, deltaI, xiI+i*dp, sigmaI);
+      y2 = ILandauGaus(i, x, deltaI, xiI+i*d2, sigmaI);
+      y3 = ILandauGaus(i, x, deltaI, xiI-i*d2, sigmaI);
+      y4 = ILandauGaus(i, x, deltaI, xiI-i*dp, sigmaI);
+      break;
+    case 2: 
+      y1 = ILandauGaus(i, x, deltaI, xiI, sigmaI+si*dp);
+      y2 = ILandauGaus(i, x, deltaI, xiI, sigmaI+si*d2);
+      y3 = ILandauGaus(i, x, deltaI, xiI, sigmaI-si*d2);
+      y4 = ILandauGaus(i, x, deltaI, xiI, sigmaI-si*dp);
+      break;
+    default:
+      return 0;
+    } 
+    
+    Double_t d0  = y1 - y4;
+    Double_t d1  = 2 * (y2 - y3);
+    
+    Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
+    
+    return g;
+  }
+  /** 
+   * Evaluate                                                                   
+   * @f[                                                                        
+   *   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)  
+   * @f]                                                                        
+   *                                                                            
+   * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a               
+   * Landau with a Gaussian (see LandauGaus).  Note that                        
+   * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,                 
+   * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.      
+   *                                                                            
+   * References:                                                                
+   *  - <a href="dx.doi.org/10.1016/0168-583X(84)90472-5">
+   *    Nucl.Instrum.Meth.B1:16</a>
+   *  - <a href="dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
+   * @param x        Where to evaluate @f$ f_N@f$                               
+   * @param delta    @f$ \Delta_1@f$                                            
+   * @param xi       @f$ \xi_1@f$                                               
+   * @param sigma    @f$ \sigma_1@f$                                            
+   * @param n        @f$ N@f$ in the sum above.                                 
+   * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for     
+   *                 @f$ i > 1@f$                                               
+   *                                                                            
+   * @return @f$ f_N(x;\Delta,\xi,\sigma)@f$                                   
+   */  
+  static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
+                             Double_t sigma, Int_t n, Double_t* a)
+  {
+    Double_t res = LandauGaus(x, delta, xi, sigma);
+    for (Int_t i = 2; i <= n; i++) 
+      res += a[i-2] * ILandauGaus(i, x, delta, xi, sigma);
+    return res;
+  }
+  /** 
+   * Calculate the the estimate of number of particle contribtutions
+   * at energy loss @a x 
+   * @f[
+   *   E(\Delta;\Delta_p,\xi,\sigma,\mathbf{a}) = 
+   *   \frac{\sum_{i=1}^N i a_i f_i(x;\Delta_p,\xi,\sigma)}{
+   *     \sum_{i=1}^N a_i f_i(x;\Delta_p,\xi,\sigma)}
+   * @f]
+   * where @f$a_1=1@f$ 
+   * 
+   * @param x      Energy loss @f$\Delta@f$
+   * @param delta  @f$\Delta_{p}@f$
+   * @param xi     @f$\xi@f$ 
+   * @param sigma  @f$\sigma@f$
+   * @param n      Maximum number of particles @f$N@f$ 
+   * @param a      Weights @f$a_i@f$ for @f$i=2,...,N@f$ 
+   * 
+   * @return @f$E(\Delta;\Delta_p,\xi,\sigma,\mathbf{a})@f$
+   */
+  Double_t NEstimate(Double_t x, Double_t delta, Double_t xi, 
+                    Double_t sigma, Int_t n, Double_t* a)
+  {
+    Double_t num = LandauGaus(x, delta, xi, sigma);
+    Double_t den = num;
+    for (Int_t i = 2; i <= n; i++) {
+      Double_t f =  ILandauGaus(i, x, delta, xi, sigma);
+      num        += a[i-2] * f * i;
+      den        += a[i-2] * f;
+    }
+    if (den < 1e-4) return 0;
+    return num / den;
+  }
+  /** 
+   * Calculate the partial derivative of the function 
+   * @f$E(\Delta;\Delta_p,\xi,\sigma,\mathbf{a})@f$ with respect to
+   * one of the parameters @f$\Delta_{p}@f$, @f$\xi@f$, @f$\sigma@f$,
+   * or @f$a_i@f$.  Note for that 
+   * @f{eqnarray*}
+   *   \frac{\partial E}{\partial C}   & \equiv & 0\\
+   *   \frac{\partial E}{\partial a_1} & \equiv & 0\\
+   *   \frac{\partial E}{\partial N}   & \equiv & 0\\
+   * @f{eqnarray*}
+   * 
+   * @param x       Where to evaluate the derivative
+   * @param ipar    Parameter to differentiate relative to 
+   * @param dPar    Variation in parameter to use in calculation
+   * @param delta   @f$\Delta_{p}@f$ 
+   * @param xi      @f$\xi@f$
+   * @param sigma   @f$\sigma@f$ 
+   * @param n       @f$N@f$ 
+   * @param a       @f$\mathbf{a} = (a_2,...,a_N)@f$ 
+   * 
+   * @return @f$\partial E/\partial p|_{x}@f$  
+   */
+  Double_t DEstimatePar(Double_t x, Int_t ipar, Double_t delta, Double_t xi, 
+                       Double_t sigma, UShort_t n, Double_t* a, 
+                       Double_t dPar=0.001) 
+  {
+    Double_t dp    = dPar;
+    Double_t d2    = dPar / 2;
+    Double_t y1    = 0;
+    Double_t y2    = 0;
+    Double_t y3    = 0;
+    Double_t y4    = 0;
+    switch (ipar) { 
+    case Function::kC: 
+    case Function::kN: 
+      return 0;
+    case Function::kDelta:
+      y1 = NEstimate(x, delta+dp, xi, sigma, n, a);
+      y2 = NEstimate(x, delta+d2, xi, sigma, n, a);
+      y3 = NEstimate(x, delta-d2, xi, sigma, n, a);
+      y4 = NEstimate(x, delta-dp, xi, sigma, n, a);
+      break;
+    case Function::kXi:
+      y1 = NEstimate(x, delta, xi+dp, sigma, n, a);
+      y2 = NEstimate(x, delta, xi+d2, sigma, n, a);
+      y3 = NEstimate(x, delta, xi-d2, sigma, n, a);
+      y4 = NEstimate(x, delta, xi-dp, sigma, n, a);
+      break;
+    case Function::kSigma:
+      y1 = NEstimate(x, delta, xi, sigma+dp, n, a);
+      y2 = NEstimate(x, delta, xi, sigma+d2, n, a);
+      y3 = NEstimate(x, delta, xi, sigma-d2, n, a);
+      y4 = NEstimate(x, delta, xi, sigma-dp, n, a);
+      break;
+    default: {
+      Int_t j = ipar-kA;
+      if (j+1 > n) return 0;
+      Double_t aa = a[j];
+      a[j] = aa+dp; y1 = NEstimate(x, delta, xi, sigma, n, a);
+      a[j] = aa+d2; y2 = NEstimate(x, delta, xi, sigma, n, a);
+      a[j] = aa-d2; y3 = NEstimate(x, delta, xi, sigma, n, a);
+      a[j] = aa-dp; y4 = NEstimate(x, delta, xi, sigma, n, a);
+      a[j] = aa;
+    }
+      break;
+    }
+    Double_t d0  = y1 - y4;
+    Double_t d1  = 2 * (y2 - y3);
+    
+    Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;    
+    return g;
+  }
+  /** 
+   * Generate a TF1 object of @f$ f_I@f$                                        
+   *                                                                            
+   * @param n        @f$ n@f$ - the number of particles                         
+   * @param c        Constant                                                   
+   * @param delta    @f$ \Delta@f$                                              
+   * @param xi       @f$ \xi_1@f$                                               
+   * @param sigma    @f$ \sigma_1@f$                                            
+   * @param a        Array of @f$a_i@f$ for @f$i > 1@f$
+   * @param xmin     Least value of range                                       
+   * @param xmax     Largest value of range                                     
+   *                                                                            
+   * @return Newly allocated TF1 object                                         
+   */  
+  static TF1* MakeFunc(Int_t n, Double_t c, Double_t delta, 
+                      Double_t xi, Double_t sigma, Double_t* a, 
+                      Double_t xmin, Double_t xmax)
+  {
+    Int_t nPar = kN + n;
+    TF1* f = new TF1(Form("landGaus%d",n), &landauGausN, xmin, xmax, nPar);
+    f->SetNpx(500);
+    f->SetParName(kC,      "C");
+    f->SetParName(kDelta,  "#Delta_{p}");
+    f->SetParName(kXi,     "#xi");
+    f->SetParName(kSigma,  "#sigma");
+    f->SetParName(kN,      "N");
+    f->SetParameter(kC,     c);
+    f->SetParameter(kDelta, delta);
+    f->SetParameter(kXi,    xi);
+    f->SetParameter(kSigma, sigma);
+    f->FixParameter(kN,     n);
+    f->SetParLimits(kDelta, xmin, 4);
+    f->SetParLimits(kXi,    0.00, 4);
+    f->SetParLimits(kSigma, 0.01, 0.11);
+
+    for (Int_t i = 2; i <= n; i++) { 
+      Int_t j = i-2;
+      f->SetParName(kA+j, Form("a_{%d}", i));
+      f->SetParameter(kA+j, a[j]);
+      f->SetParLimits(kA+j, 0, 1);
+    }
+    return f;
+  }
+  /** 
+   * Make a ROOT TF1 function object corresponding to a single 
+   * component for @f$ i@f$ particles 
+   * 
+   * @param i        Number of particles
+   * @param c        @f$C@f$ 
+   * @param delta    @f$ \Delta@f$                                              
+   * @param xi       @f$ \xi_1@f$                                               
+   * @param sigma    @f$ \sigma_1@f$                                            
+   * @param xmin     Minimum of range of @f$\Delta@f$
+   * @param xmax     Maximum of range of @f$\Delta@f$
+   * 
+   * @return Pointer to newly allocated ROOT TF1 object 
+   */
+  static TF1* MakeIFunc(Int_t i, Double_t c, Double_t delta, 
+                       Double_t xi, Double_t sigma,
+                       Double_t xmin, Double_t xmax)
+  {
+    Int_t nPar = 5;
+    TF1* f = new TF1(Form("landGausI%d",i), &landauGausI, xmin, xmax, nPar);
+    f->SetNpx(100);
+    f->SetParName(kC,      "C");
+    f->SetParName(kDelta,  "#Delta_{p}");
+    f->SetParName(kXi,     "#xi");
+    f->SetParName(kSigma,  "#sigma");
+    f->SetParName(kN,      "N");
+    f->SetParameter(kC,     c);
+    f->SetParameter(kDelta, delta);
+    f->SetParameter(kXi,    xi);
+    f->SetParameter(kSigma, sigma);
+    f->FixParameter(kN,     i);
+    return f;
+  }
+
+  // --- Object code -------------------------------------------------
+  /** 
+   * Constructor 
+   * 
+   * @param n 
+   * @param c 
+   * @param delta 
+   * @param xi 
+   * @param sigma 
+   * @param a 
+   * @param xmin 
+   * @param xmax 
+   * 
+   * @return 
+   */
+  Function(Int_t n, Double_t c, Double_t delta, Double_t xi, Double_t sigma, 
+          Double_t* a, Double_t xmin, Double_t xmax)
+    : fF(MakeFunc(n, c, delta, xi, sigma, a, xmin, xmax))
+  {}
+  Function(TF1* f) : fF(f) {}
+  /** 
+   * Get pointer to ROOT TF1 object 
+   * 
+   * @return Pointer to TF1 object
+   */
+  TF1* GetF1() { return fF; }
+  /** 
+   * Evaluate the function at @a x
+   * @f[ 
+   *  f_N(x;\Delta,\xi,\sigma) = 
+   *     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i)
+   * @f] 
+   * 
+   * @param x Where to evaluate 
+   * 
+   * @return 
+   */
+  Double_t Evaluate(Double_t x) const
+  {
+    return fF->Eval(x);
+  }
+  /** 
+   * Evaluate the function at @a x
+   * @f[ 
+   *  f_N(x;\Delta,\xi,\sigma) = 
+   *     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i)
+   * @f] 
+   * 
+   * This also calculates the error on the point like 
+   * @f[
+   *   \delta^2 f_N = 
+   *    \left(\frac{\partial f_N}{\partial\Delta_p}\right)^2\delta^2\Delta_p+ 
+   *    \left(\frac{\partial f_N}{\partial\xi     }\right)^2\delta^2\xi+ 
+   *    \left(\frac{\partial f_N}{\partial\sigma  }\right)^2\delta^2\sigma+ 
+   *    \sum_{i=2}^N\left(\frac{\partial f_N}{\partial a_i}\right)^2\delta^2a_i
+   * @f]
+   *
+   * @param x Where to evaluate 
+   * @param e On return, contains the error @f$\delta f_N@f$ on the
+   * function value
+   * 
+   * @return @f$ f_N(x;\Delta,\xi,\sigma)@f$ 
+   */
+  Double_t Evaluate(Double_t x, Double_t& e) const
+  {
+    Double_t delta     = GetDelta();
+    Double_t xi        = GetXi();
+    Double_t sigma     = GetSigma();
+    Double_t dFdDelta2 = 0;
+    Double_t dFdXi2    = 0;
+    Double_t dFdSigma2 = 0;
+    e              = 0;
+    for (Int_t i = 1; i <= GetN(); i++) { 
+      Double_t a      = GetA(i);
+      Double_t dDelta = a*IdLandauGausdPar(i, x, 1, 0.001, delta, xi, sigma); 
+      Double_t dXi    = a*IdLandauGausdPar(i, x, 2, 0.001, delta, xi, sigma);
+      Double_t dSigma = a*IdLandauGausdPar(i, x, 3, 0.001, delta, xi, sigma);
+      Double_t dFda   = a*ILandauGaus(i, x, delta, xi, sigma);
+      dFdDelta2 += dDelta * dDelta;
+      dFdXi2    += dXi    * dXi;
+      dFdSigma2 += dSigma * dSigma;
+      e += TMath::Power(dFda * GetEA(i),2);
+    }      
+    Double_t edelta = GetEDelta();
+    Double_t exi    = GetEXi();
+    Double_t esigma = GetESigma();
+    e += (dFdDelta2 * edelta * edelta + 
+         dFdXi2    * exi    * exi    + 
+         dFdSigma2 * esigma * esigma);
+    e = TMath::Sqrt(e);
+    return fF->Eval(x);
+  }
+
+  /** 
+   * Evaluate 
+   * @f[ 
+   *   E(x;\Delta,\xi,\sigma,\mathbf{a}) = 
+   *   \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma)}{
+   *     f_N(x;\Delta,\xi,\sigma)} = 
+   *   \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i)}{
+   *     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i)}
+   * @f] 
+   *
+   * If the denominator is less than @f$10^{-4}@f$, then 0 is returned.
+   * 
+   * @param x    Where to evaluate @f$ E@f$ 
+   * @param maxN Maximum number of weights 
+   * 
+   * @return @f$ E(x;\Delta,\xi,\sigma,\mathbf{a})@f$
+   */  
+  Double_t EstimateNParticles(Double_t x, Short_t maxN=-1)
+  {
+    UShort_t n     = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
+    Double_t delta = GetDelta();
+    Double_t xi    = GetXi();
+    Double_t sigma = GetSigma();
+    return NEstimate(x, delta, xi, sigma, n, GetAs());
+  }
+  /** 
+   * Evaluate 
+   * @f[ 
+   *   E(x;\Delta,\xi,\sigma,\mathbf{a}) = 
+   *   \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma)}{
+   *     f_N(x;\Delta,\xi,\sigma,\mathbf{a})} = 
+   *   \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i)}{
+   *     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i)}
+   * @f] 
+   *
+   * If @f$f_N(x;\Delta,\xi,\sigma,\mathbf{a})<10^{-4}@f$ then 0 is
+   * returned.
+   *
+   * This also calculatues the error @f$\delta E@f$ of the
+   * function value at @f$x@f$: 
+   *
+   * @f[
+   *   \delta^2 E = 
+   *    \left(\frac{\partial E}{\partial\Delta_p}\right)^2\delta^2\Delta_p+ 
+   *    \left(\frac{\partial E}{\partial\xi     }\right)^2\delta^2\xi+ 
+   *    \left(\frac{\partial E}{\partial\sigma  }\right)^2\delta^2\sigma+ 
+   *    \sum_{i=2}^N\left(\frac{\partial E}{\partial a_i}\right)^2\delta^2a_i
+   * @f]
+   * The partial derivatives are evaluated numerically. 
+   * 
+   * @param x    Where to evaluate @f$ E@f$ 
+   * @param e    On return, @f$\delta E|_{x}@f$ 
+   * @param maxN Maximum number of weights 
+   * 
+   * @return @f$ E(x;\Delta,\xi,\sigma,\mathbf{a})@f$
+   */  
+  Double_t EstimateNParticles(Double_t x, Double_t& e, Short_t maxN=-1)
+  {
+    
+    UShort_t  n      = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
+    Double_t  delta  = GetDelta();
+    Double_t  xi     = GetXi();
+    Double_t  sigma  = GetSigma();
+    Double_t* a      = GetAs();
+    Double_t  dDelta = (DEstimatePar(x,Function::kDelta,delta,xi,sigma,n,a,\
+                                    0.01 * GetEDelta()) * GetEDelta());
+    Double_t  dXi    = (DEstimatePar(x,Function::kXi,   delta,xi,sigma,n,a,
+                                    0.01 * GetEXi())    * GetEXi());
+    Double_t  dSigma = (DEstimatePar(x,Function::kSigma,delta,xi,sigma,n,a,
+                                    0.01 * GetESigma()) * GetESigma());
+    e                = dDelta * dDelta + dXi * dXi + dSigma * dSigma;
+    for (Int_t i = 2; i <= n; i++) { 
+      Double_t dAi = (DEstimatePar(x, Function::kA+i-2,delta,xi,sigma,n,a,
+                                  0.01 * GetEA(i)) * GetEA(i));
+      e += dAi * dAi;
+    }
+    e = TMath::Sqrt(e);
+    return NEstimate(x, GetDelta(), GetXi(), GetSigma(), n, GetAs());
+  }
+  /** 
+   * Estimate the number of particles by calculating the weights 
+   * @f[ 
+   *   w_i = a_i f_i(x;\Delta_p,\xi,\sigma) 
+   * @f] 
+   * 
+   * and then draw a random number @f$r@f$ in the range @f$[0,\sum_i^N
+   * w_i]@f$ and return the @f$i@f$ for which @f$w_{i-1} < @r \leq
+   * w_{i}@f$ (notem @f$w_{-1}=0@f$).
+   * 
+   * @param x     Where to evaluate the component functions 
+   * @param maxN  Maximum weight to use 
+   * 
+   * @return Estimate of the number of particles using the procedure
+   * outlined above. 
+   */
+  Double_t RandomEstimateNParticles(Double_t x, Short_t maxN=-1)
+  {
+    UShort_t n     = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
+    TArrayD  p(n);
+    Double_t delta = GetDelta();
+    Double_t xi    = GetXi();
+    Double_t sigma = GetSigma();
+    
+    if (Evaluate(x) < 1e-4) return 0;
+
+    // std::cout << "RandomEstimateNParticles(" << x << "):";
+    for (Int_t i = 1; i <= n; i++) { 
+      Double_t a = GetA(i);
+      Double_t f = ILandauGaus(i, x, delta, xi, sigma);
+      p[i-1]     = a * f;
+      // std::cout << p[i-1] << ",";
+    }
+    Double_t r = gRandom->Uniform(p.GetSum());
+    Double_t l = 0;
+    // std::cout << "s=" << p.GetSum() << ",r=" << r << ",";
+    for (Int_t i = 1; i <= n; i++) { 
+      if (r > l && r <= l + p[i-1]) {
+       // std::cout << "l=" << l << ",l+p=" << l+p[i-1] << "-> " << i 
+       //  << std::endl;
+       return i;
+      }
+      l += p[i-1];
+    }
+    return 0;
+  }
+  /** 
+   * Draw the function 
+   * 
+   * @param option Drawing option 
+   */
+  void Draw(Option_t* option="") { fF->Draw(option); }
+  /** @return @f$C@f$ */
+  Double_t GetC() const { return fF->GetParameter(kC); }
+  /** @return @f$\Delta_{mp}@f$ */
+  Double_t GetDelta() const { return fF->GetParameter(kDelta); }
+  /** @return @f$\xi@f$ */
+  Double_t GetXi() const { return fF->GetParameter(kXi); }
+  /** @return @f$\sigma@f$ */
+  Double_t GetSigma() const { return fF->GetParameter(kSigma); }
+  /** @return @f$N@f$ */
+  UShort_t GetN() const { return fF->GetParameter(kN); }
+  /** @return @f$a_{i}@f$ */
+  Double_t GetA(Int_t i) const { return i == 1 ? 1 : fF->GetParameter(kA+i-2);}
+  /** @return @f$a_2,...,a_N@f$ */
+  Double_t* GetAs() const { return &(fF->GetParameters()[kA]);}
+  /** @return @f$\Delta C@f$ */
+  Double_t GetEC() const { return fF->GetParError(kC); }
+  /** @return @f$\Delta \Delta_{mp}@f$ */
+  Double_t GetEDelta() const { return fF->GetParError(kDelta); }
+  /** @return @f$\Delta \xi@f$ */
+  Double_t GetEXi() const { return fF->GetParError(kXi); }
+  /** @return @f$\Delta \sigma@f$ */
+  Double_t GetESigma() const { return fF->GetParError(kSigma); }
+  /** @return @f$\Delta N@f$ */
+  UShort_t GetEN() const { return 0; }
+  /** @return @f$\Delta a_{i}@f$ */
+  Double_t GetEA(Int_t i) const { return i == 1 ? 0 : fF->GetParError(kA+i-2);}
+  /** @return @f$\Delta a_2,...,\Delta a_N@f$ */
+  Double_t* GetEAs() const { return &(fF->GetParErrors()[kA]);}
+  /** The ROOT object */
+  TF1* fF;
+};
+#ifndef __CINT__
+const Double_t Function::fgkMPShift    = -0.22278298;
+const Double_t Function::fgkInvRoot2Pi = 1. / TMath::Sqrt(2*TMath::Pi());
+const Double_t Function::fgkConvNSigma = 5;
+const Double_t Function::fgkConvNSteps = 100;
+#endif
+  
+//====================================================================
+struct Fitter
+{
+  // --- Object code ------------------------------------------------
+  Fitter(Int_t nMinus) 
+    : fNMinus(nMinus)
+  {
+  }
+  /** 
+   * 
+   * 
+   * @param dist 
+   * 
+   * @return 
+   */
+  TF1* FitFirst(TH1* dist)
+  {
+    fFunctions.Clear();
+    
+    // Get upper edge 
+    Int_t    midBin = dist->GetMaximumBin();
+    Double_t midE   = dist->GetBinLowEdge(midBin);
+    
+    // Get low edge 
+    Int_t    minBin = midBin - fNMinus;
+    Double_t minE  = dist->GetBinCenter(minBin);
+    Double_t maxE  = dist->GetBinCenter(midBin+2*fNMinus);
+    Double_t a[]   = { 0 };
+    TF1* f = Function::MakeFunc(1, 0.5, midE, 0.07, 0.1, a, minE, maxE);
+
+    // Do the fit 
+    dist->Fit(f, "RNQS", "", minE, maxE);
+    fFunctions.AddAtAndExpand(f, 0);
+    return f;
+  }
+  /** 
+   * 
+   * 
+   * @param dist 
+   * @param n 
+   * 
+   * @return 
+   */
+  Function* FitN(TH1* dist, Int_t n)
+  {
+    if (n == 1) return new Function(FitFirst(dist));
+    TF1*        f1 = static_cast<TF1*>(fFunctions.At(0));
+    if (!f1) { 
+      f1 = FitFirst(dist);
+      if (!f1) { 
+       ::Warning("FitN", "First fit missing");
+       return 0;
+      }
+    }
+    
+    Double_t delta1 = f1->GetParameter(Function::kDelta);
+    Double_t xi1    = f1->GetParameter(Function::kXi);
+    Double_t maxEi  = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
+    Double_t minE   = f1->GetXmin();
+    
+    TArrayD a(n-1);
+    for (UShort_t i = 2; i <= n; i++) a.fArray[i-2] = (n==2? 0.05 : 0.000001);
+    
+    TF1* fn = Function::MakeFunc(n, f1->GetParameter(Function::kC), 
+                                delta1, xi1, 
+                                f1->GetParameter(Function::kSigma), 
+                                a.fArray, minE, maxEi);
+    
+    dist->Fit(fn, "RSNQ", "", minE, maxEi);
+    dist->Fit(fn, "RSNQM", "", minE, maxEi);
+    fFunctions.AddAtAndExpand(fn, n-1);
+    
+    return new Function(fn);
+  }
+  Int_t     fNMinus;
+  TObjArray fFunctions;
+};
+
+//====================================================================
+/** 
+ * Utility function 
+ * 
+ * @param xp Pointer to independent variables
+ * @param pp Pointer to parameters 
+ * 
+ * @return Function evaluated at xp[0]
+ */  
+static Double_t landauGaus1(Double_t* xp, Double_t* pp)
+{
+  Double_t x     = xp[0];
+  Double_t c     = pp[Function::kC];
+  Double_t delta = pp[Function::kDelta];
+  Double_t xi    = pp[Function::kXi];
+  Double_t sigma = pp[Function::kSigma];
+    
+  return c * Function::LandauGaus(x, delta, xi, sigma);
+}
+/** 
+ * Utility function 
+ * 
+ * @param xp Pointer to independent variables
+ * @param pp Pointer to parameters 
+ * 
+ * @return Function evaluated at xp[0]
+ */
+static Double_t landauGausN(Double_t* xp, Double_t* pp)
+{
+  Double_t  x     = xp[0];
+  Double_t  c     = pp[Function::kC];
+  Double_t  delta = pp[Function::kDelta];
+  Double_t  xi    = pp[Function::kXi];
+  Double_t  sigma = pp[Function::kSigma];
+  Int_t     n     = Int_t(pp[Function::kN]);
+  Double_t* a     = &pp[Function::kA];
+    
+  return c * Function::NLandauGaus(x, delta, xi, sigma, n, a);
+}
+/** 
+ * Utility function 
+ * 
+ * @param xp Pointer to independent variables
+ * @param pp Pointer to parameters 
+ * 
+ * @return Function evaluated at xp[0]
+ */
+static Double_t landauGausI(Double_t* xp, Double_t* pp)
+{
+  Double_t  x     = xp[0];
+  Double_t  c     = pp[Function::kC];
+  Double_t  delta = pp[Function::kDelta];
+  Double_t  xi    = pp[Function::kXi];
+  Double_t  sigma = pp[Function::kSigma];
+  Int_t     i     = Int_t(pp[Function::kN]);
+    
+  return c * Function::ILandauGaus(i, x, delta, xi, sigma);
+}
+
+//====================================================================
+/**
+ * 
+ * 
+ */
+struct Generator
+{
+  Double_t fDelta;
+  Double_t fXi;
+  Double_t fSigma;
+  Int_t    fMaxN;
+  TH1*     fSingle;
+  TH1*     fSummed;
+  TH1*     fSimple;
+  TH2*     fNSum;
+  TH1*     fN;
+  TArrayD  fRP; // Raw probabilities
+  TArrayD  fP; // Normalized probabilities 
+  TF1*     fF;
+  Int_t    fNObs;
+
+  /** 
+   * 
+   * 
+   * 
+   * @return 
+   */
+  Generator(Double_t delta=.55, Double_t xi=0.06, Double_t sigma=0.02,
+           Int_t maxN=5, const Double_t* a=0)
+    : fDelta(delta), fXi(xi), fSigma(sigma), fMaxN(maxN),
+      fSingle(0), fSummed(0), fSimple(0), fNSum(0), fRP(), fP(fMaxN), fF(0),
+      fNObs(0)
+  {
+    fMaxN = TMath::Min(fMaxN, 5);
+    const Double_t ps[] = { 1, .05, .02, .005, .002 };
+    fRP.Set(5, ps);
+    if (a) { 
+      for (Int_t i = 0; i < maxN; i++) fRP[i] = a[i];
+    }
+
+    Double_t sum = 0; 
+    for (Int_t i = 0; i < fMaxN; i++) sum += fRP[i];
+    for (Int_t i = 0; i < fP.fN; i++) {
+      fP[i] = fRP[i] / sum;
+      if (i != 0) fP[i] += fP[i-1];
+    }
+
+    fSingle = new TH1D("single", "Single signals", 500, 0, 10);
+    fSingle->SetXTitle("Signal");
+    fSingle->SetYTitle("Events");
+    fSingle->SetFillColor(kRed+1);
+    fSingle->SetMarkerColor(kRed+1);
+    fSingle->SetLineColor(kRed+1);
+    fSingle->SetFillStyle(3001);
+    fSingle->SetMarkerStyle(21);
+    fSingle->SetDirectory(0);
+    fSingle->Sumw2();
+
+    fSummed = static_cast<TH1D*>(fSingle->Clone("summed"));
+    fSummed->SetTitle("Summed signals");
+    fSummed->SetFillColor(kBlue+1);
+    fSummed->SetMarkerColor(kBlue+1);
+    fSummed->SetLineColor(kBlue+1);
+    fSummed->SetDirectory(0);
+
+    fSimple = static_cast<TH1D*>(fSingle->Clone("summed"));
+    fSimple->SetTitle("Summed signals");
+    fSimple->SetFillColor(kGreen+1);
+    fSimple->SetMarkerColor(kGreen+1);
+    fSimple->SetLineColor(kGreen+1);
+    fSimple->SetDirectory(0);
+
+    fNSum = new TH2D("single", "Observation signals", 
+                    fMaxN, .5, fMaxN+.5, 100, 0, 10);
+    fNSum->SetMarkerColor(kMagenta+1);
+    fNSum->SetLineColor(kMagenta+1);
+    fNSum->SetFillColor(kMagenta+1);
+    fNSum->SetYTitle("Signal");
+    fNSum->SetXTitle("N");
+    fNSum->SetZTitle("Events");
+    fNSum->SetDirectory(0);
+    fNSum->Sumw2();
+
+    fN = new TH1D("n", "Number of particles", fMaxN, .5, fMaxN+.5);
+    fN->SetXTitle("N");
+    fN->SetYTitle("Events");
+    fN->SetMarkerColor(kMagenta+1);
+    fN->SetLineColor(kMagenta+1);
+    fN->SetFillColor(kMagenta+1);
+    fN->SetMarkerStyle(22);
+    fN->SetFillStyle(3001);
+    fN->SetDirectory(0);
+
+    std::cout << "Probabilities:" << std::setprecision(4) 
+             << "\n  ";
+    for (Int_t i = 0; i < fRP.fN; i++) 
+      std::cout << std::setw(7) << fRP[i] << " ";
+    std::cout << " = " << fRP.GetSum() << "\n  ";
+    for (Int_t i = 0; i < fP.fN; i++) 
+      std::cout << std::setw(7) << fP[i] << " ";
+    std::cout << std::endl;
+
+    Double_t aa[] = { 0 };
+    fF = Function::MakeFunc(1, 1, fDelta, fXi, fSigma, aa,  0, 10);
+  }
+  /** 
+   * 
+   * 
+   */
+  void Clear()
+  {
+    fSingle->Reset();
+    fSummed->Reset();
+    fNSum->Reset();
+    fN->Reset();
+  }
+  /** 
+   * 
+   * 
+   * 
+   * @return 
+   */    
+  Int_t GenerateN()
+  {
+    Double_t rndm = gRandom->Uniform();
+    Int_t    ret  = 0;
+    if      (                 rndm < fP[0]) ret = 1;
+    else if (rndm >= fP[0] && rndm < fP[1]) ret = 2;
+    else if (rndm >= fP[1] && rndm < fP[2]) ret = 3;
+    else if (rndm >= fP[2] && rndm < fP[3]) ret = 4;
+    else if (rndm >= fP[3] && rndm < fP[4]) ret = 5;
+    // Printf("%f -> %d", rndm, ret);
+    fN->Fill(ret);
+    return ret;
+  }
+  /** 
+   * 
+   * 
+   * @param mpv 
+   * @param xi 
+   * @param sigma 
+   * 
+   * @return 
+   */
+  Double_t Generate1Signal()
+  {
+    Double_t ret = 0;
+    if (fF) ret = fF->GetRandom();
+    else { 
+      Double_t rmpv = gRandom->Gaus(fDelta, fSigma);
+      ret           = gRandom->Landau(rmpv, fXi);
+    }
+    fSingle->Fill(ret);
+    fSimple->Fill(gRandom->Landau(fDelta - fXi * Function::fgkMPShift, fXi));
+    return ret;
+  }
+  /** 
+   * 
+   * 
+   * @param n 
+   * @param mpv 
+   * @param xi 
+   * @param sigma 
+   * 
+   * @return 
+   */  
+  Double_t GenerateNSignal(Int_t n)
+  {
+    Double_t ret = 0;
+    for (Int_t i = 0; i < n; i++) ret += Generate1Signal();
+    fSummed->Fill(ret);
+    fNSum->Fill(n, ret);
+    return ret;
+  }
+  /** 
+   * 
+   * 
+   * @param nObs 
+   * @param reset 
+   */  
+  void Generate(Int_t nObs=1000, Bool_t reset=true)
+  {
+    if (reset) Clear();
+    fNObs = nObs;
+
+    for (Int_t i = 0; i < nObs; i++) {
+      // if (((i+1) % (nObs/10)) == 0) 
+      //       std::cout << "Event " << std::setw(6) << i << std::endl;
+      Int_t n = GenerateN();
+      GenerateNSignal(n);
+    }
+    Double_t m = fSingle->GetMaximum();
+    fSingle->Scale(1. / m);
+    m = fSummed->GetMaximum();
+    fSummed->Scale(1. / m);
+    m = fSimple->GetMaximum();
+    fSimple->Scale(1. / m);
+
+    std::cout << "Resulting probabilities:\n  ";
+    for (Int_t i = 1; i <= fN->GetNbinsX(); i++) 
+      std::cout << std::setprecision(4) << std::setw(7) 
+               << fN->GetBinContent(i) << " ";
+    std::cout << std::endl;
+  }
+  /** 
+   * Pretty print of parameters. 
+   * 
+   * @param input 
+   * @param f 
+   * @param i 
+   */
+  void PrintParameter(Double_t input, TF1* f, Int_t i)
+  {
+    Double_t val = 0, err = 0;
+    if (i < f->GetNpar()) {
+      val = f->GetParameter(i);
+      err = f->GetParError(i);
+    }
+    std::cout << "  " << std::setw(16) << f->GetParName(i) << ": "
+             << std::fixed 
+             << std::setprecision(4) << std::setw(6) << input << " -> " 
+             << std::setprecision(4) << std::setw(6) << val   << " +/- "
+             << std::setprecision(5) << std::setw(7) << err   << " [delta:" 
+             << std::setw(3) << int(100 * (input-val) / input) << "%,err:" 
+             << std::setw(3) << int(100 * err / val) << "%]"
+             << std::endl;
+  }
+  /** 
+   * Print the result 
+   * 
+   * @param f 
+   */
+  void PrintResult(TF1* f) 
+  {
+    Double_t chi2 = f->GetChisquare();
+    Int_t    ndf  = f->GetNDF();
+    std::cout << "Chi^2/NDF = " << chi2 << "/" << ndf << "=" 
+             << chi2/ndf << std::endl;
+    PrintParameter(1,      f, Function::kC);
+    PrintParameter(fDelta, f, Function::kDelta);
+    PrintParameter(fXi,    f, Function::kXi);
+    PrintParameter(fSigma, f, Function::kSigma);
+    PrintParameter(fMaxN,  f, Function::kN);
+    for (Int_t i = 0; i < fMaxN-1; i++) 
+      PrintParameter(fRP[i+1], f, Function::kA+i);
+  }
+  /** 
+   * 
+   * 
+   * @param x 
+   * @param y 
+   * @param intput 
+   * @param f 
+   * @param i 
+   */
+  void DrawParameter(Double_t x, Double_t y, Double_t input, TF1* f, Int_t i)
+  {
+    Double_t val = 0, err = 0;
+    if (i < f->GetNpar()) {
+      val = f->GetParameter(i);
+      err = f->GetParError(i);
+    }
+
+    TLatex* ltx = new TLatex(x, y, f->GetParName(i));
+    ltx->SetTextSize(.04);
+    ltx->SetTextFont(132);
+    ltx->SetTextAlign(11);
+    ltx->SetNDC();
+    ltx->Draw();
+    
+    ltx->DrawLatex(x+.05,  y, Form("%5.3f", input));
+    ltx->DrawLatex(x+.14, y, "#rightarrow");
+    ltx->SetTextAlign(21);
+    ltx->DrawLatex(x+.3,  y, Form("%5.3f #pm %6.4f", val, err));
+    ltx->SetTextAlign(11);
+    ltx->DrawLatex(x+.48,  y, Form("[#Delta=%3d%%, #delta=%3d%%]",
+                                 int(100 * (input-val)/input), 
+                                 int(100*err/val)));
+  }
+  /** 
+   * 
+   * 
+   * @param x 
+   * @param y 
+   * @param f 
+   */
+  void DrawResult(Double_t x, Double_t y, TF1* f)
+  {
+    Double_t chi2 = f->GetChisquare();
+    Int_t    ndf  = f->GetNDF();
+    TLatex* ltx = new TLatex(x, y, Form("#chi^{2}/NDF=%5.2f/%3d=%6.3f", 
+                                       chi2, ndf, chi2/ndf));
+    ltx->SetTextSize(.04);
+    ltx->SetNDC();
+    ltx->SetTextFont(132);
+    ltx->Draw();
+
+    x += .05;
+    y -= .035; DrawParameter(x, y, 1,      f, Function::kC);
+    y -= .035; DrawParameter(x, y, fDelta, f, Function::kDelta);
+    y -= .035; DrawParameter(x, y, fXi,    f, Function::kXi);
+    y -= .035; DrawParameter(x, y, fSigma, f, Function::kSigma);
+    y -= .035; DrawParameter(x, y, fMaxN,  f, Function::kN);
+    for (Int_t i = 0; i < fMaxN-1; i++) {
+      y -= .035; DrawParameter(x, y, fRP[i+1], f, Function::kA+i);
+    }
+  }
+  /** 
+   * 
+   * 
+   */
+  void Draw(const char* type="png")
+  {
+    gStyle->SetOptTitle(0);
+    gStyle->SetOptStat(0);
+    gStyle->SetPalette(1);
+
+    TCanvas* c = new TCanvas("c", "Generator distributions", 1200, 1000);
+    c->SetFillColor(0);
+    c->SetBorderSize(0);
+    c->SetBorderMode(0);
+    c->SetTopMargin(0.01);
+    c->SetRightMargin(0.01);
+
+    c->Divide(2, 2);
+
+    // --- First pad:  Data (multi, single, landau) + Fit -------------
+    TVirtualPad* p = c->cd(1);
+    p->SetLogy();
+    p->SetBorderSize(0);
+    p->SetBorderMode(0);
+    p->SetTopMargin(0.01);
+    p->SetRightMargin(0.01);
+    fSummed->Draw();
+    fSingle->Draw("same");
+    fSimple->Draw("same");
+
+    if (fF) { 
+      Double_t max = fF->GetMaximum();
+      Double_t a[] = { 0 };
+      TF1* fcopy = Function::MakeFunc(1, 
+                                     fF->GetParameter(Function::kC) /max, 
+                                     fF->GetParameter(Function::kDelta), 
+                                     fF->GetParameter(Function::kXi), 
+                                     fF->GetParameter(Function::kSigma), 
+                                     a, 0, 10);
+      fcopy->SetLineColor(2);
+      fcopy->SetLineStyle(2);
+      fcopy->Draw("same");
+    }
+
+    Fitter*   fitter = new Fitter(4);
+    Function* fit    = fitter->FitN(fSummed, 5);
+    TF1*      f      = fit->GetF1();
+    f->Draw("same");
+    f->SetLineColor(kBlack);
+    TF1* fst = static_cast<TF1*>(fitter->fFunctions.At(0));
+    fst->SetLineWidth(3);
+    fst->SetLineStyle(2);
+    fst->SetLineColor(kMagenta);
+    fst->Draw("same");
+    PrintResult(f);
+    DrawResult(.25, .95, f);
+    TGraphErrors* gr = new TGraphErrors(fSummed->GetNbinsX());
+    gr->SetName("fitError");
+    gr->SetTitle("Fit w/errors");
+    gr->SetFillColor(kYellow+1);
+    gr->SetFillStyle(3001);
+    for (Int_t j = 1; j <= fSummed->GetNbinsX(); j++) { 
+      Double_t x = fSummed->GetBinCenter(j);
+      Double_t e = 0;
+      Double_t y = fit->Evaluate(x, e);
+      gr->SetPoint(j-1, x, y);
+      gr->SetPointError(j-1, 0, e);
+    }
+    gr->Draw("same l e4");
+
+    TLegend* l = new TLegend(.15, .11, .5, .4);
+    l->SetBorderSize(0);
+    l->SetFillColor(0);
+    l->SetFillStyle(0);
+    l->SetTextFont(132);
+    l->AddEntry(fSimple, "L: Single-L", "lp");
+    l->AddEntry(fSingle, "LG: Single-L#otimesG", "lp");
+    l->AddEntry(fSummed,"NLG: Multi-L#otimesG", "lp");
+    l->AddEntry(gr, "f_{N}: Fit to NLG");
+    l->Draw();
+
+    // --- Second pad:  Data (1, 2, 3, ...) --------------------------
+    p = c->cd(2);
+    p->SetLogy();
+    p->SetBorderSize(0);
+    p->SetBorderMode(0);
+    p->SetTopMargin(0.01);
+    p->SetRightMargin(0.01);
+    p->SetGridx();
+    THStack* stack = new THStack(fNSum, "y");
+    TIter next(stack->GetHists());
+    TH1* h = 0;
+    Int_t i = 2;
+    Double_t max = 0;
+    TLegend* l2 = new TLegend(.6, .5, .95, .95);
+    l2->SetBorderSize(0);
+    l2->SetFillStyle(0);
+    l2->SetFillColor(0);
+    l2->SetTextFont(132);
+    while ((h = static_cast<TH1*>(next()))) { 
+      if (i == 2) max = h->GetMaximum();
+      h->SetLineColor(i); 
+      h->SetFillColor(i); 
+      h->SetFillStyle(3001);
+      h->Scale(1/max);
+      l2->AddEntry(h, Form("%d particles", i - 1), "f");
+      i++;
+    }
+    stack->Draw("hist nostack");
+    stack->GetHistogram()->SetXTitle("Signal");
+    stack->GetHistogram()->SetYTitle("Events");
+    i = 2;
+    for (Int_t j = 1; j <= 5; j++) {
+      TF1* fi = Function::MakeIFunc(j, fRP[j-1], fDelta, fXi, fSigma, 0, 10);
+      if (j == 1) max = fi->GetMaximum();
+      fi->SetParameter(0, fRP[j-1]/max);
+      fi->SetLineColor(i);
+      fi->SetLineWidth(3);
+      fi->SetLineStyle(2);
+      fi->Draw("same");
+
+      TF1* fj = Function::MakeIFunc(j, fit->GetC() * fit->GetA(j), 
+                                   fit->GetDelta(), fit->GetXi(), 
+                                   fit->GetSigma(), 0, 10);
+      fj->SetLineColor(i);
+      fj->SetLineWidth(2);
+      fj->Draw("same");
+
+      std::cout << "Component " << j << " scale=" << fi->GetParameter(0) 
+               << " (" << fRP[j-1] << ")" << " " << fj->GetParameter(0) 
+               << " (" << fit->GetC() << "*" << fit->GetA(j) << ")" 
+               << std::endl;
+      i++;
+    }
+    TLegendEntry* ee = l2->AddEntry("", "Input f_{i}", "l");
+    ee->SetLineStyle(2);
+    ee->SetLineWidth(3);
+    l2->AddEntry("", "Fit f_{i}", "l");
+    l2->Draw();
+    // fNSum->Draw("lego2");
+
+    // --- Third pad:  Mean multiplicity ------------------------------
+    TH1* nhist1 = new TH1F("nhist1", "NHist", 5, 0.5, 5.5);
+    nhist1->SetFillColor(kCyan+1);
+    nhist1->SetMarkerColor(kCyan+1);
+    nhist1->SetLineColor(kCyan+1);
+    nhist1->SetFillStyle(3002);
+    TH1* nhist2 = new TH1F("nhist2", "NHist", 5, 0.5, 5.5);
+    nhist2->SetFillColor(kYellow+1);
+    nhist2->SetMarkerColor(kYellow+1);
+    nhist2->SetLineColor(kYellow+1);
+    nhist2->SetFillStyle(3002);
+
+    TH2* resp = new TH2F("resp", "Reponse", 100, 0, 10, 6, -.5, 5.5);
+    resp->SetXTitle("Signal");
+    resp->SetYTitle("# particles");
+    for (Int_t j = 0; j < fNObs; j++) { 
+      if (((j+1) % (fNObs / 10)) == 0) 
+       std::cout << "Event # " << j+1 << std::endl;
+      Double_t x = fSummed->GetRandom();
+      Double_t y = fit->RandomEstimateNParticles(x);
+      nhist1->Fill(y);
+      resp->Fill(x, y);
+    }
+    TGraphErrors* graph = new TGraphErrors(fSummed->GetNbinsX());
+    graph->SetName("evalWeighted");
+    graph->SetTitle("Evaluated function");
+    graph->SetLineColor(kBlack);
+    graph->SetFillColor(kYellow+1);
+    graph->SetFillStyle(3001);
+    for (Int_t j = 1; j <= fSummed->GetNbinsX(); j++) { 
+      Double_t x = fSummed->GetBinCenter(j);
+      Double_t e = 0;
+      Double_t y = fit->EstimateNParticles(x, e);
+      nhist2->Fill(y, fSummed->GetBinContent(j));
+      graph->SetPoint(j-1, x, y);
+      graph->SetPointError(j-1, 0, e);
+    }
+    nhist1->Scale(1. / nhist1->GetMaximum());
+    nhist2->Scale(1. / nhist2->GetMaximum());
+    fN->Scale(1. / fN->GetMaximum());
+
+    p = c->cd(3);
+    p->SetLogy();
+    p->SetBorderSize(0);
+    p->SetBorderMode(0);
+    p->SetTopMargin(0.01);
+    p->SetRightMargin(0.01);
+    fN->Draw();
+    nhist1->Draw("same");
+    nhist2->Draw("same");
+
+    TLegend* l3 = new TLegend(.3, .7, .95, .95);
+    l3->SetBorderSize(0);
+    l3->SetFillStyle(0);
+    l3->SetFillColor(0);
+    l3->SetTextFont(132);
+    l3->AddEntry(fN,  
+                Form("Input n distribution: #bar{m}=%5.3f, s_{m}=%5.3f",
+                     fN->GetMean(), fN->GetRMS()), "f");
+    l3->AddEntry(nhist1, 
+                Form("Random draw from fit: #bar{m}=%5.3f, s_{m}=%5.3f", 
+                     nhist1->GetMean(), nhist1->GetRMS()), "f");
+    l3->AddEntry(nhist2, 
+                Form("Weighted evaluation of fit: #bar{m}=%5.3f, s_{m}=%5.3f", 
+                     nhist2->GetMean(), nhist2->GetRMS()), "f");
+    l3->Draw();
+    
+    // --- Fourth pad: Reponse function ------------------------------
+    p = c->cd(4);
+    // p->SetLogy();
+    p->SetBorderSize(0);
+    p->SetBorderMode(0);
+    p->SetTopMargin(0.01);
+    p->SetRightMargin(0.01);
+    p->SetGridx();
+    TProfile* prof1 = fNSum->ProfileY();
+    prof1->SetLineColor(kMagenta+1);
+    prof1->Draw();
+    TProfile* prof2 = resp->ProfileX();
+    prof2->SetLineColor(kCyan+2);
+    prof2->Draw("same");
+    graph->SetLineWidth(2);
+    graph->Draw("same LP E4");
+
+    TLegend* l4 = new TLegend(.2, .11, .8, .4);
+    l4->SetBorderSize(0);
+    l4->SetFillStyle(0);
+    l4->SetFillColor(0);
+    l4->SetTextFont(132);
+    l4->AddEntry(prof1, "Input distribution of N particles", "l");
+    l4->AddEntry(prof2, "Random estimate of N particles", "l");
+    l4->AddEntry(graph, "E: #sum_{i}^{N}i a_{i}f_{i}/"
+                "#sum_{i}^{N}a_{i}f_{i} evaluated over NLG", "lf");
+    l4->Draw();
+
+    c->cd();
+    
+    if (type && type[0] != '\0') 
+      c->Print(Form("TestELossDist.%s", type));
+  }
+};
+
+
+void
+TestELossDist(const char* type="png")
+{
+  gRandom->SetSeed(12345);
+  const Double_t a1[] = { 1, .05, .02, .005, .002 }; // 7TeV pp
+  const Double_t a2[] = { 1, .1, .05, .01, .005 };
+  const Double_t a3[] = { 1, .2, .10, .05, .01 };
+
+  Generator* g = new Generator(0.55, 0.06, 0.06, 5, a1);
+  g->Generate(100000);
+  g->Draw(type);
+
+  
+}
+
+
diff --git a/PWG2/FORWARD/analysis2/tests/TestFitELoss.C b/PWG2/FORWARD/analysis2/tests/TestFitELoss.C
new file mode 100644 (file)
index 0000000..ae5a975
--- /dev/null
@@ -0,0 +1,305 @@
+/**
+ * Test script to fit the energy loss spectra 
+ * 
+ * @deprecated
+ * This is a simple test script 
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+#ifndef __CINT__
+# include "AliForwardUtil.h"
+# include <TFile.h>
+# include <TList.h>
+# include <TH1.h>
+# include <TF1.h>
+# include <TFitResult.h>
+# include <TMath.h>
+# include <TStyle.h>
+# include <TArrow.h>
+# include <TCanvas.h>
+#else
+class TCanvas;
+class TFile;
+class TH1;
+class TList;
+class TF1;
+#endif
+
+//____________________________________________________________________
+/** 
+ * 
+ * 
+ * @param ef 
+ * @param d 
+ * @param r 
+ * @param etabin 
+ * 
+ * @return 
+ *
+ * @deprecated
+ * This is a simple test script 
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+TH1* GetEDist(TList* ef, UShort_t d, Char_t r, UShort_t etabin)
+{
+  TList* dL = static_cast<TList*>(ef->FindObject(Form("FMD%d%c",d,r)));
+  if (!dL) {
+    Error("GetEDist", "Couldn't find list FMD%d%c",d,r);
+    ef->ls();
+    return 0;
+  }
+  if (etabin > 999) {
+    TH1* hist = static_cast<TH1*>(dL->FindObject(Form("FMD%d%c_edist",d,r)));
+    if (hist) {
+      Error("GetEDist", "Couldn't find EDists histogram for FMD%d%c",d,r);
+      return 0;
+    }
+  }
+      
+  TList* edL = static_cast<TList*>(dL->FindObject("EDists"));
+  if (!edL) {
+    Error("GetEDist", "Couldn't find list EDists for FMD%d%c",d,r);
+    return 0;
+  }
+  
+  TH1* hist = static_cast<TH1*>(edL->FindObject(Form("FMD%d%c_etabin%03d", 
+                                                    d, r, etabin)));
+  if (!hist) {
+    Error("GetEDist", "Couldn't find histogra FMD%d%c_etabin%03d",
+         d,r, etabin);
+    return 0;
+  }
+  
+  return hist;
+}
+
+//____________________________________________________________________
+/** 
+ * 
+ * 
+ * @param ef 
+ * @param d 
+ * @param r 
+ * @param eta 
+ * 
+ * @return  
+ *
+ * @deprecated
+ * This is a simple test script 
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+*/
+TH1* GetEDist(TList* ef, UShort_t d, Char_t r, Float_t eta)
+{
+  if (!ef) { 
+    Error("GetEDist", "EF not set");
+    return 0;
+  }
+  TAxis* etaAxis = static_cast<TAxis*>(ef->FindObject("etaAxis"));
+  if (!etaAxis) { 
+    Error("GetEDist", "Couldn't find eta axis in list");
+    return 0;
+  }
+
+  UShort_t bin = etaAxis->FindBin(eta);
+  if (bin <= 0 || bin > etaAxis->GetNbins()) { 
+    Error("GetEDist", "eta=%f out of range [%f,%f] - getting ring histo", 
+         eta, etaAxis->GetXmin(), etaAxis->GetXmax());
+    return GetEDist(ef, d, r, UShort_t(1000));
+  }
+
+  return GetEDist(ef, d, r, bin);
+}
+
+//____________________________________________________________________
+/** 
+ * 
+ * 
+ * @param file 
+ * 
+ * @return 
+ *
+ * @deprecated
+ * This is a simple test script 
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+TList* GetEF(TFile* file) 
+{
+  TList* forward = static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
+  if (!forward) {
+    Error("GetEF", "Failed to get list PWG2forwardDnDeta/Forward from file");
+    return 0;
+  }
+  TList* ef = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
+  if (!ef) {
+    Error("GetEF", "Failed to get energy fitter list");
+    return 0;
+  }
+  
+  return ef;
+}
+
+//____________________________________________________________________
+TList* ef = 0;
+
+//____________________________________________________________________
+/** 
+ * 
+ * 
+ * 
+ * @return 
+ *
+ * @deprecated
+ * This is a simple test script 
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+TList*  CheckEF()
+{
+  if (ef) return ef;
+  
+  TFile* file = TFile::Open("AnalysisResults.root", "READ");
+  if (!file) { 
+    Error("Fit1", "Failed to open file");
+    return 0;
+  }
+  return GetEF(file);
+}
+
+//____________________________________________________________________
+TCanvas* c = 0;
+
+//____________________________________________________________________
+/** 
+ * 
+ * 
+ * 
+ * @return 
+ *
+ * @deprecated
+ * This is a simple test script 
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+TCanvas* CheckC()
+{
+  if (c) return c;
+
+  gStyle->SetOptFit(1111);
+  gStyle->SetCanvasColor(0);
+  gStyle->SetCanvasBorderSize(0);
+  gStyle->SetCanvasBorderMode(0);
+  gStyle->SetCanvasDefW(800);
+  gStyle->SetCanvasDefH(800);
+  gStyle->SetPadTopMargin(0.05);
+  gStyle->SetPadRightMargin(0.05);
+  gStyle->SetTitleBorderSize(0);
+  gStyle->SetTitleFillColor(0);
+  gStyle->SetTitleStyle(0);
+  gStyle->SetStatBorderSize(1);
+  gStyle->SetStatColor(0);
+  gStyle->SetStatStyle(0);
+  gStyle->SetStatX(0.95);
+  gStyle->SetStatY(0.95);
+  gStyle->SetStatW(0.15);
+  gStyle->SetStatH(0.15);
+
+  c = new TCanvas("c", "c");
+  c->SetLogy();
+
+  return c;
+}
+
+//____________________________________________________________________
+/** 
+ * 
+ * 
+ * @param f 
+ *
+ * @deprecated
+ * This is a simple test script 
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void PrintFit(TF1* f)
+{
+  Int_t    nu     = f->GetNDF();
+  Double_t chi2   = f->GetChisquare();
+  Double_t chi2nu = (nu > 0 ? chi2/nu : 0);
+  Printf("%s: chi^2/nu=%f/%d=%f [%f,%f]", 
+        f->GetName(), chi2, nu, chi2nu,
+        f->GetXmin(), f->GetXmax());
+  for (Int_t i = 0; i < f->GetNpar(); i++) { 
+    Double_t v = f->GetParameter(i);
+    Double_t e = f->GetParError(i);
+    Double_t r = 100*(v == 0 ? 1 : e / v);
+    Printf("%32s = %14.7f +/- %14.7f (%5.1f%%)",f->GetParName(i),v,e,r);
+  }
+}
+  
+//____________________________________________________________________
+/** 
+ * 
+ * 
+ * @param n 
+ * @param d 
+ * @param r 
+ * @param eta 
+ *
+ * @deprecated
+ * This is a simple test script 
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void TestFitELoss(Int_t n, UShort_t d, Char_t r, Float_t eta)
+{
+  TList* ef1 = CheckEF();
+  TCanvas* c1 = CheckC();
+  if (!ef1) return;
+  if (!c1)  return;
+
+  TH1* dist = GetEDist(ef1, d, r, eta);
+  if (!dist) return;
+
+  AliForwardUtil::ELossFitter f(0.4, 10, 4);
+  TF1* landau1 = new TF1(*f.Fit1Particle(dist, 0));
+  landau1->SetName("Landau1");
+  PrintFit(landau1);
+  TF1* landaun = new TF1(*f.FitNParticle(dist, n, 0)); 
+  landau1->SetName(Form("Landau%d", n));
+  PrintFit(landaun);
+  landau1->SetRange(0,10);
+  landaun->SetRange(0,10);
+  landau1->SetLineWidth(4);
+  landaun->SetLineWidth(4);
+  landau1->SetLineColor(kBlack);
+  landaun->SetLineColor(kBlack);
+
+  dist->GetListOfFunctions()->Clear();
+  dist->GetListOfFunctions()->Add(landau1);
+  dist->GetListOfFunctions()->Add(landaun);
+  dist->GetListOfFunctions()->ls();
+  dist->Draw();
+  landau1->Draw("same");
+  landaun->Draw("same");
+
+  Double_t mp = landaun->GetParameter(1);
+  Double_t xi = landaun->GetParameter(2);
+  for (Int_t i  = 1; i <= n; i++) { 
+    Double_t x  = i * (mp + xi * TMath::Log(i));
+    Double_t y  = landaun->Eval(x);
+    Double_t y1 = y < 0.05 ? 1 : 0.01;
+    TArrow* a = new TArrow(x,y1,x,y,0.03,"|>");
+    Info("FitSteer", "Delta_{p,%d}=%f", i, x);
+    a->SetLineWidth(2);
+    a->SetAngle(30);
+    a->Draw();
+  }
+
+  c1->cd();
+}
+//
+// EOF
+//
diff --git a/PWG2/FORWARD/analysis2/tests/TestMakeELossFits.C b/PWG2/FORWARD/analysis2/tests/TestMakeELossFits.C
new file mode 100644 (file)
index 0000000..16f47f5
--- /dev/null
@@ -0,0 +1,232 @@
+#ifndef __CINT__
+#include "AliForwardCorrectionManager.h"
+#include "AliFMDCorrELossFit.h"
+#include <TAxis.h>
+#include <TFile.h>
+#include <TList.h>
+#include <TError.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TClonesArray.h>
+#else
+class TAxis;
+class AliFMDCirrELossFit;
+class TH1;
+
+#endif
+
+/**
+ * Class to make correction object and save to file 
+ *
+ * Run like                                    
+ *
+ * @verbatim 
+ * gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+ * gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/Compile.C");
+ * Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/MakeELossFit.C"); 
+ * MakeELossFit mef(sys, cms, field, mc, "AnalysisResults.root"); 
+ * mef.Run();
+ * @endverbatim 
+ * where @e sys, the collision system, is one of 
+ * - 1: pp
+ * - 2: PbPb
+ * 
+ * @e cms is the center of mass energy per nuclean in GeV (e.g., 2760
+ * for first PbPb data), @e field is (signed) L3 magnetic in kG, and 
+ * @e mc is wether this correction applies to MC data or not. 
+ * 
+ * The class generates a file like 
+ * @verbatim
+ * elossfits<sys>_<cms>GeV_<field>kG_<realmc>.root
+ * @endverbatim 
+ * in the working directory. This file can be moved to 
+ * @verbatim
+ * $(ALICE_ROOT)/PWG2/FORWARD/corrections/ELossFits
+ * @endverbatim 
+ * and stored in SVN for later use. 
+ *
+ * @depcrecated 
+ * The class AliFMDELossFitter automatically generates the
+ * AliFMDCorrELossFit object.
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+
+class MakeELossFit 
+{
+protected: 
+public:
+  TList* fFitter;
+  TAxis* fAxis;
+  TClonesArray fFits;
+  UShort_t     fSys;
+  UShort_t     fCMS;
+  Short_t      fField;
+  Bool_t       fMC;
+
+  //__________________________________________________________________
+  MakeELossFit(UShort_t    sys, 
+              UShort_t    cms, 
+              Short_t     field, 
+              Bool_t      mc, 
+              const char* filename="forward_eloss.root") 
+    : fFitter(0), 
+      fAxis(0),
+      fFits("AliFMDCorrELossFit::ELossFit"),
+      fSys(sys), 
+      fCMS(cms), 
+      fField(field), 
+      fMC(mc)
+  {
+    TFile* file = TFile::Open(filename, "READ");
+    if (!file) { 
+      Error("MakeELossFit", "Failed to open file %s", filename);
+      return;
+    }
+    TList* forward = static_cast<TList*>(file->Get("Forward"));
+    // static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
+    if (!forward) { 
+      Error("MakeELossFit", "Couldn't get forward list from %s", filename);
+      return;
+    }
+
+    fFitter = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
+    if (!fFitter) { 
+      Error("MakeELossFit", "Couldn't get fitter folder");
+      return;
+    }
+
+    fAxis = static_cast<TAxis*>(fFitter->FindObject("etaAxis"));
+    if (!fAxis) { 
+      Error("MakeELossFit", "Couldn't get eta axis");
+      return;
+    }
+    file->Close();
+
+#if 0
+    AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+    mgr.Init(sys, cms, field, mc, 0, true);
+#endif
+  }
+  //__________________________________________________________________
+  TList* GetRing(UShort_t d, Char_t r) const
+  {
+    TList* rL = static_cast<TList*>(fFitter->FindObject(Form("FMD%d%c",d,r)));
+    if (!rL) { 
+      Warning("DrawFits", "List FMD%d%c not found", d, r);
+      return 0;
+    }
+    return rL;
+  }
+  //__________________________________________________________________
+  TList* GetEDists(UShort_t d, Char_t r) const
+  {
+    TList* rList = GetRing(d, r);
+    if (!rList) 
+      return 0;
+    // rList->ls();
+
+    TList* edists = static_cast<TList*>(rList->FindObject("EDists"));
+    if (!edists) { 
+      Error("PrintFits", "Couldn't get EDists list for FMD%d%c", d,r);
+      return 0; 
+    }
+    return edists;
+  }
+  //__________________________________________________________________
+  TH1* GetEDist(UShort_t d, Char_t r, UShort_t etabin)
+  {
+    TList* eList = GetEDists(d, r);
+    if (!eList) { 
+      Warning("GetEDist", "No list for FMD%d%c", d, r);
+      return 0;
+    }
+    // eList->ls();
+
+    TString cmp(Form("FMD%d%c_etabin%03d", d, r, etabin));
+
+    TIter next(eList);
+    TObject* o = 0;
+    while ((o = next())) { 
+      if (!cmp.CompareTo(o->GetName())) {
+       return static_cast<TH1*>(o);
+      }
+    }
+    return 0;
+  }
+#ifndef __CINT__
+  //__________________________________________________________________
+  AliFMDCorrELossFit::ELossFit* FindBestFit(TH1* dist) 
+  {
+    TList* funcs = dist->GetListOfFunctions();
+    TIter  next(funcs);
+    TF1*   func = 0;
+    fFits.Clear();
+    Int_t  i = 0;
+    Info("FindBestFit", "%s", dist->GetName());
+    while ((func = static_cast<TF1*>(next()))) { 
+      AliFMDCorrELossFit::ELossFit* fit = 
+       new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
+      fit->CalculateQuality(10, .20, 1e-7);
+    }
+    fFits.Sort(false);
+    // fFits.Print();
+    return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
+  }
+#endif
+  //__________________________________________________________________
+  Bool_t Run()
+  {
+    if (!fFitter || !fAxis) { 
+      Error("Run", "Missing objects");
+      return kFALSE;
+    }
+    AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
+    obj->SetEtaAxis(*fAxis);
+
+    for (UShort_t d = 1; d <= 3; d++) { 
+      Info("Run", "detector is FMD%d", d);
+      UShort_t nQ = (d == 1 ? 1 : 2);
+      for (UShort_t q = 0; q < nQ; q++) { 
+       Char_t r = (q == 0 ? 'I' : 'O');
+       Int_t nBin = fAxis->GetNbins();
+       Info("Run", " ring is FMD%d%c - %d bins", d, r, nBin);
+       
+       Bool_t oneSeen = kFALSE;
+       for (UShort_t b = 1; b <= nBin; b++) { 
+         TH1* h = GetEDist(d, r, b);
+         if (oneSeen && !h) break;
+         if (!h)            continue;
+         if (!oneSeen)      oneSeen = true;
+
+         AliFMDCorrELossFit::ELossFit* best = FindBestFit(h);
+         best->Print("");
+         best->fDet  = d; 
+         best->fRing = r;
+         best->fBin  = b;
+         // Double_t eta = fAxis->GetBinCenter(b);
+         Info("Run", "    bin=%3d->%6.4f", b, eta);
+         obj->SetFit(d, r, b, new AliFMDCorrELossFit::ELossFit(*best));
+       }
+      }
+    }
+    AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+    TString fname(mgr.GetFileName(AliForwardCorrectionManager::kELossFits,
+                                 fSys, fCMS, fField, fMC));
+    TFile* output = TFile::Open(fname.Data(), "RECREATE");
+    if (!output) { 
+      Warning("Run", "Failed to open output file %s", fname.Data());
+      return kFALSE;
+    }
+    obj->Write(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
+    output->Write();
+    output->Close();
+    Info("Run", "File %s created.  It should be copied to %s and stored in SVN",
+        fname.Data(),mgr.GetFileDir(AliForwardCorrectionManager::kELossFits));
+    
+    return kTRUE;
+  }
+};
+//
+// EOF
+//
diff --git a/PWG2/FORWARD/analysis2/tests/TestMarkers.C b/PWG2/FORWARD/analysis2/tests/TestMarkers.C
new file mode 100644 (file)
index 0000000..69158fb
--- /dev/null
@@ -0,0 +1,82 @@
+namespace {
+  enum { 
+    kSolid        = 0x000, 
+    kHollow       = 0x001, 
+    kCircle       = 0x002,
+    kSquare       = 0x004, 
+    kUpTriangle   = 0x006, 
+    kDownTriangle = 0x008, 
+    kDiamond      = 0x00a,
+    kCross        = 0x00c,
+    kStar         = 0x00e
+  };
+  Int_t MarkerStyle(UInt_t bits)
+  {
+    Int_t  base   = bits & (0xFE);
+    Bool_t hollow = bits & kHollow;
+    switch (base) { 
+    case kCircle:       return (hollow ? 24 : 20);
+    case kSquare:       return (hollow ? 25 : 21);
+    case kUpTriangle:   return (hollow ? 26 : 22);
+    case kDownTriangle: return (hollow ? 32 : 23);
+    case kDiamond:      return (hollow ? 27 : 33); 
+    case kCross:        return (hollow ? 28 : 34); 
+    case kStar:         return (hollow ? 30 : 29); 
+    }
+    return 1;
+  }
+  UShort_t MarkerBits(Int_t style) 
+  { 
+    UShort_t bits = 0;
+    switch (style) { 
+    case 24: case 25: case 26: case 27: case 28: case 30: case 32: 
+      bits |= kHollow; break;
+    }
+    switch (style) { 
+    case 20: case 24: bits |= kCircle;       break;
+    case 21: case 25: bits |= kSquare;       break;
+    case 22: case 26: bits |= kUpTriangle;   break;
+    case 23: case 32: bits |= kDownTriangle; break;
+    case 27: case 33: bits |= kDiamond;      break;
+    case 28: case 34: bits |= kCross;        break;
+    case 29: case 30: bits |= kStar;         break;
+    }
+    return bits;
+  }
+  Int_t FlipHollow(Int_t style) 
+  {
+    UShort_t bits = MarkerBits(style);
+    Int_t ret = MarkerStyle(bits ^ kHollow);
+    Info("FlipHollow", "style=%2d -> bits=0x%02x -> mask=0x%02x -> ret=%02d", 
+        style, bits, (bits ^ kHollow), ret);
+    return ret;
+  }
+}
+
+void DrawOne(const char* what, UShort_t base, Double_t y)
+{
+  TLatex* l = new TLatex(.07, y, what);
+  l->SetTextSize(0.03);
+  l->Draw();
+
+  Int_t filled = MarkerStyle(base);
+  // Info("DrawOne", "%2d (%16s) -> %d", base, what, style);
+  TMarker* p = new TMarker(.35, y, filled);
+  p->SetMarkerSize(1.5);
+  p->Draw();
+
+  Int_t hollow = MarkerStyle(base|kHollow);
+  p = new TMarker(.60, y, hollow);
+  p->SetMarkerSize(1.5);
+  p->Draw();
+
+  p = new TMarker(.75, y, FlipHollow(filled));
+  p->SetMarkerSize(1.5);
+  p->Draw();
+
+  p = new TMarker(.85, y, FlipHollow(hollow));
+  p->SetMarkerSize(1.5);
+  p->Draw();
+    
+}
+
diff --git a/PWG2/FORWARD/analysis2/tests/TestPoisson.C b/PWG2/FORWARD/analysis2/tests/TestPoisson.C
new file mode 100644 (file)
index 0000000..1ff19e4
--- /dev/null
@@ -0,0 +1,188 @@
+void
+MakeIntegerAxis(Int_t& nBins, Double_t& min, Double_t& max)
+{
+  if (nBins <= 0) return;
+  if (max <= 0)   max = nBins;
+  
+  Double_t dBin = max / nBins;
+  min   = - dBin / 2;
+  max   = max + dBin / 2;
+  nBins = nBins + 1;
+}
+
+void
+TestPoisson(Double_t o=.3, bool useWeights=false)
+{
+  const char* load = "$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C";
+  gROOT->Macro(load);
+  gROOT->GetInterpreter()->UnloadFile(gSystem->ExpandPathName(load));
+  
+  // --- Parameters of this script -----------------------------------
+  Int_t      nBin =  5;  // Our detector matrix size 
+  Int_t      nMax = TMath::Max(Int_t(nBin * nBin * o + .5)+nBin/2,nBin);  
+  Int_t      nEv  = 10000; // Number of events
+  Double_t   mp   = o;   // The 'hit' probability 
+
+
+  TH2D* base = new TH2D("base", "Basic histogram", 
+                       nBin,-.5, nBin-.5, nBin, -.5, nBin-.5);
+  base->SetXTitle("#eta");
+  base->SetYTitle("#varphi");
+  base->SetDirectory(0);
+  base->SetOption("colz");
+
+  Int_t tN1=nMax;    Double_t tMin1; Double_t tMax1;
+  Int_t tN2=nMax*10; Double_t tMin2; Double_t tMax2=nMax;
+  MakeIntegerAxis(tN1, tMin1, tMax1);
+  MakeIntegerAxis(tN2, tMin2, tMax2);
+  TH2D* corr = new TH2D("comp", "Comparison", 
+                       tN1, tMin1, tMax1, tN2, tMin2, tMax2);
+  corr->SetXTitle("Input");
+  corr->SetYTitle("Poisson");
+  corr->SetDirectory(0);
+  corr->SetOption("colz");
+  corr->SetStats(0);
+  TLine* lcorr = new TLine(0, 0, tMax2, tMax2);
+
+  Int_t mm = TMath::Max(Int_t(nBin * o + .5),nBin/2);
+  tN2=mm*10; tMax2 = mm;
+  MakeIntegerAxis(tN2, tMin2, tMax2);
+  Info("", "Making mean w/nbins=%d,range=[%f,%f]", tN2, tMin2, tMax2);
+  TH2D* mean = new TH2D("mean", "Mean comparison", 
+                       tN2, tMin2, tMax2, tN2, tMin2, tMax2);
+  mean->SetXTitle("Input");
+  mean->SetYTitle("Poisson");
+  mean->SetDirectory(0);
+  mean->SetOption("colz");
+  mean->SetStats(0);
+  TLine* lmean = new TLine(tMin2, tMin2, tMax2, tMax2);
+
+  TH1D* dist = new TH1D("dist", "Distribution of hits", tN1, tMin1, tMax1);
+  dist->SetXTitle("s");
+  dist->SetYTitle("P(s)");
+  dist->SetFillColor(kRed+1);
+  dist->SetFillStyle(3001);
+  dist->SetDirectory(0);
+                       
+  AliPoissonCalculator* c = new AliPoissonCalculator("ignored");
+  c->Init(nBin ,nBin);
+
+  for (Int_t i = 0; i < nEv; i++) { 
+    c->Reset(base);
+    base->Reset();
+
+    for (Int_t iEta = 0; iEta < nBin; iEta++) { 
+      for (Int_t iPhi = 0; iPhi < nBin; iPhi++) { 
+       // Throw a die 
+       Int_t m = gRandom->Poisson(mp);
+       dist->Fill(m);
+
+       // Fill into our base histogram 
+       base->Fill(iEta, iPhi, m);
+
+       // Fill into poisson calculator 
+       c->Fill(iEta, iPhi, m > 0, (useWeights ? m : 1));
+      }
+    }
+    // Calculate the result 
+    TH2D* res = c->Result();
+    
+    // Now loop and compare 
+    Double_t mBase = 0;
+    Double_t mPois = 0;
+    for (Int_t iEta = 0; iEta < nBin; iEta++) { 
+      for (Int_t iPhi = 0; iPhi < nBin; iPhi++) { 
+       Double_t p = res->GetBinContent(iEta, iPhi);
+       Double_t t = base->GetBinContent(iEta, iPhi);
+
+       mBase += t;
+       mPois += p;
+       corr->Fill(t, p);
+      }
+    }
+    Int_t nn = nBin * nBin;
+    mean->Fill(mBase / nn, mPois / nn);
+  }
+
+  TCanvas* cc = new TCanvas("c", "c", 900, 900);
+  cc->SetFillColor(0);
+  cc->SetFillStyle(0);
+  cc->SetBorderMode(0);
+  cc->SetRightMargin(0.02);
+  cc->SetTopMargin(0.02);
+  cc->Divide(2,2);
+  
+  TVirtualPad* pp = cc->cd(1);
+  pp->SetFillColor(0);
+  pp->SetFillStyle(0);
+  pp->SetBorderMode(0);
+  pp->SetRightMargin(0.15);
+  pp->SetTopMargin(0.02);
+  pp->SetLogz();
+  pp->SetGridx();
+  pp->SetGridy();
+  corr->Draw();
+  lcorr->Draw();
+
+  pp = cc->cd(2);
+  pp->SetFillColor(0);
+  pp->SetFillStyle(0);
+  pp->SetBorderMode(0);
+  pp->SetRightMargin(0.02);
+  pp->SetTopMargin(0.02);
+#if 1
+  c->GetMean()->Draw();
+#elif 1
+  c->GetOccupancy()->Draw();
+#else
+  pp->SetLogy();
+  dist->SetStats(0);
+  dist->Scale(1. / dist->Integral());
+  dist->Draw();
+  TH1D* m1 = c->GetMean();
+  m1->Scale(1. / m1->Integral());
+  m1->Draw("same");
+  Double_t eI;
+  Double_t ii = 100 * dist->Integral(2, 0);
+  TLatex* ll = new TLatex(.97, .85, 
+                         Form("Input #bar{m}: %5.3f", mp));
+  ll->SetNDC();
+  ll->SetTextFont(132);
+  ll->SetTextAlign(31);
+  ll->Draw();
+  ll->DrawLatex(.97, .75, Form("Result #bar{m}: %5.3f", dist->GetMean()));
+  ll->DrawLatex(.97, .65, Form("Occupancy: #int_{1}^{#infty}P(s)ds = %6.2f%%",
+                              ii));
+                        
+#endif
+
+  pp = cc->cd(3);
+  pp->SetFillColor(0);
+  pp->SetFillStyle(0);
+  pp->SetBorderMode(0);
+  pp->SetRightMargin(0.15);
+  pp->SetTopMargin(0.02);
+  pp->SetGridx();
+  pp->SetGridy();
+  c->GetCorrection()->Draw();
+
+  pp = cc->cd(4);
+  pp->SetFillColor(0);
+  pp->SetFillStyle(0);
+  pp->SetBorderMode(0);
+  pp->SetRightMargin(0.15);
+  pp->SetTopMargin(0.02);
+  pp->SetLogz();
+  pp->SetGridx();
+  pp->SetGridy();
+  mean->Draw();
+  lmean->Draw();
+
+  cc->cd();
+}
+
+//
+// EOF
+//
+
+
diff --git a/PWG2/FORWARD/analysis2/tests/TestRunMakeELossFit.C b/PWG2/FORWARD/analysis2/tests/TestRunMakeELossFit.C
new file mode 100644 (file)
index 0000000..d32dce3
--- /dev/null
@@ -0,0 +1,43 @@
+/** 
+ * Run the energy loss fit finder and generate corrections output file 
+ * 
+ * @param sys       Collision system 
+ * @param cms       Center of mass energy per nucleon in GeV
+ * @param field     Magnetic field 
+ * @param mc        Whether this is for Monte-Carlo data
+ * @param filename  Input file name 
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ *
+ * @depcrecated 
+ * The class AliFMDELossFitter automatically generates the
+ * AliFMDCorrELossFit object.
+ *
+ */
+void
+TestRunMakeELossFit(UShort_t    sys, 
+                   UShort_t    cms, 
+                   Short_t     field, 
+                   Bool_t      mc=false,
+                   const char* filename="forward_eloss.root")
+{
+  std::cout << "Loading libraries ..." << std::endl;
+  gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+
+  std::cout << "Loading compile script ..." << std::endl;
+  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/Compile.C");
+  std::cout << "Compiling MakeELossFit.C script ..." << std::endl;
+  Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/TestMakeELossFit.C"); 
+
+  std::cout << "Making MakeELossFit object (sys=" << sys 
+           << ", cms=" << cms << ", field=" << field << ", mc=" << mc 
+           << ")" << std::endl;
+  MakeELossFit mef(sys, cms, field, mc, filename); 
+
+  std::cout << "Running maker ..." << std::endl;
+  mef.Run();
+}
+//
+// EOF
+//