Re-organised scripts
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Jul 2011 08:31:31 +0000 (08:31 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Jul 2011 08:31:31 +0000 (08:31 +0000)
45 files changed:
PWG2/FORWARD/analysis2/corrs/CompELossFits.C [moved from PWG2/FORWARD/analysis2/scripts/CompELossFits.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/CompareCentralSecMaps.C [moved from PWG2/FORWARD/analysis2/scripts/CompareCentralSecMaps.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/CompareCorrs.C [moved from PWG2/FORWARD/analysis2/scripts/CompareCorrs.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/CompareSecMaps.C [moved from PWG2/FORWARD/analysis2/scripts/CompareSecMaps.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/CompareVtxBias.C [moved from PWG2/FORWARD/analysis2/scripts/CompareVtxBias.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/DrawAnaELoss.C [moved from PWG2/FORWARD/analysis2/scripts/DrawAnaELoss.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/DrawCorrAcc.C [moved from PWG2/FORWARD/analysis2/scripts/DrawCorrAcc.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/DrawCorrAcc2.C [moved from PWG2/FORWARD/analysis2/scripts/DrawCorrAcc2.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/DrawCorrCentralSecMap2.C [moved from PWG2/FORWARD/analysis2/scripts/DrawCorrCentralSecMap2.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/DrawCorrELoss.C [moved from PWG2/FORWARD/analysis2/scripts/DrawCorrELoss.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/DrawCorrSecMap.C [moved from PWG2/FORWARD/analysis2/scripts/DrawCorrSecMap.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/DrawCorrSecMap2.C [moved from PWG2/FORWARD/analysis2/scripts/DrawCorrSecMap2.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/DrawCorrVtxBias.C [moved from PWG2/FORWARD/analysis2/scripts/DrawCorrVtxBias.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/ExtractELoss.C [moved from PWG2/FORWARD/analysis2/scripts/ExtractELoss.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/ExtractSecMap.C [moved from PWG2/FORWARD/analysis2/scripts/ExtractSecMap.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/MakeAcceptanceCorrection.C [moved from PWG2/FORWARD/analysis2/scripts/MakeAcceptanceCorrection.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/MakeCorrRepository.C [moved from PWG2/FORWARD/analysis2/scripts/MakeCorrRepository.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/MakeCorrSecMap.C [moved from PWG2/FORWARD/analysis2/scripts/MakeCorrSecMap.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/MakeELossFitsFromFile.C [moved from PWG2/FORWARD/analysis2/scripts/MakeELossFitsFromFile.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/MoveCorrections.C [moved from PWG2/FORWARD/analysis2/scripts/MoveCorrections.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/RunCopyCentralSecMap.C [moved from PWG2/FORWARD/analysis2/scripts/RunCopyCentralSecMap.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/RunCopyELossFit.C [moved from PWG2/FORWARD/analysis2/scripts/RunCopyELossFit.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/RunCopyMergeEff.C [moved from PWG2/FORWARD/analysis2/scripts/RunCopyMergeEff.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/RunCopySecMap.C [moved from PWG2/FORWARD/analysis2/scripts/RunCopySecMap.C with 100% similarity]
PWG2/FORWARD/analysis2/corrs/RunCopyVtxBias.C [moved from PWG2/FORWARD/analysis2/scripts/RunCopyVtxBias.C with 100% similarity]
PWG2/FORWARD/analysis2/qa/DrawBeforeAfter.C [moved from PWG2/FORWARD/analysis2/scripts/DrawBeforeAfter.C with 100% similarity]
PWG2/FORWARD/analysis2/qa/DrawCuts.C [moved from PWG2/FORWARD/analysis2/scripts/DrawCuts.C with 100% similarity]
PWG2/FORWARD/analysis2/qa/DrawELossPoisson.C [moved from PWG2/FORWARD/analysis2/scripts/DrawELossPoisson.C with 87% similarity]
PWG2/FORWARD/analysis2/qa/DrawMCResult.C [moved from PWG2/FORWARD/analysis2/scripts/DrawMCResult.C with 100% similarity]
PWG2/FORWARD/analysis2/qa/DrawNeighbors.C [moved from PWG2/FORWARD/analysis2/scripts/DrawNeighbors.C with 100% similarity]
PWG2/FORWARD/analysis2/qa/DrawRecAnaEloss.C [moved from PWG2/FORWARD/analysis2/scripts/DrawRecAnaEloss.C with 100% similarity]
PWG2/FORWARD/analysis2/qa/DrawSteps.C [moved from PWG2/FORWARD/analysis2/scripts/DrawSteps.C with 100% similarity]
PWG2/FORWARD/analysis2/scripts/LoadLibs.C
PWG2/FORWARD/analysis2/scripts/TestAcc.C [deleted file]
PWG2/FORWARD/analysis2/scripts/TestELossDist.C [deleted file]
PWG2/FORWARD/analysis2/scripts/TestFitELoss.C [deleted file]
PWG2/FORWARD/analysis2/scripts/TestMakeELossFits.C [deleted file]
PWG2/FORWARD/analysis2/scripts/TestMarkers.C [deleted file]
PWG2/FORWARD/analysis2/scripts/TestRunMakeELossFit.C [deleted file]
PWG2/FORWARD/analysis2/trains/BuildTrain.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/trains/MakeAODTrain.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/trains/MakeFMDELossTrain.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/trains/MakeMCCorrTrain.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/trains/MakedNdetaTrain.C [new file with mode: 0644]
PWG2/FORWARD/analysis2/trains/TrainSetup.C [new file with mode: 0644]

@@ -1,4 +1,25 @@
-
+/**
+ * @file   DrawELossPoisson.C
+ * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
+ * @date   Thu Jul  7 10:24:58 2011
+ * 
+ * @brief  A script to draw the Poisson vs Energy Loss correlation 
+ * 
+ * 
+ */
+
+/** 
+ * Draw the poisson @f$N_{ch}@f$ estimate against the @f$\Delta@f$
+ * @f$N_{ch}@f$ estimate and do a regression line between the two 
+ * 
+ * @param p            List 
+ * @param d            Detector
+ * @param r            Ring
+ * @param xmin         Minimum
+ * @param xmax         Maximum
+ * 
+ * @return The regression coefficient 
+ */
 Double_t
 DrawRingELossPoisson(TList* p, UShort_t d, Char_t r, 
                     Double_t xmin, Double_t xmax)
@@ -124,7 +145,15 @@ DrawRingELossPoisson(TList* p, UShort_t d, Char_t r,
   return corr->GetCorrelationFactor();
 }
 
-
+/** 
+ * Draw the correlation between the Poisson @f$N_{ch}@f$ estimate
+ * and the @f$\Delta@f$ @f$N_{ch}@f$ estimate and do a regression
+ * line between the two for each ring
+ * 
+ * @param filename File to read
+ * @param xmax     Minimum X
+ * @param xmin     Maximum X 
+ */
 void
 DrawELossPoisson(const char* filename="forward.root", 
                 Double_t xmax=-1,
@@ -202,7 +231,6 @@ DrawELossPoisson(const char* filename="forward.root",
   c->cd();
   c->SaveAs("elossVsPoisson.png");
 }
-
-  
-  
+//
+// EOF
+//
index 2e10a82..6353954 100644 (file)
@@ -6,6 +6,14 @@
 void
 LoadLibs()
 {
+  gROOT->LoadClass("TVirtualMC",           "libVMC");
+  gROOT->LoadClass("AliVEvent",            "libSTEERBase");
+  gROOT->LoadClass("AliESDEvent",          "libESD");
+  gROOT->LoadClass("AliAnalysisManager",   "libANALYSIS");
+  gROOT->LoadClass("AliAnalysisTaskSE",    "libANALYSISalice");
+  gROOT->LoadClass("AliAODForwardMult",    "libPWG2forward2");
+
+#if 0
   const char* test = gSystem->GetLibraries("PWG2forward2","D",false);
   if (test && test[0] != '\0') { 
     // TInterpreter* inter = gROOT->GetInterpreter();
@@ -21,6 +29,7 @@ LoadLibs()
   gSystem->Load("libANALYSISalice");
   gSystem->Load("libPWG0base");
   gSystem->Load("libPWG2forward2");
+#endif
 }
 //
 // EOF
diff --git a/PWG2/FORWARD/analysis2/scripts/TestAcc.C b/PWG2/FORWARD/analysis2/scripts/TestAcc.C
deleted file mode 100644 (file)
index b012f93..0000000
+++ /dev/null
@@ -1,244 +0,0 @@
-#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/scripts/TestELossDist.C b/PWG2/FORWARD/analysis2/scripts/TestELossDist.C
deleted file mode 100644 (file)
index d48c538..0000000
+++ /dev/null
@@ -1,1426 +0,0 @@
-#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/scripts/TestFitELoss.C b/PWG2/FORWARD/analysis2/scripts/TestFitELoss.C
deleted file mode 100644 (file)
index ae5a975..0000000
+++ /dev/null
@@ -1,305 +0,0 @@
-/**
- * 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/scripts/TestMakeELossFits.C b/PWG2/FORWARD/analysis2/scripts/TestMakeELossFits.C
deleted file mode 100644 (file)
index 16f47f5..0000000
+++ /dev/null
@@ -1,232 +0,0 @@
-#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/scripts/TestMarkers.C b/PWG2/FORWARD/analysis2/scripts/TestMarkers.C
deleted file mode 100644 (file)
index 69158fb..0000000
+++ /dev/null
@@ -1,82 +0,0 @@
-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/scripts/TestRunMakeELossFit.C b/PWG2/FORWARD/analysis2/scripts/TestRunMakeELossFit.C
deleted file mode 100644 (file)
index d32dce3..0000000
+++ /dev/null
@@ -1,43 +0,0 @@
-/** 
- * 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
-//
diff --git a/PWG2/FORWARD/analysis2/trains/BuildTrain.C b/PWG2/FORWARD/analysis2/trains/BuildTrain.C
new file mode 100644 (file)
index 0000000..6e79ff3
--- /dev/null
@@ -0,0 +1,66 @@
+Bool_t
+BuildTrain(const char* script, const char* extra="", Bool_t useTmp=false)
+{
+  // --- Load needed libraries ---------------------------------------
+  gROOT->LoadClass("TVirtualMC",           "libVMC");
+  gROOT->LoadClass("AliVEvent",            "libSTEERBase");
+  gROOT->LoadClass("AliESDEvent",          "libESD");
+  gROOT->LoadClass("AliAnalysisManager",   "libANALYSIS");
+  gROOT->LoadClass("AliAnalysisTaskSE",    "libANALYSISalice");
+  gROOT->LoadClass("AliAODForwardMult",    "libPWG2forward2");
+
+  // --- Setup script path -------------------------------------------
+  const char* aliPath   = gSystem->ExpandPathName("$ALICE_ROOT");
+  const char* fwd2Path  = 
+    gSystem->ExpandPathName("$ALICE_ROOT/PWG2/FORWARD/analysis2");
+  const char* macroPath = gROOT->GetMacroPath();
+  gROOT->SetMacroPath(".:%s:%s/trains:%s/scripts",macroPath,fwd2Path,fwd2Path);
+  
+  // --- Setup include path ------------------------------------------
+  gSystem->AddIncludePath(Form("-I%s", fwd2Path));
+  gSystem->AddIncludePath(Form("-I%s", aliPath));
+  gSystem->AddIncludePath(Form("-I%s/include", aliPath));
+
+  // --- Check that we can find the script ---------------------------
+  TString path = gSystem->Which(gROOT->GetMacroPath(), script);
+  if (path.IsNull()) { 
+    path = gSystem->Which(gROOT->GetMacroPath(), Form("%s.C", script));
+    if (path.IsNull()) { 
+      Error("BuildTrain", "Cannot find script %s or %s.C in %s", 
+           script, script, gROOT->GetMacroPath());
+      return false;
+    }
+  }
+  
+  // --- Possibly make a temporary copy ------------------------------
+  if (useTmp) { 
+    TString tmp("Train");
+    FILE*   fp = gSystem->TempFileName(tmp, ".");
+    fclose(fp);
+    gSystem->Unlink(tmp);
+    tmp.Append(".C");
+    
+    gSystem->CopyFile(path, tmp);
+    path = tmp;
+  }
+
+  // --- Now compile the thing ---------------------------------------
+  Int_t ret = gROOT->LoadMacro(Form("%s+%s", path.Data(), extra));
+
+  // --- If we did a temporary file, remove stuff --------------------
+  if (useTmp) {
+    gSystem->Unlink(path);
+    path.ReplaceAll(".C",   "_C.d");  gSystem->Unlink(path);
+    path.ReplaceAll("_C.d", "_C.so"); gSystem->Unlink(path);
+  }
+
+  // --- Warning in case of problems ---------------------------------
+  if (ret != 0) 
+    Warning("BuildSetup", "Failed to build %s (%s)", path.Data(), script);
+
+  return ret == 0;
+}
+//
+// EOF
+//
+
diff --git a/PWG2/FORWARD/analysis2/trains/MakeAODTrain.C b/PWG2/FORWARD/analysis2/trains/MakeAODTrain.C
new file mode 100644 (file)
index 0000000..87f20df
--- /dev/null
@@ -0,0 +1,170 @@
+//====================================================================
+/**
+ * Analysis train to make Forward and Central multiplicity
+ * 
+ * To run, do 
+ * @code 
+ * gROOT->LoadMacro("TrainSetup.C");
+ * // Make train 
+ * MakeAODTrain t("My Analysis");
+ * // Set variaous parameters on the train 
+ * t.SetDataDir("/home/of/data");
+ * t.AddRun(118506)
+ * // Run it 
+ * t.Run("LOCAL", "FULL", -1, false, false);
+ * @endcode 
+ *
+ * @ingroup pwg2_forward_scripts_makers
+ * @ingroup pwg2_forward_aod
+ */
+class MakeAODTrain : public TrainSetup
+{
+public:
+  /** 
+   * Constructor.  Date and time must be specified when running this
+   * in Termiante mode on Grid
+   * 
+   * @param name     Name of train (free form)
+   * @param sys      Collision system (1: pp, 2: PbPb)
+   * @param sNN      Center of mass energy [GeV]
+   * @param field    L3 magnetic field - one of {-5,0,+5} kG
+   * @param useCent  Whether to use centrality 
+   * @param dateTime Append date and time to name 
+   * @param year     Year     - if not specified, current year
+   * @param month    Month    - if not specified, current month
+   * @param day      Day      - if not specified, current day
+   * @param hour     Hour     - if not specified, current hour
+   * @param min      Minutes  - if not specified, current minutes
+   */
+  MakeAODTrain(const  char* name, 
+              UShort_t     sys      = 0, 
+              UShort_t     sNN      = 0, 
+              Short_t      field    = 0, 
+              Bool_t       useCent  = false, 
+              Bool_t       dateTime = false, 
+              UShort_t     year     = 0, 
+              UShort_t     month    = 0, 
+              UShort_t     day      = 0, 
+              UShort_t     hour     = 0, 
+              UShort_t     min      = 0) 
+    : TrainSetup(name, dateTime, 
+                year, month, day, hour, min),
+      fSys(sys), 
+      fSNN(sNN), 
+      fField(field),
+      fUseCent(useCent)
+  {}
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode - see TrainSetup::EMode
+   * @param oper     Operation - see TrainSetup::EOperation
+   * @param nEvents  Number of events (negative means all)
+   * @param mc       If true, assume simulated events 
+   * @param usePar   If true, use PARs 
+   */
+  void Run(const char* mode, const char* oper, 
+          Int_t nEvents=-1, Bool_t mc=false,
+          Bool_t usePar=false)
+  {
+    Info("Run", "Running in mode=%s, oper=%s, events=%d, MC=%d, Par=%d", 
+        mode, oper, nEvents, mc, usePar);
+    Exec("ESD", mode, oper, nEvents, mc, usePar);
+  }
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode - see TrainSetup::EMode
+   * @param oper     Operation - see TrainSetup::EOperation
+   * @param nEvents  Number of events (negative means all)
+   * @param mc       If true, assume simulated events 
+   * @param usePar   If true, use PARs 
+   */
+  void Run(EMode mode, EOper oper, Int_t nEvents=-1, Bool_t mc=false, 
+          Bool_t usePar = false)
+  {
+    Info("Run", "Running in mode=%d, oper=%d, events=%d, MC=%d, Par=%d", 
+        mode, oper, nEvents, mc, usePar);
+    Exec(kESD, mode, oper, nEvents, mc, usePar);
+  }
+protected:
+  /** 
+   * Create the tasks 
+   * 
+   * @param mode Processing mode
+   * @param par  Whether to use par files 
+   * @param mgr  Analysis manager 
+   */
+  void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)
+  {
+    // --- Output file name ------------------------------------------
+    AliAnalysisManager::SetCommonFileName("forward.root");
+
+    // --- Load libraries/pars ---------------------------------------
+    LoadLibrary("PWG2forward2", mode, par, true);
+    
+    // --- Set load path ---------------------------------------------
+    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+                            gROOT->GetMacroPath()));
+
+    // --- Check if this is MC ---------------------------------------
+    Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
+    
+    // --- Add the task ----------------------------------------------
+    gROOT->Macro(Form("AddTaskForwardMult.C(%d,%d,%d,%d)", 
+                     mc, fSys, fSNN, fField));
+    AddExtraFile(gSystem->Which(gROOT->GetMacroPath(), "ForwardAODConfig.C"));
+
+    // --- Add the task ----------------------------------------------
+    gROOT->Macro(Form("AddTaskCentralMult.C(%d,%d,%d,%d)", 
+                     mc, fSys, fSNN, fField));
+    AddExtraFile(gSystem->Which(gROOT->GetMacroPath(), "CentralAODConfig.C"));
+  }
+  //__________________________________________________________________
+  /** 
+   * Create physics selection , and add to manager
+   * 
+   * @param mc Whether this is for MC 
+   * @param mgr Manager 
+   */
+  void CreatePhysicsSelection(Bool_t mc,
+                             AliAnalysisManager* mgr)
+  {
+    TrainSetup::CreatePhysicsSelection(mc, mgr);
+
+    // --- Get input event handler -----------------------------------
+    AliInputEventHandler* ih =
+      dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
+    if (!ih) 
+      Fatal("CreatePhysicsSelection", "Couldn't get input handler (%p)", ih);
+    
+    // --- Get Physics selection -------------------------------------
+    AliPhysicsSelection* ps = 
+      dynamic_cast<AliPhysicsSelection*>(ih->GetEventSelection());
+    if (!ps) 
+      Fatal("CreatePhysicsSelection", "Couldn't get PhysicsSelection (%p)",ps);
+
+    // --- Ignore trigger class when selecting events.  This means ---
+    // --- that we get offline+(A,C,E) events too --------------------
+    // ps->SetSkipTriggerClassSelection(true);
+  }
+  //__________________________________________________________________
+  /** 
+   * Create the centrality selection only if requested
+   * 
+   * @param mc  Monte-Carlo truth flag 
+   * @param mgr Manager
+   */
+  void CreateCentralitySelection(Bool_t mc, AliAnalysisManager* mgr)
+  {
+    if (!fUseCent) return;
+    TrainSetup::CreateCentralitySelection(mc, mgr);
+  }
+  UShort_t fSys;
+  UShort_t fSNN;
+  Short_t  fField;
+  Bool_t   fUseCent;
+};
+//
+// EOF
+//
diff --git a/PWG2/FORWARD/analysis2/trains/MakeFMDELossTrain.C b/PWG2/FORWARD/analysis2/trains/MakeFMDELossTrain.C
new file mode 100644 (file)
index 0000000..c191439
--- /dev/null
@@ -0,0 +1,119 @@
+//====================================================================
+/**
+ * Analysis train to do energy loss fits
+ * 
+ * @ingroup pwg2_forward_scripts_makers
+ */
+class MakeFMDELossTrain : public TrainSetup
+{
+public:
+  /** 
+   * Constructor.  Date and time must be specified when running this
+   * in Termiante mode on Grid
+   * 
+   * @param name     Name of train 
+   * @param useCent  Whether to use centrality or not 
+   * @param dateTime Append date and time to name
+   * @param year     Year
+   * @param month    Month 
+   * @param day      Day
+   * @param hour     Hour 
+   * @param min      Minutes
+   */
+  MakeFMDELossTrain(const char* name  = "FMD Energy Loss",
+                   Bool_t   useCent  = false,
+                   Bool_t   dateTime = false, 
+                   UShort_t year     = 0, 
+                   UShort_t month    = 0, 
+                   UShort_t day      = 0, 
+                   UShort_t hour     = 0, 
+                   UShort_t min      = 0) 
+    : TrainSetup(name, dateTime, year, month, day, hour, min), 
+      fUseCent(useCent)
+  {}
+  //__________________________________________________________________
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode
+   * @param oper     Operation
+   * @param nEvents  Number of events (negative means all)
+   * @param mc       If true, assume simulated events 
+   * @param par      IF true, use par files 
+   */
+  void Run(const char* mode, const char* oper, 
+          Int_t nEvents=-1, Bool_t mc=false, Bool_t par=false)
+  {
+    Exec("ESD", mode, oper, nEvents, mc, par);
+  }
+  //__________________________________________________________________
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode
+   * @param oper     Operation
+   * @param nEvents  Number of events (negative means all)
+   * @param mc       If true, assume simulated events 
+   * @param par      IF true, use par files 
+   */
+  void Run(EMode mode, EOper oper, Int_t nEvents=-1, Bool_t mc=false,
+          Bool_t par=false)
+  {
+    Exec(kESD, mode, oper, nEvents, mc, par);
+  }
+protected:
+  //__________________________________________________________________
+  /** 
+   * Create the tasks 
+   * 
+   * @param mode Processing mode
+   * @param par  Whether to use par files 
+   * @param mgr  Analysis manager 
+   */
+  void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)
+  {
+    // --- Output file name ------------------------------------------
+    AliAnalysisManager::SetCommonFileName("forward_eloss.root");
+
+    // --- Load libraries/pars ---------------------------------------
+    LoadLibrary("PWG2forward2", mode, par, true);
+    
+    // --- Set load path ---------------------------------------------
+    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+                            gROOT->GetMacroPath()));
+
+    // --- Check if this is MC ---------------------------------------
+    Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
+
+    // --- Add the task ----------------------------------------------
+    gROOT->Macro(Form("AddTaskFMDELoss.C(%d,%d)", mc, fUseCent));
+  }
+  /** 
+   * Create entrality selection if enabled 
+   * 
+   * @param mc   Whether this is MC or not
+   * @param mgr  Analysis manager 
+   */
+  virtual void CreateCentralitySelection(Bool_t mc, AliAnalysisManager* mgr)
+  {
+    if (!fUseCent) return;
+
+    gROOT->Macro("AddTaskCentrality.C");
+    AliCentralitySelectionTask* ctask = 
+      dynamic_cast<AliCentralitySelectionTask*>(mgr->GetTask("CentralitySelection"));
+    if (!ctask) return;
+    ctask->SetPass(fESDPass);
+    if (mc) ctask->SetMCInput();
+  }
+  /** 
+   * Crete output handler - we don't want one here. 
+   * 
+   * @return 0
+   */
+  AliVEventHandler* CreateOutputHandler(EType) { return 0; }
+  Bool_t fUseCent; // Whether to use centrality or not 
+};
+
+//
+// EOF
+//
diff --git a/PWG2/FORWARD/analysis2/trains/MakeMCCorrTrain.C b/PWG2/FORWARD/analysis2/trains/MakeMCCorrTrain.C
new file mode 100644 (file)
index 0000000..768b363
--- /dev/null
@@ -0,0 +1,144 @@
+#include "TrainSetup.C"
+
+//====================================================================
+/**
+ * Analysis train to make Forward and Central MC corrections
+ * 
+ * To run, do 
+ * @code 
+ * gROOT->LoadMacro("TrainSetup.C");
+ * // Make train 
+ * MakeMCCorrTrain t("My Analysis");
+ * // Set variaous parameters on the train 
+ * t.SetDataDir("/home/of/data");
+ * t.AddRun(118506)
+ * // Run it 
+ * t.Run("LOCAL", "FULL", -1, false, false);
+ * @endcode 
+ *
+ * @ingroup pwg2_forward_scripts_makers
+ * @ingroup pwg2_forward_mc
+ */
+class MakeMCCorrTrain : public TrainSetup
+{
+public:
+  /** 
+   * Constructor.  Date and time must be specified when running this
+   * in Termiante mode on Grid
+   * 
+   * @param name     Name of train (free form)
+   * @param dateTime Append date and time to name 
+   * @param year     Year     - if not specified, current year
+   * @param month    Month    - if not specified, current month
+   * @param day      Day      - if not specified, current day
+   * @param hour     Hour     - if not specified, current hour
+   * @param min      Minutes  - if not specified, current minutes
+   */
+  MakeMCCorrTrain(const  char* name, 
+                 Bool_t       dateTime = false,
+                 UShort_t     year     = 0, 
+                 UShort_t     month    = 0, 
+                 UShort_t     day      = 0, 
+                 UShort_t     hour     = 0, 
+                 UShort_t     min      = 0) 
+    : TrainSetup(name, dateTime, 
+                year, month, day, hour, min)
+  {}
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode - see TrainSetup::EMode
+   * @param oper     Operation - see TrainSetup::EOperation
+   * @param nEvents  Number of events (negative means all)
+   * @param usePar   If true, use PARs 
+   */
+  void Run(const char* mode, const char* oper, 
+          Int_t nEvents=-1, Bool_t usePar=false)
+  {
+    Info("Run", "Running in mode=%s, oper=%s, events=%d, Par=%d", 
+        mode, oper, nEvents, usePar);
+    Exec("ESD", mode, oper, nEvents, true, usePar);
+  }
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode - see TrainSetup::EMode
+   * @param oper     Operation - see TrainSetup::EOperation
+   * @param nEvents  Number of events (negative means all)
+   * @param usePar   If true, use PARs 
+   */
+  void Run(EMode mode, EOper oper, Int_t nEvents=-1,
+          Bool_t usePar = false)
+  {
+    Info("Run", "Running in mode=%d, oper=%d, events=%d, Par=%d", 
+        mode, oper, nEvents, usePar);
+    Exec(kESD, mode, oper, nEvents, true, usePar);
+  }
+protected:
+  /** 
+   * Create the tasks 
+   * 
+   * @param mode Processing mode
+   * @param par  Whether to use par files 
+   * @param mgr  Analysis manager 
+   */
+  void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)
+  {
+    // --- Output file name ------------------------------------------
+    AliAnalysisManager::SetCommonFileName("forward_mccorr.root");
+
+    // --- Load libraries/pars ---------------------------------------
+    LoadLibrary("PWG2forward2", mode, par, true);
+    
+    // --- Set load path ---------------------------------------------
+    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+                            gROOT->GetMacroPath()));
+
+    // --- Check if this is MC ---------------------------------------
+    if (!mgr->GetMCtruthEventHandler()) return;
+    
+    // --- Add the task ----------------------------------------------
+    gROOT->Macro("AddTaskForwardMCCorr.C"); 
+
+    // --- Add the task ----------------------------------------------
+    gROOT->Macro("AddTaskCentralMCCorr.C");
+  }
+  //__________________________________________________________________
+  /** 
+   * Create physics selection , and add to manager
+   * 
+   * @param mc Whether this is for MC 
+   * @param mgr Manager 
+   */
+  void CreatePhysicsSelection(Bool_t mc,
+                             AliAnalysisManager* mgr)
+  {
+    TrainSetup::CreatePhysicsSelection(mc, mgr);
+
+    // --- Get input event handler -----------------------------------
+    AliInputEventHandler* ih =
+      dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
+    if (!ih) 
+      Fatal("CreatePhysicsSelection", "Couldn't get input handler (%p)", ih);
+    
+    // --- Get Physics selection -------------------------------------
+    AliPhysicsSelection* ps = 
+      dynamic_cast<AliPhysicsSelection*>(ih->GetEventSelection());
+    if (!ps) 
+      Fatal("CreatePhysicsSelection", "Couldn't get PhysicsSelection (%p)", ps);
+
+    // --- Ignore trigger class when selecting events.  This means ---
+    // --- that we get offline+(A,C,E) events too --------------------
+    // ps->SetSkipTriggerClassSelection(true);
+  }
+  /** 
+   * Do not the centrality selection
+   */
+  void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {}
+  Double_t fVzMin;     // Least v_z
+  Double_t fVzMax;     // Largest v_z
+};
+
+//
+// EOF
+//
diff --git a/PWG2/FORWARD/analysis2/trains/MakedNdetaTrain.C b/PWG2/FORWARD/analysis2/trains/MakedNdetaTrain.C
new file mode 100644 (file)
index 0000000..8c527ed
--- /dev/null
@@ -0,0 +1,159 @@
+#include "TrainSetup.C"
+
+//====================================================================
+/**
+ * Analysis train to make @f$ dN/d\eta@f$
+ * 
+ * To run, do 
+ * @code 
+ * gROOT->LoadMacro("TrainSetup.C");
+ * // Make train 
+ * MakedNdetaTrain t("My Analysis");
+ * // Set variaous parameters on the train 
+ * t.SetDataDir("/home/of/data");
+ * t.AddRun(118506)
+ * // Run it 
+ * t.Run("LOCAL", "FULL", -1, false, false);
+ * @endcode 
+ *
+ * @ingroup pwg2_forward_scripts_makers
+ * @ingroup pwg2_forward_dndeta
+ */
+class MakedNdetaTrain : public TrainSetup
+{
+public:
+  /** 
+   * Constructor.  Date and time must be specified when running this
+   * in Termiante mode on Grid
+   * 
+   * @param name     Name of train (free form)
+   * @param trig     Trigger to use 
+   * @param vzMin    Least @f$ v_z@f$
+   * @param vzMax    Largest @f$ v_z@f$
+   * @param scheme   Normalisation scheme 
+   * @param useCent  Whether to use centrality 
+   * @param dateTime Append date and time to name 
+   * @param year     Year     - if not specified, current year
+   * @param month    Month    - if not specified, current month
+   * @param day      Day      - if not specified, current day
+   * @param hour     Hour     - if not specified, current hour
+   * @param min      Minutes  - if not specified, current minutes
+   */
+  MakedNdetaTrain(const char* name, 
+                 const char* trig="INEL", 
+                 Double_t    vzMin=-10, 
+                 Double_t    vzMax=10, 
+                 const char* scheme="FULL", 
+                 Bool_t      useCent=false,
+                 Bool_t      dateTime=false,
+                 UShort_t    year  = 0, 
+                 UShort_t    month = 0, 
+                 UShort_t    day   = 0, 
+                 UShort_t    hour  = 0, 
+                 UShort_t    min   = 0) 
+    : TrainSetup(name, dateTime, year, month, day, hour, min),
+      fTrig(trig), 
+      fVzMin(vzMin), 
+      fVzMax(vzMax),
+      fScheme(scheme),
+      fUseCent(useCent)
+  {}
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode - see TrainSetup::EMode
+   * @param oper     Operation - see TrainSetup::EOperation
+   * @param nEvents  Number of events (negative means all)
+   * @param usePar   If true, use PARs 
+   */
+  void Run(const char* mode, const char* oper, 
+          Int_t nEvents=-1, Bool_t usePar=false)
+  {
+    Exec("AOD", mode, oper, nEvents, false, usePar);
+  }
+  /** 
+   * Run this analysis 
+   * 
+   * @param mode     Mode - see TrainSetup::EMode
+   * @param oper     Operation - see TrainSetup::EOperation
+   * @param nEvents  Number of events (negative means all)
+   * @param usePar   If true, use PARs 
+   */
+  void Run(EMode mode, EOper oper, Int_t nEvents=-1, 
+          Bool_t usePar=false)
+  {
+    Exec(kAOD, mode, oper, nEvents, false, usePar);
+  }
+  /** 
+   * Set the trigger to use (INEL, INEL>0, NSD)
+   * 
+   * @param trig Trigger to use 
+   */
+  void SetTrigger(const char* trig) { fTrig = trig; }
+  /** 
+   * Set the vertex range to accept 
+   * 
+   * @param min Minimum 
+   * @param max Maximum 
+   */
+  void SetVertexRange(Double_t min, Double_t max) { fVzMin=min; fVzMax=max; }
+  /** 
+   * Set the normalisation scheme 
+   * 
+   * @param scheme Normalisation scheme options 
+   */
+  void SetScheme(const char* scheme) { fScheme = scheme; }
+  /** 
+   * Whether to use centrality or not 
+   * 
+   * @param use To use the centrality 
+   */
+  void SetUseCentrality(Bool_t use) { fUseCent = use; }
+protected:
+  /** 
+   * Create the tasks 
+   * 
+   * @param mode Processing mode
+   * @param par  Whether to use par files 
+   */
+  void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager*)
+  {
+    // --- Output file name ------------------------------------------
+    AliAnalysisManager::SetCommonFileName("forward_dndeta.root");
+
+    // --- Load libraries/pars ---------------------------------------
+    LoadLibrary("PWG2forward2", mode, par, true);
+    
+    // --- Set load path ---------------------------------------------
+    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+                            gROOT->GetMacroPath()));
+
+    // --- Add the task ----------------------------------------------
+    gROOT->Macro(Form("AddTaskForwarddNdeta.C(\"%s\",%f,%f,%d,\"%s\")",
+                     fTrig.Data(), fVzMin, fVzMax, fUseCent, fScheme.Data()));
+
+    gROOT->Macro(Form("AddTaskCentraldNdeta.C(\"%s\",%f,%f,%d,\"%s\")",
+                     fTrig.Data(), fVzMin, fVzMax, fUseCent, fScheme.Data()));
+
+    gROOT->Macro(Form("AddTaskMCTruthdNdeta.C(\"%s\",%f,%f,%d,\"%s\")",
+                     fTrig.Data(), fVzMin, fVzMax, fUseCent, fScheme.Data()));
+  }
+  /** 
+   * Do not the centrality selection
+   */
+  void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {}
+  /** 
+   * Crete output handler - we don't want one here. 
+   * 
+   * @return 0
+   */
+  AliVEventHandler* CreateOutputHandler(EType) { return 0; }
+  TString  fTrig;      // Trigger to use 
+  Double_t fVzMin;     // Least v_z
+  Double_t fVzMax;     // Largest v_z
+  TString  fScheme;    // Normalisation scheme 
+  Bool_t   fUseCent;   // Use centrality 
+};
+//
+// EOF
+//
diff --git a/PWG2/FORWARD/analysis2/trains/TrainSetup.C b/PWG2/FORWARD/analysis2/trains/TrainSetup.C
new file mode 100644 (file)
index 0000000..86af4df
--- /dev/null
@@ -0,0 +1,1579 @@
+/**
+ * @file   TrainSetup.C
+ * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
+ * @date   Wed Mar 23 12:12:00 2011
+ * 
+ * @brief  
+ * 
+ * @ingroup pwg2_forward_scripts_makers
+ * 
+ */
+
+#ifndef __CINT__
+#include <fstream>
+#include <iostream>
+
+#include <TAlienCollection.h>
+#include <TArrayI.h>
+#include <TChain.h>
+#include <TDatime.h>
+#include <TEnv.h>
+#include <TFile.h>
+#include <TGrid.h>
+#include <TList.h>
+#include <TObjString.h>
+#include <TProof.h>
+#include <TString.h>
+#include <TSystem.h>
+#include <TSystemDirectory.h>
+#include <TSystemFile.h>
+#include <TROOT.h>
+
+#include <AliAODHandler.h>
+#include <AliAODInputHandler.h>
+#include <AliAnalysisDataContainer.h>
+#include <AliAnalysisManager.h>
+#include <AliAnalysisAlien.h>
+#include <AliESDInputHandler.h>
+#include <AliMCEventHandler.h>
+#include <AliVEventHandler.h>
+#include <AliPhysicsSelection.h>
+#include <AliCentralitySelectionTask.h>
+#else
+class TArrayI;
+class TChain;
+class AliAnalysisManager;
+#endif
+
+//====================================================================
+/** 
+ * Generic set-up of an analysis train using the grid-handler (AliEn plugin). 
+ * 
+ * Users should define a class that derives from this.  The class
+ * should implement the member function CreateTasks to add needed
+ * tasks to the train
+ * 
+ * @code 
+ * // MyTrain.C 
+ * class MyTrain : public TrainSetup
+ * {
+ * public:
+ *   MyTrain(Bool_t   dateTime = false, 
+ *           UShort_t year     = 0, 
+ *           UShort_t month    = 0, 
+ *           UShort_t day      = 0, 
+ *           UShort_t hour     = 0, 
+ *           UShort_t min      = 0) 
+ *     : TrainSetup("My train", dateTime, year, month, day, hour, min)
+ *   {}
+ *   void Run(const char* type, const char* mode, const char* oper, 
+ *            Int_t nEvents=-1, Bool_t mc=false,
+ *            Bool_t usePar=false)
+ *   {
+ *     Exec(type, mode, oper, nEvents, mc, usePar);
+ *   }
+ * protected:
+ *   void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)
+ *   {
+ *     AliAnalysisManager::SetCommonFileName("my_analysis.root");
+ *     LoadLibrary("MyAnalysis", mode, par, true);
+ *     Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
+ *     gROOT->Macro("MyAnalysis.C");
+ *   }
+ * };
+ * @endcode 
+ * 
+ * This can then be run like 
+ * 
+ * @verbatim 
+ * > aliroot 
+ * Root> .L TrainSetup.C 
+ * Root> .L MyTrain.C 
+ * Root> MyTrain t;
+ * Root> t.Run();
+ * @endverbatim 
+ * 
+ * or as a script 
+ * 
+ * @code 
+ * {
+ *   gROOT->LoadMacro("TrainSetup.C");
+ *   gROOT->LoadMacro("MyTrain.C");
+ *   MyTrain t;
+ *   t.Run();
+ * }
+ * @endcode 
+ * 
+ * To byte compile this, you need to 
+ * - load the ROOT AliEn library
+ * - load the analysis libraries 
+ * - add $ALICE_ROOT/include to header search 
+ * first 
+ *
+ * @verbatim 
+ * > aliroot 
+ * Root> gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
+ * Root>                          "$ALICE_ROOT/ANALYSIS/macros",
+ * Root>                         gROOT->GetMacroPath()));
+ * Root> gSystem->AddIncludePath("-I${ALICE_ROOT}/include");
+ * Root> gSystem->Load("libRAliEn");
+ * Root> gSystem->Load("libANALYSIS");
+ * Root> gSystem->Load("libANALYSISalice");
+ * Root> gROOT->LoadMacro("TrainSetup.C+");
+ * @endverbatim 
+ * 
+ * @ingroup pwg2_forward_scripts_makers
+ * 
+ */
+struct TrainSetup
+{
+  /** 
+   * Data type to process 
+   */
+  enum EType { 
+    /** Event Summary Data */
+    kESD, 
+    /** Analysis Object Data */
+    kAOD
+  };
+  /**
+   * How to run the analysis
+   * 
+   */
+  enum EMode {
+    /** Locally */ 
+    kLocal = 1, 
+    /** In PROOF(-Lite) cluster */
+    kProof, 
+    /** On the Grid */
+    kGrid 
+  };
+  /**
+   * What stage of the analysis to run 
+   * 
+   */
+  enum EOper { 
+    /** Testing.  Local processing, a single copied from Grid */
+    kTest, 
+    /** Off-line */
+    kOffline, 
+    /** Submit to queue */
+    kSubmit, 
+    /** Merge and terminate */
+    kTerminate, 
+    /** Full run */
+    kFull
+  };
+
+  //__________________________________________________________________
+  /** 
+   * Constructor 
+   * 
+   * @param name         Name of analysis (free-form)
+   * @param useDateTime  Whether to append date and time to the name 
+   * @param year         Year - if not specified, taken from current date
+   * @param month        Month - if not specified, taken from current date 
+   * @param day          Day - if not specified, taken from current date 
+   * @param hour         Hour - if not specified, taken from current time  
+   * @param min          Minute - if not specified, taken from current time  
+   */
+  TrainSetup(const char* name, Bool_t useDateTime=true, 
+            UShort_t year=0, UShort_t month=0, 
+            UShort_t day=0, UShort_t hour=0, UShort_t min=0) 
+    : fName(name),
+      fRootVersion("v5-28-00a"),
+      fAliRootVersion("v4-21-18-AN"),
+      fAliEnAPIVersion("V1.1x"),
+      fProofServer("alicecaf.cern.ch"),
+      fDataDir("/alice/data/2010/LHC10c"),
+      fDataSet("/COMMON/COMMON/LHC09a4_run8100X#/esdTree"),
+      fXML(""), 
+      fRunNumbers(0),
+      fListOfPARs(),
+      fListOfSources(),
+      fListOfLibraries(),
+      fListOfExtras(),
+      fNReplica(4),
+      fESDPass(3),
+      fPassPostfix(""),
+      fEscapedName(name),
+      fAllowOverwrite(kFALSE)
+  {
+    char  c[] = { ' ', '/', '@', 0 };
+    char* p   = c;
+    while (*p) { 
+      fEscapedName.ReplaceAll(Form("%c", *p), "_");
+      p++;
+    }
+
+    if (useDateTime) { 
+      if (year == 0 || month == 0 || day == 0) {
+       TDatime now;
+       year  = now.GetYear();
+       month = now.GetMonth();
+       day   = now.GetDay();
+       hour  = now.GetHour();
+       min   = now.GetMinute();
+      }
+      fEscapedName.Append(Form("_%04d%02d%02d_%02d%02d", 
+                              year, month, day, hour, min));
+    }
+
+  }
+
+  //__________________________________________________________________
+  /** 
+   * Parse a string into a type enum
+   * 
+   * @param type String to pass
+   * 
+   * @return Enumaration value 
+   */
+  static EType ParseType(const char* type, Bool_t& /*mc*/)
+  {
+    // mc = false;
+    TString sType(type);
+    sType.ToUpper();
+    EType eType = kESD;
+    // if      (sType.Contains("MC"))    mc    = true;
+    if      (sType.Contains("ESD"))   eType = kESD; 
+    else if (sType.Contains("AOD"))   eType = kAOD;
+    else 
+      Fatal("Run", "Unknown type '%s'", type);
+    
+    return eType;
+  }
+  //__________________________________________________________________
+  /** 
+   * Return a string that reflects the passed mode
+   * 
+   * @param eMode Mode 
+   * 
+   * @return String representation of mode 
+   */
+  static const char* ModeString(EMode eMode) 
+  {
+    switch (eMode) {
+    case kLocal:       return "LOCAL";
+    case kProof:       return "PROOF";
+    case kGrid:                return "GRID";
+    }
+    return 0;
+  }
+  //__________________________________________________________________
+  /** 
+   * Parse a string for mode specifier 
+   * 
+   * @param mode Mode string
+   * 
+   * @return EMode value
+   */
+  static EMode ParseMode(const char* mode)
+  {
+    TString sMode(mode);
+    sMode.ToUpper();
+    EMode eMode = kLocal;
+    if      (sMode == "LOCAL") eMode = kLocal;
+    else if (sMode == "PROOF") eMode = kProof;
+    else if (sMode == "GRID")  eMode = kGrid;
+    else 
+      Fatal("Run", "Unknown mode '%s'", mode);
+    return eMode;
+  }
+
+  //__________________________________________________________________
+  /** 
+   * Return a string that reflects the passed operation
+   * 
+   * @param eOper Operation
+   * 
+   * @return String representation of operation 
+   */
+  static const char* OperString(EOper eOper) 
+  {
+    switch (eOper) {
+    case kTest:                return "TEST";
+    case kOffline:     return "OFFLINE";
+    case kSubmit:      return "SUBMIT";
+    case kTerminate:   return "TERMINATE";
+    case kFull:                return "FULL";
+    }
+    return 0;
+  }
+  //__________________________________________________________________
+  /** 
+   * Parse an operation string 
+   * 
+   * @param oper Operation 
+   * 
+   * @return An EOper value
+   */
+  static EOper ParseOperation(const char* oper)
+  {
+    TString sOper(oper);
+    sOper.ToUpper();
+    EOper eOper = kFull;
+    if      (sOper == "TEST")      eOper = kTest;
+    else if (sOper == "OFFLINE")   eOper = kOffline;
+    else if (sOper == "SUBMIT")    eOper = kSubmit;
+    else if (sOper == "TERMINATE") eOper = kTerminate;
+    else if (sOper == "FULL")      eOper = kFull;
+    else 
+      Fatal("Run", "unknown operation '%s'", oper);
+    return eOper;
+  }
+
+  //__________________________________________________________________
+  /** 
+   * Set ROOT version to use 
+   * 
+   * @param v Version string of ROOT 
+   */
+  void SetROOTVersion(const char* v)    { fRootVersion = v; }
+  //__________________________________________________________________
+  /** 
+   * Set AliROOT version to use 
+   * 
+   * @param v Version string of AliROOT 
+   */
+  void SetAliROOTVersion(const char* v) { fAliRootVersion = v; }
+  //__________________________________________________________________
+  /** 
+   * Set the PROOF server URL
+   * 
+   * @param s PROOF server URL 
+   */
+  void SetProofServer(const char* s)    { fProofServer = s; }
+  //__________________________________________________________________
+  /** 
+   * Set the GRID/Local data dir 
+   * 
+   * @param d Directory with data 
+   */
+  void SetDataDir(const char* d) { fDataDir = d; }
+  //__________________________________________________________________
+  /** 
+   * Set the PROOF data set 
+   * 
+   * @param d PROOF registered data set 
+   */
+  void SetDataSet(const char* d) { fDataSet = d; }
+  //__________________________________________________________________
+  /** 
+   * Set the XML file to use 
+   * 
+   * @param x XML file 
+   */
+  void SetXML(const char* x) { fXML = x; }
+  //__________________________________________________________________
+  /** 
+   * Set how many replicas of the output we want 
+   * 
+   * @param n Number of replicas requested 
+   */
+  void SetNReplica(Int_t n) { fNReplica = n; }
+  /** 
+   * Set the ESD pass to use 
+   * 
+   * @param pass Pass number 
+   */
+  void SetESDPass(Int_t pass) { fESDPass = pass; }
+  /** 
+   * Set the ESD pass to use 
+   * 
+   * @param postfix Post fix to pass number 
+   */
+  void SetPassPostfix(const char* postfix) { fPassPostfix = postfix; }
+  //__________________________________________________________________
+  /** 
+   * Add a source file to be copied and byte compiled on slaves 
+   * 
+   * @param src          Sources 
+   * @param addToExtra   If false, do not copy 
+   */
+  void AddSource(const char* src, bool addToExtra=true) 
+  { 
+    fListOfSources.Add(new TObjString(src)); 
+    if (addToExtra) AddExtraFile(src); // Source code isn't copied!
+  }
+  //__________________________________________________________________
+  /** 
+   * Add binary data to be uploaded to slaves 
+   * 
+   * @param lib Name of binary file 
+   */
+  void AddLibrary(const char* lib) { fListOfLibraries.Add(new TObjString(lib));}
+  //__________________________________________________________________
+  /** 
+   * Add a run to be analysed
+   *  
+   * @param run Run number
+   */
+  void AddRun(Int_t run) 
+  {
+    Int_t i = fRunNumbers.fN; fRunNumbers.Set(i+1); fRunNumbers[i] = run;
+  }
+  //__________________________________________________________________
+  /** 
+   * Read run numbers from a file 
+   * 
+   * @param filename File name 
+   */
+  void ReadRunNumbers(const char* filename)
+  {
+    std::ifstream file(filename);
+    if (!file) 
+      Fatal("ReadRunNumbers", "Cannot read from %s", filename);
+    
+    while (!file.eof()) { 
+      Int_t run;
+      file >> run;
+      AddRun(run);
+      Char_t c;
+      file >> c;
+      if (file.bad()) break;
+    }
+    file.close();
+  }
+  //__________________________________________________________________
+  /** 
+   * Add an extra file to be uploaded to slave 
+   * 
+   * @param file Extra file to be uploaded 
+   */
+  void AddExtraFile(const char* file)
+  {
+    if (!file || file[0] == '\0') return;
+    fListOfExtras.Add(new TObjString(file));
+  }
+  //__________________________________________________________________
+  /** 
+   * Set whether to allow overwritting existing files/directories 
+   * 
+   * @param allow If true, allow overwritting files/directories
+   */
+  void SetAllowOverwrite(Bool_t allow) { fAllowOverwrite = allow; }
+  //__________________________________________________________________
+  /** 
+   * Print the setup 
+   * 
+   */
+  void Print() const 
+  {
+    bool mc = AliAnalysisManager::GetAnalysisManager()
+      ->GetMCtruthEventHandler();
+    std::cout << fName << " train setup\n"
+             << std::boolalpha
+             << "  ROOT version:         " << fRootVersion    << "\n"
+             << "  AliROOT version:      " << fAliRootVersion << "\n"
+             << "  Name of proof server: " << fProofServer    << "\n"
+             << "  Grid Input directory: " << fDataDir        << "\n"
+             << "  Proof data set name:  " << fDataSet        << "\n"
+             << "  XML collection:       " << fXML            << "\n"
+             << "  Monte-Carlo input:    " << mc              << "\n"
+             << "  Storage replication:  " << fNReplica       << "\n"
+             << "  Run numbers:          " << std::flush;
+    for (Int_t i = 0; i < fRunNumbers.GetSize(); i++) 
+      std::cout << (i == 0 ? "" : ", ") << fRunNumbers.At(i);
+
+    std::cout << "\n"
+             << "  PAR files:            " << std::flush;
+    Bool_t first = true;
+    TObject* obj = 0;
+    TIter nextPar(&fListOfPARs);
+    while ((obj = nextPar())) {
+      std::cout << (first ? "" : ", ") << obj->GetName();
+      first = false;
+    }
+
+    std::cout << "\n"
+             << "  Script sources:       " << std::flush;
+    first = true;
+    TIter nextSrc(&fListOfSources);
+    while ((obj = nextSrc())) {
+      std::cout << (first ? "" : ", ") << obj->GetName();
+      first = false;
+    }
+
+    std::cout << "\n"
+             << "  Libraries to load:    " << std::flush;
+    first = true;
+    TIter nextLib(&fListOfLibraries);
+    while ((obj = nextLib())) {
+      std::cout << (first ? "" : ", ") << obj->GetName();
+      first = false;
+    }
+    std::cout << std::noboolalpha << std::endl;
+
+    AliAnalysisGrid* plugin = 
+      AliAnalysisManager::GetAnalysisManager()->GetGridHandler();
+    if (!plugin) return;
+    
+  }
+
+protected:
+  //__________________________________________________________________
+  /** 
+   * Copy constructor 
+   * 
+   * @param o Object to copy from
+   */
+  TrainSetup(const TrainSetup& o)
+    : fName(o.fName),
+      fRootVersion(o.fRootVersion),
+      fAliRootVersion(o.fAliRootVersion),
+      fProofServer(o.fProofServer),
+      fDataDir(o.fDataDir),    
+      fDataSet(o.fDataSet),    
+      fXML(o.fXML),    
+      fRunNumbers(o.fRunNumbers),
+      fListOfPARs(),
+      fListOfSources(),
+      fListOfLibraries(),
+      fListOfExtras(),
+      fNReplica(o.fNReplica),
+      fESDPass(o.fESDPass)
+  {
+    if (isdigit(fName[0])) { 
+      Warning("TrainSetup", "Name starts with a digit, prepending 'a' to name");
+      fName = Form("a%s", fName.Data());
+    }
+    TObject* obj = 0;
+    TIter nextPar(&o.fListOfPARs);
+    while ((obj = nextPar())) fListOfPARs.Add(obj->Clone());
+    TIter nextSrc(&o.fListOfSources);
+    while ((obj = nextSrc())) fListOfSources.Add(obj->Clone());
+    TIter nextLib(&o.fListOfLibraries);
+    while ((obj = nextLib())) fListOfLibraries.Add(obj->Clone());
+    TIter nextExa(&o.fListOfExtras);
+    while ((obj = nextExa())) fListOfExtras.Add(obj->Clone());
+  }
+  //__________________________________________________________________
+  /** 
+   * Assignment operator 
+   * 
+   * @param o Object to assign from 
+   * 
+   * @return Reference to this object. 
+   */
+  TrainSetup& operator=(const TrainSetup& o)
+  {
+    fName              = o.fName;
+    fRootVersion       = o.fRootVersion;
+    fAliRootVersion    = o.fAliRootVersion;
+    fProofServer       = o.fProofServer;
+    fDataDir           = o.fDataDir;   
+    fDataSet           = o.fDataSet;   
+    fXML               = o.fXML;       
+    fNReplica          = o.fNReplica;  
+    fESDPass            = o.fESDPass;
+    fRunNumbers         = o.fRunNumbers;
+    TObject* obj = 0;
+    TIter nextPar(&o.fListOfPARs);
+    while ((obj = nextPar())) fListOfPARs.Add(obj->Clone());
+    TIter nextSrc(&o.fListOfSources);
+    while ((obj = nextSrc())) fListOfSources.Add(obj->Clone());
+    TIter nextLib(&o.fListOfLibraries);
+    while ((obj = nextLib())) fListOfLibraries.Add(obj->Clone());
+    TIter nextExa(&o.fListOfExtras);
+    while ((obj = nextExa())) fListOfExtras.Add(obj->Clone());
+
+    return *this;
+  }
+
+  //__________________________________________________________________
+  /** 
+   * Run this analysis 
+   * 
+   * @param type    Type of input for analysis  (kESD, kAOD)
+   * @param mode    Mode of job (kLocal, kProof, kGrid)
+   * @param oper    Operation 
+   * @param nEvents Number of events to analyse (<0 means all)
+   * @param mc      Whether to connect MC data 
+   * @param usePar  Whether to use PARs  
+   * @param dbg     Debug level
+   */
+  void Exec(const char*  type, 
+           const char*  mode="GRID", 
+           const char*  oper="FULL", 
+           Int_t        nEvents=-1, 
+           Bool_t       mc=false, 
+           Bool_t       usePar=false, 
+           Int_t        dbg=0)
+  {
+    Info("Exec", "Doing exec with type=%s, mode=%s, oper=%s, events=%d "
+        "mc=%d, usePar=%d", type, mode, oper, nEvents, mc, usePar);
+    EType eType = ParseType(type, mc);
+    EMode eMode = ParseMode(mode);
+    EOper eOper = ParseOperation(oper);
+
+    Exec(eType, eMode, eOper, nEvents, mc, usePar, dbg);
+  }
+
+  //__________________________________________________________________
+  /** 
+   * Run this analysis 
+   * 
+   * @param type    Type of input for analysis  (kESD, kAOD)
+   * @param mode    Mode of job (kLocal, kProof, kGrid)
+   * @param oper    Operation 
+   * @param nEvents Number of events to analyse (<0 means all)
+   * @param mc      Whether to connect MC data 
+   * @param usePar  Whether to use PARs  
+   * @param dbg     Debug level
+   */
+  void Exec(EType  type, 
+           EMode  mode, 
+           EOper  oper, 
+           Int_t  nEvents, 
+           Bool_t mc, 
+           Bool_t usePar, 
+           Int_t  dbg=0)
+  {
+    Info("Exec", "Doing exec with type=%d, mode=%d, oper=%d, events=%d "
+        "mc=%d, usePar=%d", type, mode, oper, nEvents, mc, usePar);
+
+    if (mode == kProof) usePar    = true;
+
+    if (!Connect(mode)) return;
+
+    TString cwd = gSystem->WorkingDirectory();
+    TString nam = EscapedName();
+    if (oper != kTerminate) { 
+      if (!fAllowOverwrite && !gSystem->AccessPathName(nam.Data())) {
+       Error("Exec", "File/directory %s already exists", nam.Data());
+       return;
+      }
+      if (gSystem->AccessPathName(nam.Data())) {
+       if (gSystem->MakeDirectory(nam.Data())) {
+         Error("Exec", "Failed to make directory %s", nam.Data());
+         return;
+       }
+      }
+    }
+    else {
+      if (gSystem->AccessPathName(nam.Data())) {
+       Error("Exec", "File/directory %s does not exists", nam.Data());
+       return;
+      }
+    }
+      
+    if (!gSystem->ChangeDirectory(nam.Data())) { 
+      Error("Exec", "Failed to change directory to %s", nam.Data());
+      return;
+    }
+    Info("Exec", "Made subdirectory %s, and cd'ed there", nam.Data());
+      
+    if (!LoadCommonLibraries(mode, usePar)) return;
+    
+    // --- Create analysis manager -----------------------------------
+    AliAnalysisManager *mgr  = new AliAnalysisManager(fName,"Analysis Train");
+
+    // In test mode, collect system information on every event 
+    // if (oper == kTest)  mgr->SetNSysInfo(1); 
+    if (dbg  >  0)      mgr->SetDebugLevel(dbg);
+    if (mode == kLocal) mgr->SetUseProgressBar(kTRUE, 100);
+   
+    // --- ESD input handler ------------------------------------------
+    AliVEventHandler*  inputHandler = CreateInputHandler(type);
+    if (inputHandler) mgr->SetInputEventHandler(inputHandler);
+    
+    // --- Monte-Carlo ------------------------------------------------
+    AliVEventHandler*  mcHandler = CreateMCHandler(type,mc);
+    if (mcHandler) mgr->SetMCtruthEventHandler(mcHandler);
+    
+    // --- AOD output handler -----------------------------------------
+    AliVEventHandler*  outputHandler = CreateOutputHandler(type);
+    if (outputHandler) mgr->SetOutputEventHandler(outputHandler);
+    
+    // --- Include analysis macro path in search path ----------------
+    gROOT->SetMacroPath(Form("%s:%s:$ALICE_ROOT/ANALYSIS/macros",
+                            cwd.Data(), gROOT->GetMacroPath()));
+
+    // --- Physics selction - only for ESD ---------------------------
+    if (type == kESD) CreatePhysicsSelection(mc, mgr);
+    
+    // --- Create centrality task ------------------------------------
+    CreateCentralitySelection(mc, mgr);
+
+    // --- Create tasks ----------------------------------------------
+    CreateTasks(mode, usePar, mgr);
+
+    // --- Create Grid handler ----------------------------------------
+    // _must_ be done after all tasks has been added
+    AliAnalysisAlien* gridHandler = CreateGridHandler(type, mode, oper);
+    if (gridHandler) mgr->SetGridHandler(gridHandler);
+    
+    // --- Create the chain ------------------------------------------
+    TChain* chain = CreateChain(type, mode, oper, mc);
+    if (mode == kLocal && !chain) {
+      Error("Exec", "No chain defined in local mode!");
+      return;
+    }
+
+    // --- Print setup -----------------------------------------------
+    Print();
+
+    // --- Initialise the train --------------------------------------
+    if (!mgr->InitAnalysis())  {
+      gSystem->ChangeDirectory(cwd.Data());
+      Error("Run","Failed to initialise train");
+      return;
+    }
+
+    // --- Show status -----------------------------------------------
+    mgr->PrintStatus();
+
+    Long64_t ret = StartAnalysis(mgr, mode, chain, nEvents);
+
+    // Make sure we go back 
+    gSystem->ChangeDirectory(cwd.Data());
+
+    if (ret < 0) Error("Exec", "Analysis failed");
+  }
+  //__________________________________________________________________
+  /** 
+   * Start the analysis 
+   * 
+   * @param mgr       Analysis manager
+   * @param mode      Run mode
+   * @param chain     Input data (local and proof only)
+   * @param nEvents   Number of events to analyse 
+   */
+  Long64_t StartAnalysis(AliAnalysisManager* mgr, 
+                        EMode               mode, 
+                        TChain*             chain,
+                        Int_t               nEvents)
+  {
+    // --- Run the analysis ------------------------------------------
+    switch (mode) { 
+    case kLocal: 
+      if (!chain) {
+       Error("StartAnalysis", "No chain defined");
+       return -1;
+      }
+      if (nEvents < 0) nEvents = chain->GetEntries();
+      return mgr->StartAnalysis(ModeString(mode), chain, nEvents);
+    case kProof: 
+      if (fDataSet.IsNull()) {
+       if (!chain) { 
+         Error("StartAnalysis", "No chain defined");
+         return -1;
+       }
+       if (nEvents < 0) nEvents = chain->GetEntries();
+       return mgr->StartAnalysis(ModeString(mode), chain, nEvents);
+      }
+      return mgr->StartAnalysis(ModeString(mode), fDataSet);
+    case kGrid: 
+      if (nEvents < 0)
+       return mgr->StartAnalysis(ModeString(mode));
+      return mgr->StartAnalysis(ModeString(mode), nEvents);
+    }
+    // We should never get  here 
+    return -1;
+  }
+  //__________________________________________________________________
+  /** 
+   * Return the escaped name 
+   * 
+   * 
+   * @return Escaped name 
+   */
+  const TString& EscapedName() const 
+  {
+    return fEscapedName;
+  }
+  //__________________________________________________________________
+  /** 
+   * Create a grid handler 
+   * 
+   * @param type Data type
+   * @param mode Run mode 
+   * @param oper Operation 
+   * 
+   * @return Grid handler 
+   */
+  virtual AliAnalysisAlien* 
+  CreateGridHandler(EType type, EMode mode, EOper oper)
+  {
+    if (mode != kGrid) return 0;
+
+    TString name = EscapedName();
+
+    // Create the plug-in object, and set run mode 
+    AliAnalysisAlien* plugin = new AliAnalysisAlien();
+    plugin->SetRunMode(OperString(oper));
+    
+    // Production mode - not used here 
+    // plugin->SetProductionMode();
+    
+    // Set output to be per run 
+    plugin->SetOutputToRunNo();
+
+    // Set the job tag 
+    plugin->SetJobTag(fName);
+
+    // Set number of test files - used in test mode only 
+    plugin->SetNtestFiles(1);
+    
+    // Set required version of software 
+    plugin->SetAPIVersion(fAliEnAPIVersion);
+    plugin->SetROOTVersion(fRootVersion);
+    plugin->SetAliROOTVersion(fAliRootVersion);
+
+    // Keep log files 
+    plugin->SetKeepLogs();
+
+    // Declare root of input data directory 
+    plugin->SetGridDataDir(fDataDir);
+
+    // Data search patterns 
+    TString pat;
+    if (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) {
+      pat = "*/";
+      plugin->SetRunPrefix("");
+    }
+    else {
+      pat = Form("*ESDs/pass%d%s/*/", fESDPass, fPassPostfix.Data());
+      plugin->SetRunPrefix("000");
+    }
+    pat.Append(Form("*%s.root", type == kESD ? "ESDs" : "AOD"));
+    plugin->SetDataPattern(pat);
+
+    // Add the run numbers 
+    for (Int_t i = 0; i < fRunNumbers.fN; i++) {
+      if (fRunNumbers[i] < 0) continue; 
+      plugin->AddRunNumber(fRunNumbers[i]);
+    }
+    
+    // Set the working directory to be the trains name (with special
+    // characters replaced by '_' and the date appended), and also set
+    // the output directory (relative to working directory)
+    plugin->SetGridWorkingDir(name.Data());
+    plugin->SetGridOutputDir("output");
+
+    // Enable configured PARs 
+    TIter nextPar(&fListOfPARs);
+    TObject* parName;
+    while ((parName = nextPar()))
+      plugin->EnablePackage(parName->GetName());
+    
+    // Add sources that need to be compiled on the workers using
+    // AcLIC. 
+    TString addSources = SetupSources();
+    if (!addSources.IsNull()) plugin->SetAnalysisSource(addSources.Data());
+
+    // Add binary libraries that should be uploaded to the workers 
+    TString addLibs = SetupLibraries();
+    if (!addLibs.IsNull()) plugin->SetAdditionalLibs(addLibs.Data());
+    
+    // Disable default outputs 
+    plugin->SetDefaultOutputs(true);
+
+    // Merge parameters 
+    plugin->SetMaxMergeFiles(20);
+    plugin->SetMergeExcludes("AliAOD.root "
+                           "*EventStat*.root "
+                           "*event_stat*.root");
+
+    // Set number of runs per master - set to one to per run
+    plugin->SetNrunsPerMaster(1);
+
+    // Loop over defined containers in the analysis manager, 
+    // and declare these as outputs 
+    TString listOfAODs  = "";
+    TString listOfHists = "";
+    AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+    AliAnalysisDataContainer* cont = 0;
+    TIter nextCont(mgr->GetOutputs());
+    while ((cont = static_cast<AliAnalysisDataContainer*>(nextCont()))) {
+      TString outName(cont->GetFileName());
+      TString& list = (outName == "default" ? listOfAODs : listOfHists);
+      if (outName == "default") { 
+       if (!mgr->GetOutputEventHandler()) continue; 
+
+       outName = mgr->GetOutputEventHandler()->GetOutputFileName();
+      }
+      if (list.Contains(outName)) continue;
+      if (!list.IsNull()) list.Append(",");
+      list.Append(outName);
+    }
+    if (!mgr->GetExtraFiles().IsNull()) { 
+      if (!listOfAODs.IsNull()) listOfAODs.Append("+");
+      TString extra = mgr->GetExtraFiles();
+      extra.ReplaceAll(" ", ",");
+      listOfAODs.Append(extra);
+   }
+    TString outArchive = Form("stderr, stdout@disk=%d", fNReplica);
+    if (!listOfHists.IsNull()) 
+      outArchive.Append(Form(" hist_archive.zip:%s@disk=%d", 
+                            listOfHists.Data(), fNReplica));
+    if (!listOfAODs.IsNull()) 
+      outArchive.Append(Form(" aod_archive.zip:%s@disk=%d", 
+                            listOfAODs.Data(), fNReplica));
+    if (listOfAODs.IsNull() && listOfHists.IsNull()) 
+      Fatal("CreateGridHandler", "No outputs defined");
+    // Disabled for now 
+    // plugin->SetOutputArchive(outArchive);
+    
+    // Set name of generated analysis macro 
+    plugin->SetAnalysisMacro(Form("%s.C", name.Data()));
+    
+    // Maximum number of sub-jobs 
+    // plugin->SetSplitMaxInputFileNumber(25);
+    
+    // Set the Time-To-Live 
+    plugin->SetTTL(70000);
+    
+    // Re-submit failed jobs as long as the ratio of failed jobs is
+    // below this percentage. 
+    plugin->SetMasterResubmitThreshold(95);
+
+    // Set the input format
+    plugin->SetInputFormat("xml-single");
+
+    // Set the name of the generated jdl 
+    plugin->SetJDLName(Form("%s.jdl", name.Data()));
+
+    // Set the name of the generated executable 
+    plugin->SetExecutable(Form("%s.sh", name.Data()));
+    
+    // Set the job price !?
+    plugin->SetPrice(1);
+
+    // Set whether to merge via JDL 
+    plugin->SetMergeViaJDL(true);
+    
+    // Fast read otion 
+    plugin->SetFastReadOption(false);
+
+    // Whether to overwrite existing output 
+    plugin->SetOverwriteMode(true);
+
+    // Set the executable binary name and options 
+    plugin->SetExecutableCommand("aliroot -b -q -x");
+    
+    // Split by storage element - must be lower case!
+    plugin->SetSplitMode("se");
+
+    return plugin;
+  }
+  //__________________________________________________________________
+  /** 
+   * Create input handler 
+   * 
+   * @param type 
+   * 
+   * @return 
+   */
+  virtual AliVEventHandler* CreateInputHandler(EType type)
+  {
+    switch (type) {
+    case kESD: return new AliESDInputHandler(); 
+    case kAOD: return new AliAODInputHandler(); 
+    }
+    return 0;
+  }
+  //__________________________________________________________________
+  /** 
+   * Create input handler 
+   * 
+   * @param type  Run type (ESD or AOD)
+   * @param mc    Assume monte-carlo input 
+   * 
+   * @return 
+   */
+  virtual AliVEventHandler* CreateMCHandler(EType type, bool mc)
+  {
+    if (!mc)          return 0;
+    if (type != kESD) return 0;
+    Info("CreateMCHandler", "Making MC handler");
+    AliMCEventHandler* mcHandler = new AliMCEventHandler();
+    mcHandler->SetReadTR(true); 
+    return mcHandler;
+  }
+  //__________________________________________________________________
+  /** 
+   * Create output event handler 
+   * 
+   * @param type 
+   * 
+   * @return 
+   */
+  virtual AliVEventHandler* CreateOutputHandler(EType type)
+  {
+    AliAODHandler* ret = new AliAODHandler();
+    switch (type) { 
+    case kESD: 
+      ret->SetOutputFileName("AliAOD.root");
+      break;
+    case kAOD: 
+      ret->SetOutputFileName("AliAOD.pass2.root");
+      break;
+    }
+    return ret;
+  }
+  //__________________________________________________________________
+  /** 
+   * Create physics selection , and add to manager
+   * 
+   * @param mc Whether this is for MC 
+   * @param mgr Manager
+   */
+  virtual void CreatePhysicsSelection(Bool_t mc,
+                                     AliAnalysisManager* mgr)
+  {
+    gROOT->Macro(Form("AddTaskPhysicsSelection.C(%d)", mc));
+    mgr->RegisterExtraFile("event_stat.root");
+  }
+  //__________________________________________________________________
+  /** 
+   * Create physics selection , and add to manager
+   * 
+   * @param mc Whether this is for MC 
+   * @param mgr Manager
+   */
+  virtual void CreateCentralitySelection(Bool_t mc, AliAnalysisManager* mgr)
+  {
+    gROOT->Macro("AddTaskCentrality.C");
+    AliCentralitySelectionTask* ctask = 
+      dynamic_cast<AliCentralitySelectionTask*>(mgr->GetTask("CentralitySelection"));
+    if (!ctask) return;
+    ctask->SetPass(fESDPass);
+    if (mc) ctask->SetMCInput();
+  }
+  //__________________________________________________________________
+  /** 
+   * Create analysis tasks 
+   * 
+   * @param mode Run mode
+   * @param mgr  Manager
+   * @param par  Whether to use pars 
+   */
+  virtual void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager* mgr)=0;
+  //__________________________________________________________________
+  /** 
+   * Connect to external services (Proof and/or grid)
+   * 
+   * @param mode Running mode 
+   * 
+   * @return true on success 
+   */
+  virtual Bool_t Connect(EMode mode)
+  {
+    if (mode == kLocal) return true;
+                         
+    // --- Set-up connections to Proof cluster and alien -------------
+    if (mode == kProof) { 
+      // --- Find user name ------------------------------------------
+      TString userName(gSystem->Getenv("alien_API_USER"));
+      if (userName.IsNull()) {
+       userName = gSystem->GetUserInfo()->fUser;
+       Warning("Connect", 
+               "environment variable 'alien_API_USER' not set, using %s", 
+               userName.Data());
+      }
+
+      // --- Set prefered GSI method ---------------------------------
+      gEnv->SetValue("XSec.GSI.DelegProxy", "2");
+      
+      // --- Now open connection to PROOF cluster --------------------
+      TString serv = "";
+      Bool_t  lite = false;
+      if (fProofServer.BeginsWith("workers=") || fProofServer.IsNull()) {
+       lite = true;
+       serv = fProofServer;
+      }
+      else 
+       serv = Form("%s@%s", userName.Data(), fProofServer.Data());
+      TProof::Open(serv);
+      if (!gProof) { 
+       Error("Connect", "Failed to connect to Proof cluster %s as %s",
+             fProofServer.Data(), userName.Data());
+       return false;
+      }
+      if (lite) return true;
+    }
+
+    // --- Open a connection to the grid -----------------------------
+    TGrid::Connect("alien://");
+    if (!gGrid || !gGrid->IsConnected()) { 
+      // This is only fatal in grid mode 
+      Error("Connect", "Failed to connect to AliEN");
+      if (mode == kGrid) return false; 
+      return true;
+    }
+    if (mode == kGrid) return true;
+
+    
+    // --- Set and make output directory -----------------------------
+    TString name = EscapedName();
+    TString homeDir(gGrid->GetHomeDirectory());
+    TString workDir(homeDir);
+    workDir.Append("/");
+    workDir.Append(name);
+    
+    // Make working directory 
+    if (!gGrid->Cd(workDir)) { 
+      gGrid->Cd(homeDir);
+      if (gGrid->Mkdir(workDir)) {
+       gGrid->Cd(name);
+       Info("Connect", "Directory %s created", workDir.Data());
+      }
+    }
+    // Make output directory 
+    gGrid->Mkdir("proof_output");
+    gGrid->Cd("proof_output");
+
+    return true;
+  }      
+  //__________________________________________________________________
+  /** 
+   * Load common libraries 
+   * 
+   * @param mode Running mode                  
+   * @param par  If true, load as PARs 
+   * 
+   * @return true on success 
+   */
+  Bool_t LoadCommonLibraries(EMode mode, Bool_t par) 
+  {
+    if (!gSystem->Getenv("ALICE_ROOT")) { 
+      Error("LoadCommonLibraries", "Local AliROOT not available");
+      return false;
+    }
+    gSystem->Load("libTree.so");
+    gSystem->Load("libGeom.so");
+    gSystem->Load("libVMC.so");
+    gSystem->Load("libPhysics.so");
+    gSystem->Load("libMinuit.so");
+
+    Bool_t ret   = true;
+    Bool_t basic = mode == kGrid ? false : par;
+    
+    ret &= LoadLibrary("STEERBase",     mode, basic, false);
+    ret &= LoadLibrary("ESD",           mode, basic, false);
+    ret &= LoadLibrary("AOD",           mode, basic, false);
+    ret &= LoadLibrary("ANALYSIS",      mode, basic, true);
+    ret &= LoadLibrary("ANALYSISalice", mode, basic, true);
+
+    return ret;
+  }
+  //__________________________________________________________________
+  /** 
+   * Load a library 
+   * 
+   * @param what What library to load
+   * @param mode Mode (local, proof, grid)
+   * @param par  If true, load as PAR
+   * @param rec  If true, also load on slaves
+   * 
+   * @return true on success 
+   */
+  Bool_t LoadLibrary(const char* what, EMode mode, Bool_t par, Bool_t rec=false)
+  {
+    if (!what || what[0] == '\0') return true;
+    
+    TString module(what);
+    TString libName(what);
+    if (!libName.BeginsWith("lib")) libName = Form("lib%s", libName.Data());
+    if (!libName.EndsWith(".so"))   libName.Append(".so");
+
+    Int_t ret = 0;
+
+    switch (mode) { 
+    case kLocal: // Just load and exit 
+      gSystem->Load(libName.Data());
+      break;
+    case kGrid: 
+      if (par) { 
+       ret = SetupPAR(what) ? 0 : -1;
+       if (rec) fListOfPARs.Add(new TObjString(what));
+      } else  {
+       ret = gSystem->Load(libName.Data());
+       if (rec) fListOfLibraries.Add(new TObjString(libName));
+      }
+      break;
+    case kProof: 
+      ret = gProof->UploadPackage(what);
+      if (ret < 0)  {  
+         ret = gProof->UploadPackage(gSystem->ExpandPathName(Form("../%s.par",
+                                                                  what)));
+       if (ret < 0) {  
+         ret = 
+           gProof->UploadPackage(gSystem
+                                 ->ExpandPathName(Form("$ALICE_ROOT/%s.par", 
+                                                       what)));
+         if (ret < 0) {
+           Error("LoadLibrary", 
+                 "Could not find module %s.par in current directory nor "
+                 "in $ALICE_ROOT", module.Data());
+           return false;
+         }
+       }
+      }
+      ret = gProof->EnablePackage(what);
+      break;
+    }
+    if (ret < 0) { 
+      Error("LoadLibrary", "Couldn't load %s", what);
+      return false;
+    }
+    return true;
+  }
+          
+  //__________________________________________________________________
+  Bool_t SetupPAR(const char* what)
+  {
+    if (!what || what[0] == '\0') return -1;
+    
+    TString parFile(Form("%s.par", what));
+    if (gSystem->AccessPathName(parFile.Data())) { 
+      if (gSystem->AccessPathName(Form("../%s.par", what))) { 
+       // If not found 
+       TString aliParFile = 
+         gSystem->ExpandPathName(Form("$(ALICE_ROOT)/%s.par", what));
+       if (gSystem->AccessPathName(aliParFile.Data())) { 
+         Error("SetupPAR", "PAR file %s not found in current directory or "
+               "$(ALICE_ROOT)", what);
+         return false;
+       }
+       // Copy to current directory 
+       TFile::Cp(aliParFile, parFile);
+      }
+      else 
+       gSystem->Exec(Form("ln -s ../%s.par .", what));
+    }
+    
+    // Extract archive 
+    gSystem->Exec(Form("tar xvzf %s", parFile.Data()));
+    
+    // Change directory into par archive
+    TString cwd = gSystem->WorkingDirectory();
+    
+    if (!gSystem->ChangeDirectory(what)) { 
+      Error("SetupPAR", "Failed to change directory to %s", what);
+      return false;
+    }
+    
+    // Test the build 
+    if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
+      Info("SetupPar", "Building in PAR archive %s", what);
+      if (gSystem->Exec("PROOF-INF/BUILD.sh")) { 
+       Error("SetupPar", "Failed to build in PAR directory %s", what);
+       gSystem->ChangeDirectory(cwd.Data());
+       return false;
+      }
+    }
+    
+    // Check for setup script
+    if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
+      Info("SetupPAR", "Setting up for PAR %s", what);
+      gROOT->Macro("PROOF-INF/SETUP.C");
+    }
+    if (!gSystem->ChangeDirectory(cwd.Data())) return false;
+
+    return true;
+  }
+  //__________________________________________________________________
+  TString SetupExtras()
+  {
+    TString ret;
+    TIter next(&fListOfExtras);
+    TObjString* obj = 0;
+    while ((obj = static_cast<TObjString*>(next()))) {
+      TString path = gSystem->ExpandPathName(obj->GetName());
+      if (!path.BeginsWith("/")) 
+       // If not an absolute path, prepend to up-one
+       path = Form("../%s", path.Data());
+      if (gSystem->AccessPathName(path.Data())) { 
+       // File not accessible 
+       Warning("SetupExtras", "File %s not accessible", path.Data());
+       continue;
+      }
+      ret.Append(Form("%s ", gSystem->BaseName(obj->GetName())));
+      gSystem->Exec(Form("ln -s %s .", path.Data()));
+    }
+    ret = ret.Strip();
+    return ret;
+  }
+  //__________________________________________________________________
+  TString SetupSources()
+  {
+    TString nam = EscapedName();
+    TString ret;
+    TIter next(&fListOfSources); 
+    TObject* src;
+    while ((src = next())) {
+      TString path = gSystem->ExpandPathName(src->GetName());
+      if (!path.BeginsWith("/")) 
+       // If not an absolute path, prepend to up-one
+       path = Form("../%s", path.Data());
+      if (gSystem->AccessPathName(path.Data())) { 
+       // File not accessible 
+       Warning("SetupSources", "File %s not accessible", path.Data());
+       continue;
+      }
+      ret.Append(Form("%s ", gSystem->BaseName(src->GetName())));
+      gSystem->Exec(Form("ln -s %s .", path.Data()));
+    }
+    ret = ret.Strip();
+    return ret;
+  }
+  //__________________________________________________________________
+  TString SetupLibraries()
+  {
+    TString ret;
+    TIter next(&fListOfLibraries); 
+    TObject* lib;
+    while ((lib = next())) {
+      ret.Append(lib->GetName());
+      ret.Append(" ");
+    }
+    // Also add extra files to this variable 
+    ret.Append(SetupExtras());
+    ret = ret.Strip();
+    return ret;
+  }
+  //__________________________________________________________________
+  /** 
+   * Scan directory @a dir (possibly recursive) for tree files to add
+   * to the chain.    This does not follow sym-links
+   * 
+   * @param dir        Directory to scan
+   * @param chain      Chain to add to
+   * @param type       Type of tree (ESD or AOD)
+   * @param recursive  Whether to scan recursively 
+   * @param mc         Look also for MC files if true 
+   *
+   * @return true if any files where added 
+   */
+  Bool_t ScanDirectory(TSystemDirectory* dir, TChain* chain, 
+                      EType type, bool recursive, bool mc)
+  {
+    TString fnPattern;
+    switch (type) { 
+    case kESD:  fnPattern = "AliESD"; break;
+    case kAOD:  fnPattern = "AliAOD"; break;
+    }
+
+    // Assume failure 
+    Bool_t ret = false;
+
+    // Get list of files, and go back to old working directory
+    TString oldDir(gSystem->WorkingDirectory());
+    TList*  files = dir->GetListOfFiles();
+    if (!gSystem->ChangeDirectory(oldDir)) { 
+      Error("ScanDirectory", "Failed to go back to %s", oldDir.Data());
+      return false;
+    }
+    if (!files) return false;
+
+    TList toAdd;
+    toAdd.SetOwner();
+    Bool_t hasGAlice = (!mc ? true : false);
+    Bool_t hasKine   = (!mc ? true : false);
+    Bool_t hasTrRef  = (!mc ? true : false);
+    
+    // Sort list of files and check if we should add it 
+    files->Sort();
+    TIter next(files);
+    TSystemFile* file = 0;
+    while ((file = static_cast<TSystemFile*>(next()))) {
+      TString name(file->GetName());
+      TString title(file->GetTitle());
+      TString full(gSystem->ConcatFileName(file->GetTitle(), name.Data()));
+      if (dynamic_cast<TSystemDirectory*>(file)) full = title;
+      // Ignore special links 
+      if (name == "." || name == "..") { 
+       // Info("ScanDirectory", "Ignoring %s", name.Data());
+       continue;
+      }
+
+      FileStat_t fs;
+      if (gSystem->GetPathInfo(full.Data(), fs)) {
+       Warning("ScanDirectory", "Cannot stat %s (%s)", full.Data(),
+                gSystem->WorkingDirectory());
+       continue;
+      }
+      // Check if this is a directory 
+      if (file->IsDirectory(full)) { 
+       if (recursive) {
+         // if (title[0] == '/') 
+         TSystemDirectory* d = new TSystemDirectory(file->GetName(),
+                                                    full.Data());
+          if (ScanDirectory(d,chain,type,recursive,mc))
+           ret = true;
+         delete d;
+       }
+        continue;
+      }
+    
+      // If this is not a root file, ignore 
+      if (!name.EndsWith(".root")) continue;
+
+      // If this file does not contain AliESDs, ignore 
+      if (!name.Contains(fnPattern)) { 
+       // Info("ScanDirectory", "%s does not match pattern %s", 
+       //      name.Data(), fnPattern.Data());
+       if (mc) { 
+         if (name.CompareTo("galice.root") == 0)     hasGAlice = true;
+         if (name.CompareTo("Kinematics.root") == 0) hasKine   = true;
+         if (name.CompareTo("TrackRefs.root")  == 0) hasTrRef = true;
+       }
+       continue;
+      }
+    
+      // Add 
+      // Info("ScanDirectory", "Adding %s", full.Data());
+      toAdd.Add(new TObjString(full));
+    }
+
+    if (mc && toAdd.GetEntries() > 0 && 
+       (!hasGAlice || !hasKine || !hasTrRef)) { 
+      Warning("ScanDirectory", 
+             "one or more of {galice,Kinematics,TrackRefs}.root missing from "
+             "%s, not adding anything from this directory", 
+             dir->GetTitle());
+      toAdd.Delete();
+    }
+
+    TIter nextAdd(&toAdd);
+    TObjString* s = 0;
+    while ((s = static_cast<TObjString*>(nextAdd()))) {
+      // Info("ScanDirectory", "Adding %s", s->GetString().Data());
+      chain->Add(s->GetString());
+    }
+    if (toAdd.GetEntries() > 0) ret = true;
+
+    gSystem->ChangeDirectory(oldDir);
+    return ret;
+  }
+  //__________________________________________________________________
+  /** 
+   * Create a chain from an XML containing an collection
+   * 
+   * @param treeName Name of tree's 
+   * @param xmlFile  XML collection
+   * 
+   * @return Newly allocated chain or null
+   */
+  TChain* CreateChainFromXML(const char* treeName, 
+                            const char* xmlFile) 
+  {
+    TGridCollection* collection = TAlienCollection::Open(xmlFile);
+    if (!collection) { 
+      Error("CreateChainFromXML", "Cannot create AliEn collection from "
+           "XML file %s", xmlFile);
+      return 0;
+    }
+
+    TChain* chain = new TChain(treeName);
+    collection->Reset();
+    while (collection->Next()) chain->Add(collection->GetTURL(""));
+    
+    return chain;
+  }
+  //__________________________________________________________________
+  /** 
+   * Create a chain of data 
+   * 
+   * @param type Type of data
+   * @param mode Operation mode 
+   * @param mc   Assume MC input if true
+   *
+   * @return TChain of data 
+   */    
+  TChain* CreateChain(EType type, EMode mode, EOper /* oper */, Bool_t mc)
+  {
+    TString treeName;
+    switch (type) { 
+    case kESD:  treeName = "esdTree"; break;
+    case kAOD:  treeName = "aodTree"; break;
+    }
+
+    TChain* chain = 0;
+    switch (mode) { 
+    case kProof: 
+      if (!fDataSet.IsNull()) break; 
+      // Otherwise fall through
+    case kLocal:
+      if (fXML.IsNull()) {
+       chain = new TChain(treeName.Data());
+       TString dir(fDataDir);
+       if (dir == ".") dir = "";
+       if (!dir.BeginsWith("/")) dir = Form("../%s", dir.Data());
+       TString savdir(gSystem->WorkingDirectory());
+       TSystemDirectory d(gSystem->BaseName(dir.Data()), dir.Data());
+       if (!ScanDirectory(&d, chain, type, true, mc)) { 
+         delete chain;
+         chain = 0;
+       }
+       gSystem->ChangeDirectory(savdir);
+      }
+      else 
+       chain = CreateChainFromXML(treeName.Data(), fXML.Data());
+      break;
+    case kGrid:  break; // Do nothing - we use plugin
+    }
+    
+    if (chain && chain->GetNtrees() <= 0) { 
+      delete chain;
+      return 0;
+    }
+    return chain;
+  }
+  //__________________________________________________________________
+  TString fName;             // Name of analysis
+  TString fRootVersion;      // ROOT version to use 
+  TString fAliRootVersion;   // AliROOT version to use 
+  TString fAliEnAPIVersion;  // AliEn API version to use 
+  TString fProofServer;      // Name of proof server
+  TString fDataDir;          // Grid Input directory 
+  TString fDataSet;          // Proof data set name 
+  TString fXML;              // XML collection for local/proof mode
+  TArrayI fRunNumbers;       // List of run number 
+  TList   fListOfPARs;       // List of PAR files to use 
+  TList   fListOfSources;    // List of sources to upload and AcLIC
+  TList   fListOfLibraries;  // List of libraries to load
+  TList   fListOfExtras;     // List of extra files to upload
+  Int_t   fNReplica;         // Storage replication
+  Int_t   fESDPass;
+  TString fPassPostfix;      // Possible pass postfix
+  TString fEscapedName;
+  Bool_t  fAllowOverwrite;
+};
+
+
+void
+BuildTrainSetup()
+{
+  gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
+                          "$ALICE_ROOT/ANALYSIS/macros",
+                          gROOT->GetMacroPath()));
+  gSystem->AddIncludePath("-I${ALICE_ROOT}/include");
+  gSystem->Load("libRAliEn");
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  TString path = gSystem->Which(gROOT->GetMacroPath(), "TrainSetup.C");
+  Info("BuildTrainSetup", "Path=%s", path.Data());
+  TString tmp("TrainSetup");
+  FILE* fp = gSystem->TempFileName(tmp, ".");
+  fclose(fp);
+  gSystem->Unlink(tmp);
+  tmp.Append(".C");
+  Info("BuildTrainSetup", "Copy %s -> %s", path.Data(), tmp.Data());
+  gSystem->CopyFile(path, tmp);
+  gROOT->LoadMacro(Form("%s+g", tmp.Data()));
+  gSystem->Unlink(tmp);
+  tmp.ReplaceAll(".C", "_C.so");
+  gSystem->Unlink(tmp);
+  tmp.ReplaceAll("_C.so", "_C.d");
+  gSystem->Unlink(tmp);
+}
+
+  
+//____________________________________________________________________
+//
+// EOF
+//