Refactoring for AliFMDEnergyFitter to allow sub-classing
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Nov 2013 13:32:36 +0000 (13:32 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Nov 2013 13:32:36 +0000 (13:32 +0000)
for use with MC hit information.  Also decode ELoss from
AliTrackReference, etc.  Some code clean-up to remove
virtual'ness of members that shouldn't be virtual, etc.

Some fixes so that the cache histograms of
AliForwardMultiplicityBase are initialized always.

25 files changed:
PWGLF/FORWARD/analysis2/AddTaskFMDELoss.C
PWGLF/FORWARD/analysis2/AddTaskFMDMCHit.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliBaseESDTask.cxx
PWGLF/FORWARD/analysis2/AliBaseESDTask.h
PWGLF/FORWARD/analysis2/AliFMDEnergyFitter.cxx
PWGLF/FORWARD/analysis2/AliFMDEnergyFitter.h
PWGLF/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx
PWGLF/FORWARD/analysis2/AliFMDEnergyFitterTask.h
PWGLF/FORWARD/analysis2/AliFMDEventInspector.cxx
PWGLF/FORWARD/analysis2/AliFMDEventInspector.h
PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitter.cxx [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitter.h [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitterTask.cxx [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitterTask.h [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.h
PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx
PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h
PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.h
PWGLF/FORWARD/analysis2/AliForwardQATask.cxx
PWGLF/FORWARD/analysis2/AliForwardQATask.h
PWGLF/FORWARD/analysis2/AliMCAuxHandler.cxx
PWGLF/FORWARD/analysis2/ForwardAODConfig.C
PWGLF/FORWARD/analysis2/scripts/LoadLibs.C
PWGLF/FORWARD/analysis2/trains/MakeFMDELossTrain.C
PWGLF/FORWARD/analysis2/trains/MakeFMDMCHitTrain.C [new file with mode: 0644]

index f957dfe..c49763d 100644 (file)
@@ -95,7 +95,6 @@ AddTaskFMDELoss(Bool_t        mc,
     else // Anything else gives plain difference and errors in errors
       rm = AliFMDEnergyFitter::kResidualDifference;
   }
-  Printf("Got residual: \"%s\" -> %d", resi.Data(), rm);
   task->GetEnergyFitter().SetStoreResiduals(rm);
 
   // --- Set limits on fits the energy -------------------------------
@@ -110,5 +109,9 @@ AddTaskFMDELoss(Bool_t        mc,
   // --- Make the output container and connect it --------------------
   task->Connect(0,0);
 
+  Printf("Returning task %p", task);
   return task;
 }
+//
+// EOF
+//
diff --git a/PWGLF/FORWARD/analysis2/AddTaskFMDMCHit.C b/PWGLF/FORWARD/analysis2/AddTaskFMDMCHit.C
new file mode 100644 (file)
index 0000000..d2628ab
--- /dev/null
@@ -0,0 +1,111 @@
+/**
+ * @file   AddTaskFMDELoss.C
+ * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
+ * @date   Wed Mar 23 12:14:03 2011
+ * 
+ * @brief  
+ * 
+ * 
+ * @ingroup pwglf_forward_scripts_tasks
+ */
+/**
+ * @defgroup pwglf_forward_eloss Energy Loss Fits
+ *
+ * Fitting the energy loss @f$\Delta/\Delta_{mip}@f$ spectra 
+ *
+ * @ingroup pwglf_forward_topical
+ */
+
+/**
+ * This is the macro to include the FMD energy fitter in a train.  
+ * 
+ * @param useTuple  If true, create NTuple of hits
+ * @param debug     Debug level
+ *
+ * @return Newly created task 
+ *
+ * @ingroup pwglf_forward_eloss
+ */
+AliAnalysisTask*
+AddTaskFMDMCHit(Bool_t useTuple=false, 
+               Int_t  debug=0)
+{
+  // --- Load libraries ----------------------------------------------
+  gROOT->LoadClass("AliFMDDigit",             "FMDbase");
+  gROOT->LoadClass("AliFMDHit",               "FMDsim");
+  gROOT->LoadClass("AliAODForwardMult",       "PWGLFforward2");
+  gROOT->LoadClass("AliFMDMCHitEnergyFitter", "PWGLFforwardhit");
+
+  // --- Get analysis manager ----------------------------------------
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("AddTaskFMDELoss", "No analysis manager to connect to.");
+    return NULL;
+  }   
+
+  // --- Make the task and add it to the manager ---------------------
+  AliFMDMCHitEnergyFitterTask* task = 
+    new AliFMDMCHitEnergyFitterTask("ForwardHitELoss", 
+                                   useTuple);
+  // --- Set parameters on the algorithms ----------------------------
+  // Set the number of SPD tracklets for which we consider the event a
+  // low flux event
+  task->GetEventInspector().SetLowFluxCut(1000); 
+  // Set the maximum error on v_z [cm]
+  task->GetEventInspector().SetMaxVzErr(0.2);
+  // Set the eta axis to use - note, this overrides whatever is used
+  // by the rest of the algorithms - but only for the energy fitter
+  // algorithm. 
+  task->GetEnergyFitter().SetEtaAxis(100, -4, 6);
+  // Set maximum energy loss to consider 
+  task->GetEnergyFitter().SetMaxE(15); 
+  // Set number of energy loss bins 
+  task->GetEnergyFitter().SetNEbins(500);
+  // Set whether to use increasing bin sizes 
+  task->GetEnergyFitter().SetUseIncreasingBins(true);
+  // Set whether to do fit the energy distributions 
+  task->GetEnergyFitter().SetDoFits(kTRUE);
+  // Set whether to make the correction object 
+  // task->GetEnergyFitter().SetDoMakeObject(kTRUE);
+  // Set the low cut used for energy
+  task->GetEnergyFitter().SetLowCut(0.05);
+  // Set the number of bins to subtract from maximum of distributions
+  // to get the lower bound of the fit range
+  task->GetEnergyFitter().SetFitRangeBinWidth(4);
+  // Set the maximum number of landaus to try to fit (max 5)
+  task->GetEnergyFitter().SetNParticles(1);
+  // Set the minimum number of entries in the distribution before
+  // trying to fit to the data - 10K seems the absolute minimum
+  task->GetEnergyFitter().SetMinEntries(3000 /*10000*/);
+  // If set, only collect statistics for MB.  This is to prevent a
+  // bias when looping over data where the MB trigger is downscaled.
+  // task->SetOnlyMB(onlyMB);
+  // Debug
+  task->SetDebug(debug);
+
+  // --- Set limits on fits the energy -------------------------------
+  // DO NOT CHANGE THESE UNLESS YOU KNOW WHAT YOU ARE DOING
+  // Maximum relative error on parameters 
+  // AliFMDCorrELossFit::ELossFit::fgMaxRelError = .12;
+  // Least weight to use 
+  // AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
+  // Maximum value of reduced chi^2 
+  // AliFMDCorrELossFit::ELossFit::fgMaxChi2nu   = 20;
+  
+  // --- Make the output container and connect it --------------------
+  task->Connect(0,0);
+  if (useTuple) { 
+    AliAnalysisDataContainer* tuple = 
+      mgr->CreateContainer("tuple", TNtuple::Class(), 
+                          AliAnalysisManager::kOutputContainer,
+                          "forward_tuple.root"
+                          /*AliAnalysisManager::GetCommonFileName()*/);
+    mgr->ConnectOutput(task, 3, tuple);
+  }
+
+  Printf("Returning task %p", task);
+  return task;
+}
+//
+// EOF
+//
index a36aaa6..edd88e6 100644 (file)
@@ -67,18 +67,30 @@ AliBaseESDTask::Connect(const char* sumFile,
   if      (sumFile && sumFile[0] != '\0') sumOut = sumFile;
   if      (resFile && resFile[0] != '\0') resOut = resFile;
   else if (sumFile && sumFile[0] != '\0') resOut = sumFile;
-  if (sumOut.IsNull()) sumOut = AliAnalysisManager::GetCommonFileName();
-  if (resOut.IsNull()) resOut = AliAnalysisManager::GetCommonFileName();
-
-  AliAnalysisDataContainer* sumCon = 
-    mgr->CreateContainer(Form("%sSums", GetName()), TList::Class(), 
-                        AliAnalysisManager::kOutputContainer, sumOut);
-  AliAnalysisDataContainer* resCon = 
-    mgr->CreateContainer(Form("%sResults", GetName()), TList::Class(), 
-                        AliAnalysisManager::kParamContainer, resOut);
+  // If the string is null or 'default' connect to standard output file 
+  if (sumOut.IsNull() || sumOut.EqualTo("default", TString::kIgnoreCase)) 
+    sumOut = AliAnalysisManager::GetCommonFileName();
+  // If the string is null or 'default' connect to standard output file 
+  if (resOut.IsNull() || resOut.EqualTo("default", TString::kIgnoreCase)) 
+    resOut = AliAnalysisManager::GetCommonFileName();
+
+  // Always connect input 
   mgr->ConnectInput(this, 0, mgr->GetCommonInputContainer());
-  mgr->ConnectOutput(this, 1, sumCon);
-  mgr->ConnectOutput(this, 2, resCon);
+
+  // Connect sum list unless the output 'none' is specified
+  if (!sumOut.EqualTo("none", TString::kIgnoreCase)) {
+    AliAnalysisDataContainer* sumCon = 
+      mgr->CreateContainer(Form("%sSums", GetName()), TList::Class(), 
+                          AliAnalysisManager::kOutputContainer, sumOut);
+    mgr->ConnectOutput(this, 1, sumCon);
+  }
+  // Connect the result list unless the output 'none' is specified
+  if (!resOut.EqualTo("none", TString::kIgnoreCase)) {
+    AliAnalysisDataContainer* resCon = 
+      mgr->CreateContainer(Form("%sResults", GetName()), TList::Class(), 
+                          AliAnalysisManager::kParamContainer, resOut);
+    mgr->ConnectOutput(this, 2, resCon);
+  }
   
   return true;
 }
@@ -122,6 +134,15 @@ AliBaseESDTask::Configure(const char* macro)
  
  return true;
 }
+
+//____________________________________________________________________
+void 
+AliBaseESDTask::LocalInit() 
+{ 
+  fFirstEvent = true; 
+  Setup(); 
+}
+
 //____________________________________________________________________
 void
 AliBaseESDTask::UserCreateOutputObjects()
@@ -144,7 +165,8 @@ AliBaseESDTask::UserCreateOutputObjects()
   GetEventInspector().CreateOutputObjects(fList);
 
   if (!Book()) AliFatalF("Failed to book output objects for %s", GetName());
-  
+
+  // gSystem->Exec("root-config --version --prefix");
   PostData(1, fList);
 }
 
index de2eb80..16a0070 100644 (file)
@@ -105,6 +105,10 @@ public:
    */
   virtual Bool_t Connect(const char* sumFile=0, const char* resFile=0);
   /** 
+   * Called when initializing the train 
+   */
+  virtual Bool_t Setup() { return true; }
+  /** 
    * Book output objects. Derived class should define this to book
    * output objects on the processing output list @c fList before the
    * actual event processing.  This is called on the master and on
@@ -270,23 +274,23 @@ protected:
   /** 
    * Initialize the task 
    */
-  virtual void Init() { fFirstEvent = true; }
+  void LocalInit();
   /** 
    * Create output objects 
    */
-  virtual void UserCreateOutputObjects();
+  void UserCreateOutputObjects();
   /** 
    * Process each event 
    *
    * @param option Not used
    */  
-  virtual void UserExec(Option_t* option);
+  void UserExec(Option_t* option);
   /** 
    * End of job
    * 
    * @param option Not used 
    */
-  virtual void Terminate(Option_t* option);
+  void Terminate(Option_t* option);
   /** 
    * @} 
    */
index affcc81..e504866 100644 (file)
@@ -7,6 +7,7 @@
 #include <TAxis.h>
 #include <TList.h>
 #include <TH1.h>
+#include <TH2.h>
 #include <TF1.h>
 #include <TMath.h>
 #include <AliLog.h>
@@ -90,46 +91,6 @@ AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
   fEtaAxis.SetTitle("#eta");
   fCentralityAxis.SetName("centralityAxis");
   fCentralityAxis.SetTitle("Centrality [%]");
-  fRingHistos.Add(new RingHistos(1, 'I'));
-  fRingHistos.Add(new RingHistos(2, 'I'));
-  fRingHistos.Add(new RingHistos(2, 'O'));
-  fRingHistos.Add(new RingHistos(3, 'I'));
-  fRingHistos.Add(new RingHistos(3, 'O'));
-}
-
-//____________________________________________________________________
-AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
-  : TNamed(o), 
-    fRingHistos(), 
-    fLowCut(o.fLowCut),
-    fNParticles(o.fNParticles),
-    fMinEntries(o.fMinEntries),
-    fFitRangeBinWidth(o.fFitRangeBinWidth),
-    fDoFits(o.fDoFits),
-    fDoMakeObject(o.fDoMakeObject),
-    fEtaAxis(o.fEtaAxis),
-    fCentralityAxis(o.fCentralityAxis),
-    fMaxE(o.fMaxE),
-    fNEbins(o.fNEbins), 
-    fUseIncreasingBins(o.fUseIncreasingBins),
-    fMaxRelParError(o.fMaxRelParError),
-    fMaxChi2PerNDF(o.fMaxChi2PerNDF), 
-    fMinWeight(o.fMinWeight),
-    fDebug(o.fDebug),
-    fResidualMethod(o.fResidualMethod),
-    fSkips(o.fSkips),
-    fRegularizationCut(o.fRegularizationCut)
-{
-  // 
-  // Copy constructor 
-  // 
-  // Parameters:
-  //    o Object to copy from 
-  //
-  DGUARD(fDebug, 1, "Copy CTOR of AliFMDEnergyFitter");
-  TIter    next(&o.fRingHistos);
-  TObject* obj = 0;
-  while ((obj = next())) fRingHistos.Add(obj);
 }
 
 //____________________________________________________________________
@@ -143,54 +104,10 @@ AliFMDEnergyFitter::~AliFMDEnergyFitter()
 }
 
 //____________________________________________________________________
-AliFMDEnergyFitter&
-AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
+AliFMDEnergyFitter::RingHistos*
+AliFMDEnergyFitter::CreateRingHistos(UShort_t d, Char_t r) const
 {
-  // 
-  // Assignment operator 
-  // 
-  // Parameters:
-  //    o Object to assign from 
-  // 
-  // Return:
-  //    Reference to this 
-  //
-  fDebug = o.fDebug;
-  DGUARD(fDebug, 3, "Assignment of AliFMDEnergyFitter");
-  if (&o == this) return *this; 
-  TNamed::operator=(o);
-
-  fLowCut        = o.fLowCut;
-  fNParticles       = o.fNParticles;
-  fMinEntries    = o.fMinEntries;
-  fFitRangeBinWidth= o.fFitRangeBinWidth;
-  fDoFits        = o.fDoFits;
-  fDoMakeObject  = o.fDoMakeObject;
-  fEtaAxis.Set(o.fEtaAxis.GetNbins(),
-              o.fEtaAxis.GetXmin(),
-              o.fEtaAxis.GetXmax());
-  if (o.fCentralityAxis.GetXbins()) {
-    const TArrayD* bins = o.fCentralityAxis.GetXbins();
-    fCentralityAxis.Set(bins->GetSize()-1,bins->GetArray());
-  }
-  else 
-    fCentralityAxis.Set(o.fCentralityAxis.GetNbins(),
-                       o.fCentralityAxis.GetXmin(),
-                       o.fCentralityAxis.GetXmax());
-  fDebug             = o.fDebug;
-  fMaxRelParError    = o.fMaxRelParError;
-  fMaxChi2PerNDF     = o.fMaxChi2PerNDF;
-  fMinWeight         = o.fMinWeight;
-  fResidualMethod    = o.fResidualMethod;
-  fSkips             = o.fSkips;
-  fRegularizationCut = o.fRegularizationCut;
-
-  fRingHistos.Clear();
-  TIter    next(&o.fRingHistos);
-  TObject* obj = 0;
-  while ((obj = next())) fRingHistos.Add(obj);
-  
-  return *this;
+  return new RingHistos(d,r);
 }
 
 //____________________________________________________________________
@@ -220,10 +137,87 @@ AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
 
 //____________________________________________________________________
 void
+AliFMDEnergyFitter::Init()
+{
+  // Create the ring histograms.  
+  // 
+  // Should be called from task initialization. 
+  DGUARD(1,fDebug, "Creating histogram caches for each ring");
+  fRingHistos.Add(CreateRingHistos(1, 'I'));
+  fRingHistos.Add(CreateRingHistos(2, 'I'));
+  fRingHistos.Add(CreateRingHistos(2, 'O'));
+  fRingHistos.Add(CreateRingHistos(3, 'I'));
+  fRingHistos.Add(CreateRingHistos(3, 'O'));
+  TIter    next(&fRingHistos);
+  RingHistos* o = 0;
+  while ((o = static_cast<RingHistos*>(next()))) {
+    o->fDebug = fDebug;
+  }
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::CreateOutputObjects(TList* dir)
+{
+  // 
+  // Define the output histograms.  These are put in a sub list of the
+  // passed list.   The histograms are merged before the parent task calls 
+  // AliAnalysisTaskSE::Terminate.  Called on task initialization on slaves.  
+  // 
+  // Parameters:
+  //    dir Directory to add to 
+  //
+  DGUARD(fDebug, 1, "Define output in AliFMDEnergyFitter");
+  TList* d = new TList;
+  d->SetName(GetName());
+  d->SetOwner(true);
+  dir->Add(d);
+
+  // Store eta axis as a histogram, since that can be merged!
+  TH1* hEta = 0;
+  if (fEtaAxis.GetXbins()->GetArray()) 
+    hEta = new TH1I(fEtaAxis.GetName(), fEtaAxis.GetTitle(), 
+                   fEtaAxis.GetNbins(), fEtaAxis.GetXbins()->GetArray());
+  else 
+    hEta = new TH1I(fEtaAxis.GetName(), fEtaAxis.GetTitle(), 
+                   fEtaAxis.GetNbins(),fEtaAxis.GetXmin(),fEtaAxis.GetXmax());
+  hEta->SetXTitle("#eta");
+  hEta->SetYTitle("Nothing");
+  hEta->SetDirectory(0);
+
+  d->Add(hEta);
+  d->Add(AliForwardUtil::MakeParameter("lowCut",        fLowCut));
+  d->Add(AliForwardUtil::MakeParameter("nParticles",    fNParticles));
+  d->Add(AliForwardUtil::MakeParameter("minEntries",    fMinEntries));
+  d->Add(AliForwardUtil::MakeParameter("subtractBins",  fFitRangeBinWidth));
+  d->Add(AliForwardUtil::MakeParameter("doFits",        fDoFits));
+  d->Add(AliForwardUtil::MakeParameter("doObject",      fDoMakeObject));
+  d->Add(AliForwardUtil::MakeParameter("maxE",          fMaxE));
+  d->Add(AliForwardUtil::MakeParameter("nEbins",        fNEbins));
+  d->Add(AliForwardUtil::MakeParameter("increasingBins",fUseIncreasingBins));
+  d->Add(AliForwardUtil::MakeParameter("maxRelPerError",fMaxRelParError));
+  d->Add(AliForwardUtil::MakeParameter("maxChi2PerNDF", fMaxChi2PerNDF));
+  d->Add(AliForwardUtil::MakeParameter("minWeight",     fMinWeight));
+  d->Add(AliForwardUtil::MakeParameter("regCut",        fRegularizationCut));
+
+  if (fRingHistos.GetEntries() <= 0) { 
+    AliFatal("No ring histograms where defined - giving up!");
+    return;
+  }
+  TIter    next(&fRingHistos);
+  RingHistos* o = 0;
+  while ((o = static_cast<RingHistos*>(next()))) {
+    o->fDebug = fDebug;
+    o->CreateOutputObjects(d);
+  }
+}
+
+//____________________________________________________________________
+void
 AliFMDEnergyFitter::SetupForData(const TAxis& eAxis)
 {
   // 
-  // Initialise the task
+  // Initialise the task - called at first event 
   // 
   // Parameters:
   //    etaAxis The eta axis to use.  Note, that if the eta axis
@@ -304,7 +298,7 @@ AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
   // 
   // Return:
   //    True on success, false otherwise 
-  DGUARD(fDebug, 3, "Accumulate statistics in AliFMDEnergyFitter - cholm");
+  DGUARD(fDebug, 5, "Accumulate statistics in AliFMDEnergyFitter - cholm");
   Int_t icent = fCentralityAxis.FindBin(cent);
   if (icent < 1 || icent > fCentralityAxis.GetNbins()) icent = 0;
 
@@ -317,6 +311,7 @@ AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
       UShort_t    nstr = (q == 0 ? 512 : 256);
 
       RingHistos* histos = GetRingHistos(d, r);
+      if (!histos) continue;
       
       for(UShort_t s = 0; s < nsec;  s++) {
        for(UShort_t t = 0; t < nstr; t++) {
@@ -328,10 +323,11 @@ AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
 
          // Get the pseudo-rapidity 
          Double_t eta1 = input.Eta(d,r,s,t);
-         Int_t ieta = fEtaAxis.FindBin(eta1);
-         if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue; 
+         // Int_t ieta = fEtaAxis.FindBin(eta1);
+         // if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue; 
 
-         histos->Fill(empty, ieta-1, icent, mult);
+         // histos->Fill(empty, ieta-1, icent, mult);
+         histos->Fill(empty, eta1, icent, mult);
          nFills++;
        } // for strip
       } // for sector
@@ -388,7 +384,7 @@ AliFMDEnergyFitter::Fit(const TList* dir)
       continue;
     }
     
-    TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNParticles,
+    TObjArray* l = o->Fit(d, fLowCut, fNParticles,
                          fMinEntries, fFitRangeBinWidth,
                          fMaxRelParError, fMaxChi2PerNDF,
                          fMinWeight, fRegularizationCut,
@@ -428,53 +424,13 @@ AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
       continue;
     }
     
-    o->FindBestFits(d, *obj, fEtaAxis, fMaxRelParError, 
-                   fMaxChi2PerNDF, fMinWeight);
+    o->FindBestFits(d, *obj, fEtaAxis);
   }
   d->Add(obj, "elossFits");
 }
 
 //____________________________________________________________________
 void
-AliFMDEnergyFitter::CreateOutputObjects(TList* dir)
-{
-  // 
-  // Define the output histograms.  These are put in a sub list of the
-  // passed list.   The histograms are merged before the parent task calls 
-  // AliAnalysisTaskSE::Terminate 
-  // 
-  // Parameters:
-  //    dir Directory to add to 
-  //
-  DGUARD(fDebug, 1, "Define output in AliFMDEnergyFitter");
-  TList* d = new TList;
-  d->SetName(GetName());
-  d->SetOwner(true);
-  dir->Add(d);
-  
-  d->Add(&fEtaAxis);
-  d->Add(AliForwardUtil::MakeParameter("lowCut",        fLowCut));
-  d->Add(AliForwardUtil::MakeParameter("nParticles",    fNParticles));
-  d->Add(AliForwardUtil::MakeParameter("minEntries",    fMinEntries));
-  d->Add(AliForwardUtil::MakeParameter("subtractBins",  fFitRangeBinWidth));
-  d->Add(AliForwardUtil::MakeParameter("doFits",        fDoFits));
-  d->Add(AliForwardUtil::MakeParameter("doObject",      fDoMakeObject));
-  d->Add(AliForwardUtil::MakeParameter("maxE",          fMaxE));
-  d->Add(AliForwardUtil::MakeParameter("nEbins",        fNEbins));
-  d->Add(AliForwardUtil::MakeParameter("increasingBins",fUseIncreasingBins));
-  d->Add(AliForwardUtil::MakeParameter("maxRelPerError",fMaxRelParError));
-  d->Add(AliForwardUtil::MakeParameter("maxChi2PerNDF", fMaxChi2PerNDF));
-  d->Add(AliForwardUtil::MakeParameter("minWeight",     fMinWeight));
-  d->Add(AliForwardUtil::MakeParameter("regCut",        fRegularizationCut));
-  
-  TIter    next(&fRingHistos);
-  RingHistos* o = 0;
-  while ((o = static_cast<RingHistos*>(next()))) {
-    o->CreateOutputObjects(d);
-  }
-}
-//____________________________________________________________________
-void
 AliFMDEnergyFitter::SetDebug(Int_t dbg)
 {
   // 
@@ -521,7 +477,7 @@ AliFMDEnergyFitter::ReadParameters(const TCollection* col)
   else {
     if (axis->GetXbins()->GetArray()) 
       fEtaAxis.Set(axis->GetNbins(), axis->GetXbins()->GetArray());
-    else 
+   else 
       fEtaAxis.Set(axis->GetNbins(), axis->GetXmin(), axis->GetXmax());
   } 
   GetParam(ret,col,"lowCut",        fLowCut);
@@ -598,7 +554,8 @@ AliFMDEnergyFitter::RingHistos::RingHistos()
   : AliForwardUtil::RingHistos(), 
     fEDist(0), 
     fEmpty(0),
-    fEtaEDists(0), 
+    // fEtaEDists(0), 
+    fHist(0),
     fList(0),
     fBest(0),
     fFits("AliFMDCorrELossFit::ELossFit", 200),
@@ -616,7 +573,8 @@ AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
   : AliForwardUtil::RingHistos(d,r), 
     fEDist(0), 
     fEmpty(0),
-    fEtaEDists(0), 
+    // fEtaEDists(0), 
+    fHist(0),
     fList(0),
     fBest(0),
     fFits("AliFMDCorrELossFit::ELossFit", 200),
@@ -634,102 +592,6 @@ AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
   fBest.Expand(0);
 }
 //____________________________________________________________________
-AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
-  : AliForwardUtil::RingHistos(o), 
-    fEDist(o.fEDist), 
-    fEmpty(o.fEmpty),
-    fEtaEDists(0), 
-    fList(0),
-    fBest(0),
-    fFits("AliFMDCorrELossFit::ELossFit", 200),
-    fDebug(0)
-{
-  // 
-  // Copy constructor 
-  // 
-  // Parameters:
-  //    o Object to copy from 
-  //
-  DGUARD(fDebug, 3, "Copy CTOR AliFMDEnergyFitter::RingHistos");
-  fFits.Clear();
-  if (o.fEtaEDists) {
-    fEtaEDists = new TList;
-    fEtaEDists->SetOwner();
-    fEtaEDists->SetName(o.fEtaEDists->GetName());
-    TIter next(o.fEtaEDists);
-    TObject* obj = 0;
-    while ((obj = next())) fEtaEDists->Add(obj->Clone());
-  }
-  if (o.fList) {
-    fList = new TList;
-    fList->SetOwner();
-    fList->SetName(fName);
-    TIter next(o.fList);
-    TObject* obj = 0;
-    while ((obj = next())) {
-      if (obj == o.fEtaEDists) {
-       fList->Add(fEtaEDists);
-       continue;
-      }
-      fList->Add(obj->Clone());
-    }
-  }
-}
-
-//____________________________________________________________________
-AliFMDEnergyFitter::RingHistos&
-AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
-{
-  // 
-  // Assignment operator 
-  // 
-  // Parameters:
-  //    o Object to assign from 
-  // 
-  // Return:
-  //    Reference to this 
-  //
-  fDebug = o.fDebug;
-  DGUARD(fDebug, 3, "Assignment of AliFMDEnergyFitter::RingHistos");
-  if (&o == this) return *this; 
-  AliForwardUtil::RingHistos::operator=(o);
-  
-  if (fEDist)     delete fEDist;
-  if (fEmpty)     delete fEmpty;
-  if (fEtaEDists) fEtaEDists->Delete();
-  if (fList)      fList->Delete();
-  fFits.Clear();
-
-  fEDist = static_cast<TH1D*>(o.fEDist->Clone());
-  fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
-  
-  if (o.fEtaEDists) {
-    fEtaEDists = new TList;
-    fEtaEDists->SetOwner();
-    fEtaEDists->SetName(o.fEtaEDists->GetName());
-    TIter next(o.fEtaEDists);
-    TObject* obj = 0;
-    while ((obj = next())) fEtaEDists->Add(obj->Clone());
-  }
-
-  if (o.fList) {
-    fList = new TList;
-    fList->SetOwner();
-    fList->SetName(fName);
-    TIter next(o.fList);
-    TObject* obj = 0;
-    while ((obj = next())) {
-      if (obj == o.fEtaEDists) {
-       fList->Add(fEtaEDists);
-       continue;
-      }
-      fList->Add(obj->Clone());
-    }
-  }
-
-  return *this;
-}
-//____________________________________________________________________
 AliFMDEnergyFitter::RingHistos::~RingHistos()
 {
   // 
@@ -739,43 +601,6 @@ AliFMDEnergyFitter::RingHistos::~RingHistos()
   // if (fEDist) delete fEDist;
   // fEtaEDists.Delete();
 }
-
-//____________________________________________________________________
-void
-AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, 
-                                    Int_t /* icent */,  Double_t mult)
-{
-  // 
-  // Fill histogram 
-  // 
-  // Parameters:
-  //    empty  True if event is empty
-  //    ieta   Eta bin (0-based)
-  //    icent  Centrality bin (1-based)
-  //    mult   Signal 
-  //
-  DGUARD(fDebug, 10, "Filling in AliFMDEnergyFitter::RingHistos");
-  if (empty) { 
-    fEmpty->Fill(mult);
-    return;
-  }
-  fEDist->Fill(mult);
-
-  if (!fEtaEDists) { 
-    Warning("Fill", "No list of E dists defined");
-    return;
-  }
-  // Here, we should first look up in a centrality array, and then in
-  // that array look up the eta bin - or perhaps we should do it the
-  // other way around?
-  if (ieta < 0 || ieta >= fEtaEDists->GetEntries()) return;
-  
-  TH1D* h = static_cast<TH1D*>(fEtaEDists->At(ieta));
-  if (!h) return;
-
-  h->Fill(mult);
-}
-
 //__________________________________________________________________
 TArrayD
 AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min, 
@@ -808,50 +633,181 @@ AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
 }
 
 //____________________________________________________________________
-void
-AliFMDEnergyFitter::RingHistos::Make(Int_t    ieta, 
-                                    Double_t emin, 
-                                    Double_t emax,
-                                    Double_t deMax, 
-                                    Int_t    nDeBins, 
-                                    Bool_t   incr)
+TH2*
+AliFMDEnergyFitter::RingHistos::Make(const char*  name, 
+                                    const char*  title, 
+                                    const TAxis& eAxis, 
+                                    Double_t     deMax, 
+                                    Int_t        nDeBins, 
+                                    Bool_t       incr)
 {
   // 
   // Make E/E_mip histogram 
   // 
   // Parameters:
-  //    ieta    Eta bin
-  //    eMin    Least signal
-  //    eMax    Largest signal 
   //    deMax   Maximum energy loss 
   //    nDeBins Number energy loss bins 
   //    incr    Whether to make bins of increasing size
   //
-  TH1D* h = 0;
-  TString name  = Form(fgkEDistFormat, fName.Data(), ieta);
-  TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
-                      "(#eta bin %d)", fName.Data(), emin, emax, ieta);
-  if (!incr) 
-    h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
-  else { 
+  TH2* h = 0;
+  TAxis mAxis;
+  if (incr) {
     TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
-    h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
+    mAxis.Set(deAxis.GetSize()-1, deAxis.GetArray());
+  }
+  else 
+    mAxis.Set(nDeBins, 0, deMax);
+  
+  if (mAxis.GetXbins()->GetArray()) {
+    // Variable bin since in Delta 
+    if (eAxis.GetXbins()->GetArray()) {
+      // Variadic bin size in eta 
+      h = new TH2D(name, title, 
+                  eAxis.GetNbins(), eAxis.GetXbins()->GetArray(),
+                  mAxis.GetNbins(), mAxis.GetXbins()->GetArray());
+    }
+    else { 
+      // Evenly spaced bins in eta
+      h = new TH2D(name, title, 
+                  eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(), 
+                  mAxis.GetNbins(), mAxis.GetXbins()->GetArray());
+    }
+  }
+  else { 
+    // Make increasing bins axis 
+    if (eAxis.GetXbins()->GetArray()) {
+      // Variable size eta bins 
+      h = new TH2D(name, title, 
+                  eAxis.GetNbins(), eAxis.GetXbins()->GetArray(),
+                  mAxis.GetNbins(), mAxis.GetXmin(), mAxis.GetXmax());
+    }
+    else {
+      // Fixed size eta bins 
+      h = new TH2D(name, title, 
+                  eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(), 
+                  mAxis.GetNbins(), mAxis.GetXmin(), mAxis.GetXmax());
+    }
   }
-    
   h->SetDirectory(0);
-  h->SetXTitle("#DeltaE/#DeltaE_{mip}");       
+  h->SetYTitle("#sum#DeltaE/#DeltaE_{mip}");   
+  h->SetXTitle("#eta");        
   h->SetFillColor(Color());
   h->SetMarkerColor(Color());
   h->SetLineColor(Color());
   h->SetFillStyle(3001);
+  h->SetMarkerStyle(20);
   h->Sumw2();
+
+  return h;
+}
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::RingHistos::CreateOutputObjects(TList* dir)
+{
+  // 
+  // Define outputs
+  // 
+  // Parameters:
+  //    dir 
+  //
+  DGUARD(fDebug, 2, "Define output in AliFMDEnergyFitter::RingHistos");
+  fList = DefineOutputList(dir);
+  
+  // fEtaEDists = new TList;
+  // fEtaEDists->SetOwner();
+  // fEtaEDists->SetName("EDists");
+  // 
+  // fList->Add(fEtaEDists);
+}
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::RingHistos::SetupForData(const TAxis& eAxis, 
+                                            const TAxis& /* cAxis */,
+                                            Double_t     maxDE, 
+                                            Int_t        nDEbins, 
+                                            Bool_t       useIncrBin)
+{
+  // 
+  // Initialise object 
+  // 
+  // Parameters:
+  //    eAxis      Eta axis
+  //    maxDE      Max energy loss to consider 
+  //    nDEbins    Number of bins 
+  //    useIncrBin Whether to use an increasing bin size 
+  //
+  DGUARD(fDebug, 2, "Initialize in AliFMDEnergyFitter::RingHistos");
+  fEDist = new TH1D(Form("%s_edist", fName.Data()), 
+                   Form("#sum#DeltaE/#DeltaE_{mip} for %s", fName.Data()), 
+                   200, 0, 6);
+  fEDist->SetXTitle("#sum#Delta/#Delta_{mip}");
+  fEDist->SetFillColor(Color());
+  fEDist->SetLineColor(Color());
+  fEDist->SetMarkerColor(Color());
+  fEDist->SetFillStyle(3001);
+  fEDist->SetMarkerStyle(20);
+  fEDist->Sumw2();
+  fEDist->SetDirectory(0);
+
+  fEmpty = static_cast<TH1D*>(fEDist->Clone(Form("%s_empty", fName.Data())));
+  fEmpty->SetTitle(Form("#sum#DeltaE/#DeltaE_{mip} for %s (empty events)", 
+                       fName.Data()));
+  fEmpty->SetDirectory(0);
+
+  fList->Add(fEDist);
+  fList->Add(fEmpty);
+  fHist = Make("eloss", "#sum#Delta/#Delta_{mip}", eAxis, 
+              maxDE, nDEbins, useIncrBin);
+  fList->Add(fHist);
+  // fList->ls();
+  // fEtaEDists.ls();
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::RingHistos::Fill(Bool_t   empty, 
+                                    Double_t eta, 
+                                    Int_t /* icent */,  
+                                    Double_t mult)
+{
+  // 
+  // Fill histogram 
+  // 
+  // Parameters:
+  //    empty  True if event is empty
+  //    ieta   Eta bin (0-based)
+  //    icent  Centrality bin (1-based)
+  //    mult   Signal 
+  //
+  DGUARD(fDebug, 10, "Filling in AliFMDEnergyFitter::RingHistos");
+  if (empty) { 
+    fEmpty->Fill(mult);
+    return;
+  }
+  fEDist->Fill(mult);
+
+  // if (!fEtaEDists) { 
+  if (!fHist) {
+    Warning("Fill", "No list of E dists defined");
+    return;
+  }
+  fHist->Fill(eta, mult);
+  // Here, we should first look up in a centrality array, and then in
+  // that array look up the eta bin - or perhaps we should do it the
+  // other way around?
+  // if (ieta < 0 || ieta >= fEtaEDists->GetEntries()) return;
   
-  fEtaEDists->AddAt(h, ieta-1);
+  // TH1D* h = static_cast<TH1D*>(fEtaEDists->At(ieta));
+  // if (!h) return;
+
+  // h->Fill(mult);
 }
+
 //____________________________________________________________________
-TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name, 
-                                             const char* title, 
-                                             const TAxis& eta) const
+TH1* 
+AliFMDEnergyFitter::RingHistos::MakePar(const char* name, 
+                                       const char* title, 
+                                       const TAxis& eta) const
 {
   // 
   // Make a parameter histogram
@@ -879,7 +835,7 @@ TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name,
   return h;
 }
 //____________________________________________________________________
-TH1D*
+TH1*
 AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name, 
                                          const char* title, 
                                          const TAxis& eta, 
@@ -924,95 +880,8 @@ AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name,
 }
                     
 //____________________________________________________________________
-void
-AliFMDEnergyFitter::RingHistos::SetupForData(const TAxis& eAxis, 
-                                    const TAxis& /* cAxis */,
-                                    Double_t maxDE, Int_t nDEbins, 
-                                    Bool_t useIncrBin)
-{
-  // 
-  // Initialise object 
-  // 
-  // Parameters:
-  //    eAxis      Eta axis
-  //    maxDE      Max energy loss to consider 
-  //    nDEbins    Number of bins 
-  //    useIncrBin Whether to use an increasing bin size 
-  //
-  DGUARD(fDebug, 2, "Initialize in AliFMDEnergyFitter::RingHistos");
-  fEDist = new TH1D(Form("%s_edist", fName.Data()), 
-                   Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()), 
-                   200, 0, 6);
-  fEmpty = new TH1D(Form("%s_empty", fName.Data()), 
-                   Form("#DeltaE/#DeltaE_{mip} for %s (empty events)", 
-                        fName.Data()), 200, 0, 6);
-  fList->Add(fEDist);
-  fList->Add(fEmpty);
-  // Here, we should make an array of centrality bins ...
-  // In that case, the structure will be 
-  // 
-  //    -+- Ring1 -+- Centrality1 -+- Eta1
-  //     |         |               +- Eta2
-  //     ...       ...             ...
-  //     |         |               +- EtaM
-  //     |         +- Centrality2 -+- Eta1
-  //     ...       ...             ...
-  //     |         |               +- EtaM
-  //     ...       ...
-  //     |         +- CentralityN -+- Eta1
-  //     ...       ...             ...
-  //     |                         +- EtaM
-  //    -+- Ring2 -+- Centrality1 -+- Eta1
-  //     |         |               +- Eta2
-  //     ...       ...             ...
-  //     |         |               +- EtaM
-  //     |         +- Centrality2 -+- Eta1
-  //     ...       ...             ...
-  //     |         |               +- EtaM
-  //     ...       ...
-  //     |         +- CentralityN -+- Eta1
-  //     ...       ...             ...
-  //     |                         +- EtaM
-  //     ...       ...             ...
-  //    -+- RingO -+- Centrality1 -+- Eta1
-  //               |               +- Eta2
-  //               ...             ...
-  //               |               +- EtaM
-  //               +- Centrality2 -+- Eta1
-  //               ...             ...
-  //               |               +- EtaM
-  //               ...
-  //               +- CentralityN -+- Eta1
-  //               ...             ...
-  //                               +- EtaM
-  // 
-  // fEtaEDists.Expand(eAxis.GetNbins());
-  for (Int_t i = 1; i <= eAxis.GetNbins(); i++) { 
-    Double_t min = eAxis.GetBinLowEdge(i);
-    Double_t max = eAxis.GetBinUpEdge(i);
-    // Or may we should do it here?  In that case we'd have the
-    // following structure:
-    //     Ring1 -+- Eta1 -+- Centrality1 
-    //            |        +- Centrality2
-    //            ...      ...
-    //            |        +- CentralityN
-    //            +- Eta2 -+- Centrality1 
-    //            |        +- Centrality2
-    //            ...      ...
-    //            |        +- CentralityN
-    //            ...
-    //            +- EtaM -+- Centrality1 
-    //                     +- Centrality2
-    //                     ...
-    //                     +- CentralityN
-    Make(i, min, max, maxDE, nDEbins, useIncrBin);
-  }
-  // fEtaEDists.ls();
-}
-//____________________________________________________________________
 TObjArray*
 AliFMDEnergyFitter::RingHistos::Fit(TList*           dir, 
-                                   const TAxis&     eta,
                                    Double_t         lowCut, 
                                    UShort_t         nParticles,
                                    UShort_t         minEntries,
@@ -1023,6 +892,27 @@ AliFMDEnergyFitter::RingHistos::Fit(TList*           dir,
                                    Double_t         regCut,
                                    EResidualMethod  residuals) const
 {
+  return FitSlices(dir, "eloss", lowCut, nParticles, minEntries, 
+                  minusBins, relErrorCut, chi2nuCut, minWeight, regCut, 
+                  residuals, true, &fBest);
+}
+
+//____________________________________________________________________
+TObjArray*
+AliFMDEnergyFitter::RingHistos::FitSlices(TList*           dir, 
+                                         const char*      name, 
+                                         Double_t         lowCut, 
+                                         UShort_t         nParticles,
+                                         UShort_t         minEntries,
+                                         UShort_t         minusBins, 
+                                         Double_t         relErrorCut, 
+                                         Double_t         chi2nuCut,
+                                         Double_t         minWeight,
+                                         Double_t         regCut,
+                                         EResidualMethod  residuals,
+                                         Bool_t           scaleToPeak,
+                                         TObjArray*       best) const
+{
   // 
   // Fit each histogram to up to @a nParticles particle responses.
   // 
@@ -1047,26 +937,39 @@ AliFMDEnergyFitter::RingHistos::Fit(TList*           dir,
   if (!l) return 0; 
 
   // Get the energy distributions from the output container 
-  TList* dists = static_cast<TList*>(l->FindObject("EDists"));
-  if (!dists) { 
-    AliWarning(Form("Didn't find EtaEDists (%s) in %s", 
-                   fName.Data(), l->GetName()));
+  // TList* dists = static_cast<TList*>(l->FindObject("EDists"));
+  // if (!dists) { 
+  //   AliWarning(Form("Didn't find EtaEDists (%s) in %s", 
+  //               fName.Data(), l->GetName()));
+  //   l->ls();
+  //   return 0;
+  // }
+  TH2* h = static_cast<TH2*>(l->FindObject(name));
+  if (!h) { 
+    AliWarningF("Didn't find 2D histogram '%s' in %s", name, l->GetName());
     l->ls();
     return 0;
   }
+  const TAxis& eta = *(h->GetXaxis());
+
+  // Create an output list for the fitted distributions 
+  TList* out = new TList;
+  out->SetOwner();
+  out->SetName(Form("%sDists", name));
+  l->Add(out);
 
   // Optional container for residuals 
   TList* resi = 0;
   if (residuals != kNoResiduals) {
     resi = new TList();
-    resi->SetName("residuals");
+    resi->SetName(Form("%sResiduals", name));
     resi->SetOwner();
     l->Add(resi);
   }
 
   // Container of the fit results histograms 
   TObjArray* pars  = new TObjArray(3+nParticles+1);
-  pars->SetName("FitResults");
+  pars->SetName(Form("%sResults", name));
   l->Add(pars);
 
   // Result objects stored in separate list on the output 
@@ -1089,64 +992,72 @@ AliFMDEnergyFitter::RingHistos::Fit(TList*           dir,
     pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
 
   
-  Int_t nDists = dists->GetEntries();
+  Int_t nDists = h->GetNbinsX(); // dists->GetEntries();
   Int_t low    = nDists;
   Int_t high   = 0;
   Int_t nEmpty = 0;
   Int_t nLow   = 0;
   Int_t nFitted= 0;
-  fBest.Expand(nDists+1);
-  fBest.Clear();
-  fBest.SetOwner(false);
+  if (best) {
+    best->Expand(nDists+1);
+    best->Clear();
+    best->SetOwner(false);
+  }
   for (Int_t i = 0; i < nDists; i++) { 
-    TH1D* dist = static_cast<TH1D*>(dists->At(i));
+    // TH1D* dist = static_cast<TH1D*>(dists->At(i));
     // Ignore empty histograms altoghether 
-    if (!dist || dist->GetEntries() <= 0) { 
+    Int_t b    = i+1;
+    TH1D* dist = h->ProjectionY(Form(fgkEDistFormat,GetName(),b),b,b,"e");
+    if (!dist) { 
+      // If we got the null pointer, return 0
       nEmpty++;
       continue;
     }
+    // Then releasing the histogram from the it's directory
+    dist->SetDirectory(0);
+    // Set a meaningful title
+    dist->SetTitle(Form("#Delta/#Delta_{mip} for %s in %6.2f<#eta<%6.2f",
+                       GetName(), eta.GetBinLowEdge(b),
+                       eta.GetBinUpEdge(b)));
 
-    // Scale to the bin-width
-    dist->Scale(1., "width");
-    
-    // Narrow search window for the peak 
-    Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(lowCut),3);
-    Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(10),
-                                 dist->GetNbinsX());
-    dist->GetXaxis()->SetRange(cutBin, maxBin);
-    
-    // Get the bin with maximum 
-    Int_t    peakBin = dist->GetMaximumBin();
-    
-    // Normalize peak to 1 
-    // Double_t max = dist->GetMaximum(); 
-    Double_t max = dist->GetBinContent(peakBin); // Maximum(); 
-    if (max <= 0) continue;
-    dist->Scale(1/max);
-    
-    // Check that we have enough entries 
-    Int_t nEntries = Int_t(dist->GetEntries());
-    if (nEntries <= minEntries) { 
-      AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)",
-                     i, dist->GetName(), nEntries, minEntries));
-      nLow++;
+    // Now fit 
+    UShort_t    status1 = 0;
+    ELossFit_t* res     = FitHist(dist,
+                                 lowCut, 
+                                 nParticles,
+                                 minEntries,
+                                 minusBins,   
+                                 relErrorCut,
+                                 chi2nuCut,
+                                 minWeight,
+                                 regCut,
+                                 scaleToPeak,
+                                 status1);
+    if (!res) {
+      switch (status1) { 
+      case 1: nEmpty++; break;
+      case 2: nLow++;   break;
+      }
+      delete dist;
       continue;
     }
-
-    // Now fit 
-    ELossFit_t* res = FitHist(dist,i+1,lowCut, nParticles,minusBins,
-                             relErrorCut,chi2nuCut,minWeight,regCut);
-    if (!res) continue;
+      
+    // Now count as fitted, store as best fits, and add distribution
+    // to the dedicated output list
     nFitted++;
-    fBest.AddAt(res, i+1);
+    res->fBin = b; // We only have the bin information here 
+    if (best) best->AddAt(res, b);
+    out->Add(dist);
     // dist->GetListOfFunctions()->Add(res);
     
+    // If asked to calculate residuals, do so, and store result on the
+    // dedicated output list
     if (residuals != kNoResiduals && resi) 
       CalculateResiduals(residuals, lowCut, dist, res, resi);
 
     // Store eta limits 
-    low   = TMath::Min(low,i+1);
-    high  = TMath::Max(high,i+1);
+    low   = TMath::Min(low,b);
+    high  = TMath::Max(high,b);
 
     // Get the reduced chi square
     Double_t chi2 = res->fChi2; // GetChisquare();
@@ -1154,24 +1065,24 @@ AliFMDEnergyFitter::RingHistos::Fit(TList*           dir,
     
     // Store results of best fit in output histograms 
     // res->SetLineWidth(4);
-    hChi2   ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
-    hC      ->SetBinContent(i+1, res->fC); // GetParameter(kC));   
-    hDelta  ->SetBinContent(i+1, res->fDelta); // GetParameter(kDelta)); 
-    hXi     ->SetBinContent(i+1, res->fXi); // GetParameter(kXi));   
-    hSigma  ->SetBinContent(i+1, res->fSigma); // GetParameter(kSigma));   
-    hSigmaN ->SetBinContent(i+1, res->fSigmaN); // GetParameter(kSigmaN));   
-    hN      ->SetBinContent(i+1, res->fN); // GetParameter(kN));   
-
-    hC     ->SetBinError(i+1, res->fEC); // GetParError(kC));
-    hDelta ->SetBinError(i+1, res->fEDelta); // GetParError(kDelta));
-    hXi    ->SetBinError(i+1, res->fEXi); // GetParError(kXi));
-    hSigma ->SetBinError(i+1, res->fESigma); // GetParError(kSigma));
-    hSigmaN->SetBinError(i+1, res->fESigmaN); // GetParError(kSigmaN));
-    // hN     ->SetBinError(i+1, res->fGetParError(kN));
+    hChi2   ->SetBinContent(b, ndf > 0 ? chi2 / ndf : 0);
+    hC      ->SetBinContent(b, res->fC); 
+    hDelta  ->SetBinContent(b, res->fDelta); 
+    hXi     ->SetBinContent(b, res->fXi); 
+    hSigma  ->SetBinContent(b, res->fSigma); 
+    hSigmaN ->SetBinContent(b, res->fSigmaN); 
+    hN      ->SetBinContent(b, res->fN); 
+
+    hC     ->SetBinError(b, res->fEC); 
+    hDelta ->SetBinError(b, res->fEDelta);
+    hXi    ->SetBinError(b, res->fEXi); 
+    hSigma ->SetBinError(b, res->fESigma); 
+    hSigmaN->SetBinError(b, res->fESigmaN); 
+    // hN     ->SetBinError(b, res->fGetParError(kN));
 
     for (UShort_t j = 0; j < nParticles-1; j++) {
-      hA[j]->SetBinContent(i+1, res->fA[j]); // GetParameter(kA+j));
-      hA[j]->SetBinError(i+1, res->fEA[j]); // GetParError(kA+j));
+      hA[j]->SetBinContent(b, res->fA[j]); 
+      hA[j]->SetBinError(b, res->fEA[j]); 
     }
   }
   printf("%s: Out of %d histograms, %d where empty, %d had too little data,"
@@ -1179,18 +1090,20 @@ AliFMDEnergyFitter::RingHistos::Fit(TList*           dir,
         GetName(), nDists, nEmpty, nLow, nDists-nEmpty-nLow, nFitted);
 
   // Fit the full-ring histogram 
-  TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
-  if (total && total->GetEntries() >= minEntries) { 
-
-    // Scale to the bin-width
-    total->Scale(1., "width");
-    
-    // Normalize peak to 1 
-    Double_t max = total->GetMaximum(); 
-    if (max > 0) total->Scale(1/max);
-    
-    ELossFit_t* resT = FitHist(total,nDists+1,lowCut, nParticles,minusBins,
-                              relErrorCut,chi2nuCut,minWeight,regCut);
+  TH1*        total   = GetOutputHist(l, Form("%s_edist", fName.Data()));
+  if (total) {
+    UShort_t    statusT = 0;
+    ELossFit_t* resT    = FitHist(total,
+                                 lowCut, 
+                                 nParticles,
+                                 minEntries, 
+                                 minusBins,
+                                 relErrorCut,
+                                 chi2nuCut,
+                                 minWeight,
+                                 regCut,
+                                 scaleToPeak,
+                                 statusT);
     if (resT) { 
       // Make histograms for the result of this fit 
       Double_t chi2 = resT->GetChi2();
@@ -1216,7 +1129,7 @@ AliFMDEnergyFitter::RingHistos::Fit(TList*           dir,
     }
   }
 
-  TH1* status = new TH1I("status", "Status of Fits", 5, 0, 5);
+  TH1* status = new TH1I(Form("%sStatus",name), "Status of Fits", 5, 0, 5);
   status->GetXaxis()->SetBinLabel(1, "Total");
   status->GetXaxis()->SetBinLabel(2, "Empty");
   status->GetXaxis()->SetBinLabel(3, Form("<%d", minEntries));
@@ -1235,48 +1148,23 @@ AliFMDEnergyFitter::RingHistos::Fit(TList*           dir,
   status->SetDirectory(0);
   status->SetStats(0);
   pars->Add(status);
-    
-  // Clean up list of histogram.  Histograms with no entries or 
-  // no functions are deleted.  We have to do this using the TObjLink 
-  // objects stored in the list since ROOT cannot guaranty the validity 
-  // of iterators when removing from a list - tsck.  Should just implement
-  // TIter::Remove(). 
-  TObjLink* lnk = dists->FirstLink();
-  while (lnk) {
-    TH1* h = static_cast<TH1*>(lnk->GetObject());
-    bool remove = false;
-    if (h->GetEntries() <= 0) { 
-      // AliWarning(Form("No entries in %s - removing", h->GetName()));
-      remove = true;
-    }
-    else if (h->GetListOfFunctions()->GetEntries() <= 0) {
-      // AliWarning(Form("No fuctions associated with %s - removing",
-      //            h->GetName()));
-      remove = true;
-    }
-    if (remove) {
-      TObjLink* keep = lnk->Next();
-      dists->Remove(lnk);
-      lnk = keep;
-      continue;
-    }
-    lnk = lnk->Next();
-  }
-
   return pars;
 }
 
+
 //____________________________________________________________________
 AliFMDEnergyFitter::RingHistos::ELossFit_t*
-AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
-                                       UShort_t bin, 
-                                       Double_t lowCut, 
-                                       UShort_t nParticles, 
-                                       UShort_t minusBins, 
-                                       Double_t relErrorCut, 
-                                       Double_t chi2nuCut,
-                                       Double_t minWeight,
-                                       Double_t regCut) const
+AliFMDEnergyFitter::RingHistos::FitHist(TH1*      dist,
+                                       Double_t  lowCut, 
+                                       UShort_t  nParticles, 
+                                       UShort_t  minEntries,
+                                       UShort_t  minusBins, 
+                                       Double_t  relErrorCut, 
+                                       Double_t  chi2nuCut,
+                                       Double_t  minWeight,
+                                       Double_t  regCut,
+                                       Bool_t    scaleToPeak,
+                                       UShort_t& status) const
 {
   // 
   // Fit a signal histogram.  First, the bin @f$ b_{min}@f$ with
@@ -1303,10 +1191,47 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
   // Return:
   //    The best fit function 
   //
-  DGUARD(fDebug, 3, "Fit histogram in AliFMDEnergyFitter::RingHistos: %s",
+  DGUARD(fDebug, 2, "Fit histogram in AliFMDEnergyFitter::RingHistos: %s",
         dist->GetName());
   Double_t maxRange = 10;
 
+
+  if (dist->GetEntries() <= 0) { 
+    status = 1; // `empty'
+    return 0;
+  }
+
+  // Scale to the bin-width
+  dist->Scale(1., "width");
+    
+  // Narrow search window for the peak 
+  Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(lowCut),3);
+  Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(10),
+                               dist->GetNbinsX());
+  dist->GetXaxis()->SetRange(cutBin, maxBin);
+    
+  // Get the bin with maximum 
+  Int_t    peakBin = dist->GetMaximumBin();
+    
+  // Normalize peak to 1 
+  // Double_t max = dist->GetMaximum(); 
+  Double_t max = dist->GetBinContent(peakBin); // Maximum(); 
+  if (max <= 0) {
+    status = 1; // `empty'
+    return 0;
+  }
+  if (scaleToPeak) dist->Scale(1/max);
+  DMSG(fDebug,5,"max(%s) -> %f", dist->GetName(), max);
+
+  // Check that we have enough entries 
+  Int_t nEntries = Int_t(dist->GetEntries());
+  if (nEntries <= minEntries) { 
+    AliWarning(Form("Histogram at %s has too few entries (%d <= %d)",
+                   dist->GetName(), nEntries, minEntries));
+    status = 2;
+    return 0;
+  }
+
   // Create a fitter object 
   AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins); 
   f.Clear();
@@ -1326,16 +1251,21 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
   // no matter what. 
   if (nParticles == 1) {
     TF1* r = f.Fit1Particle(dist, 0);
-    if (!r) return 0;
+    if (!r) {
+      status = 3; // No-fit
+      return 0;
+    }
     TF1* ff = new TF1(*r);
     dist->GetListOfFunctions()->Add(ff);
     ELossFit_t* ret = new ELossFit_t(0, *ff);
     ret->CalculateQuality(chi2nuCut, relErrorCut, minWeight);
+    status = 0; // OK
     return ret;
   }
 
   // Fit from 2 upto n particles  
   for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0);
+
   // Now, we need to select the best fit 
   Int_t nFits = f.GetFitResults().GetEntriesFast();
   for (Int_t i = nFits-1; i >= 0; i--) { 
@@ -1343,87 +1273,21 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
     // ff->SetRange(0, maxRange);
     dist->GetListOfFunctions()->Add(new TF1(*ff));
   }
+  status = 0; // OK
 
   // Here, we use the real quality assesor instead of the old
   // `CheckResult' to ensure consitency in all output.
-  ELossFit_t* ret = FindBestFit(bin, dist, relErrorCut, chi2nuCut, minWeight);
+  ELossFit_t* ret = FindBestFit(dist, relErrorCut, chi2nuCut, minWeight);
+  if (!ret) status = 3;
   return ret;
 }
 
 //__________________________________________________________________
-void
-AliFMDEnergyFitter::RingHistos::FindBestFits(const TList*        d, 
-                                            AliFMDCorrELossFit& obj,
-                                            const TAxis&        eta,     
-                                            Double_t           ,//relErrorCut, 
-                                            Double_t           ,//chi2nuCut,
-                                            Double_t           )//minWeightCut)
-{
-  // 
-  // Find the best fits 
-  // 
-  // Parameters:
-  //    d              Parent list
-  //    obj            Object to add fits to
-  //    eta            Eta axis 
-  //    relErrorCut    Cut applied to relative error of parameter. 
-  //                       Note, for multi-particle weights, the cut 
-  //                       is loosend by a factor of 2 
-  //    chi2nuCut      Cut on @f$ \chi^2/\nu@f$ - 
-  //                       the reduced @f$\chi^2@f$ 
-  //    minWeightCut   Least valid @f$ a_i@f$ 
-  //
-
-  // Get our ring list from the output container 
-  TList* l = GetOutputList(d);
-  if (!l) return; 
-
-  // Get the energy distributions from the output container 
-  TList* dists = static_cast<TList*>(l->FindObject("EDists"));
-  if (!dists) { 
-    AliWarning(Form("Didn't find %s_EtaEDists in %s", 
-                   fName.Data(), l->GetName()));
-    l->ls();
-    return;
-  }
-  Int_t nBin = eta.GetNbins();
-  if (fBest.GetEntriesFast() <= 0) { 
-    AliWarningF("No fits found for %s", GetName());
-    return;
-  }
-  
-  for (Int_t b = 1; b <= nBin; b++) { 
-    TString n(Form(fgkEDistFormat, fName.Data(), b));
-    TH1D*   dist = static_cast<TH1D*>(dists->FindObject(n));
-    // Ignore empty histograms altoghether 
-    if (!dist || dist->GetEntries() <= 0) continue; 
-    
-    ELossFit_t* best = static_cast<ELossFit_t*>(fBest.At(b));
-    if (!best) { 
-      AliErrorF("No best fit found @ %d for %s", b, GetName());
-      continue;
-    }
-    // FindBestFit(b, dist, relErrorCut, chi2nuCut, minWeightCut);
-    best->fDet  = fDet; 
-    best->fRing = fRing;
-    best->fBin  = b; // 
-    if (fDebug > 0) {
-      printf("Bin # %3d: ", b); 
-      best->Print("s");
-    }
-    // Double_t eta = fAxis->GetBinCenter(b);
-    obj.SetFit(fDet, fRing, b, best); 
-    // new ELossFit_t(*best));
-  }
-}
-
-//__________________________________________________________________
 AliFMDEnergyFitter::RingHistos::ELossFit_t* 
-AliFMDEnergyFitter::RingHistos::FindBestFit(UShort_t b, 
-                                           const TH1* dist,
-                                           Double_t relErrorCut, 
-                                           Double_t chi2nuCut,
-                                           Double_t minWeightCut)  const
+AliFMDEnergyFitter::RingHistos::FindBestFit(const TH1* dist,
+                                           Double_t   relErrorCut, 
+                                           Double_t   chi2nuCut,
+                                           Double_t   minWeightCut)  const
 {
   // 
   // Find the best fit 
@@ -1438,37 +1302,43 @@ AliFMDEnergyFitter::RingHistos::FindBestFit(UShort_t b,
   //    minWeightCut   Least valid @f$ a_i@f$ 
   // 
   // Return:
-  //    Best fit 
+  //    Best fit or null
   //
   TList* funcs = dist->GetListOfFunctions();
+  TF1*   func  = 0;
+  Int_t  i     = 0;
   TIter  next(funcs);
-  TF1*   func = 0;
-  fFits.Clear();
-  Int_t  i = 0;
-  if (fDebug) 
-    printf("Find best fit for %s ... ", dist->GetName());
-  // Info("FindBestFit", "%s", dist->GetName());
+  fFits.Clear(); // This is only ever used here
+
+  if (fDebug) printf("Find best fit for %s ... ", dist->GetName());
+  if (fDebug > 2) printf("\n");
+
+  // Loop over all functions stored in distribution, 
+  // and calculate the quality 
   while ((func = static_cast<TF1*>(next()))) { 
     ELossFit_t* fit = new(fFits[i++]) ELossFit_t(0,*func);
     fit->fDet  = fDet;
     fit->fRing = fRing;
-    fit->fBin  = b;
+    // fit->fBin  = b;
     fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
     if (fDebug > 2) 
       Printf("%10s: %3d (chi^2/nu: %6.3f)", 
             func->GetName(), fit->fQuality, 
             (fit->fNu > 0 ? fit->fChi2 / fit->fNu : 999));
   }
+
+  // Sort all the found fit objects in increasing quality 
   fFits.Sort();
   if (fDebug > 1) fFits.Print("s");
+
+  // Get the top-most fit
   ELossFit_t* ret = static_cast<ELossFit_t*>(fFits.At(i-1));
   if (!ret) {
-    AliWarningF("No fit found for %s, bin %3d", GetName(), b);
+    AliWarningF("No fit found for %s", GetName());
     return 0;
   }
   if (ret && fDebug > 0) {
-    if (fDebug > 1)
-      printf(" %d: ", i-1);
+    if (fDebug > 1) printf(" %d: ", i-1);
     ret->Print("s");
   }
   // We have to make a copy here, because other wise the clones array
@@ -1547,26 +1417,66 @@ AliFMDEnergyFitter::RingHistos::CalculateResiduals(EResidualMethod mode,
   }  
 }
 
-//____________________________________________________________________
+//__________________________________________________________________
 void
-AliFMDEnergyFitter::RingHistos::CreateOutputObjects(TList* dir)
+AliFMDEnergyFitter::RingHistos::FindBestFits(const TList*        d, 
+                                            AliFMDCorrELossFit& obj,
+                                            const TAxis&        eta)
 {
   // 
-  // Define outputs
+  // Find the best fits 
   // 
   // Parameters:
-  //    dir 
+  //    d              Parent list
+  //    obj            Object to add fits to
+  //    eta            Eta axis 
+  //    relErrorCut    Cut applied to relative error of parameter. 
+  //                       Note, for multi-particle weights, the cut 
+  //                       is loosend by a factor of 2 
+  //    chi2nuCut      Cut on @f$ \chi^2/\nu@f$ - 
+  //                       the reduced @f$\chi^2@f$ 
+  //    minWeightCut   Least valid @f$ a_i@f$ 
   //
-  DGUARD(fDebug, 2, "Define output in AliFMDEnergyFitter::RingHistos");
-  fList = DefineOutputList(dir);
 
-  fEtaEDists = new TList;
-  fEtaEDists->SetOwner();
-  fEtaEDists->SetName("EDists");
+  // Get our ring list from the output container 
+  TList* l = GetOutputList(d);
+  if (!l) return; 
 
-  fList->Add(fEtaEDists);
+  // Get the energy distributions from the output container 
+  TList* dists = static_cast<TList*>(l->FindObject("elossDists"));
+  if (!dists) { 
+    AliWarningF("Didn't find elossDists in %s", l->GetName());
+    l->ls();
+    return;
+  }
+  Int_t nBin = eta.GetNbins();
+  if (fBest.GetEntriesFast() <= 0) { 
+    AliWarningF("No fits found for %s", GetName());
+    return;
+  }
+  
+  for (Int_t b = 1; b <= nBin; b++) { 
+    ELossFit_t* best = static_cast<ELossFit_t*>(fBest.At(b));
+    if (!best) { 
+      // AliErrorF("No best fit found @ %d for %s", b, GetName());
+      continue;
+    }
+    // FindBestFit(b, dist, relErrorCut, chi2nuCut, minWeightCut);
+    best->fDet  = fDet; 
+    best->fRing = fRing;
+    best->fBin  = b; // 
+    if (fDebug > 0) {
+      printf("Bin # %3d: ", b); 
+      best->Print("s");
+    }
+    // Double_t eta = fAxis->GetBinCenter(b);
+    obj.SetFit(fDet, fRing, b, best); 
+    // new ELossFit_t(*best));
+  }
 }
 
+
+
 //____________________________________________________________________
 //
 // EOF
index 35859d4..99982ba 100644 (file)
  * @ingroup pwglf_forward_eloss
  */
 #include <TNamed.h>
-#include <TH1D.h>
+// #include <TH1D.h>
 #include <TAxis.h>
 #include <TList.h>
 #include <TObjArray.h>
 #include <TClonesArray.h>
 #include "AliFMDCorrELossFit.h"
 #include "AliForwardUtil.h"
+class TH1;
+class TH2;
 class AliESDFMD;
 class TFitResult;
 class TF1;
@@ -118,29 +120,12 @@ public:
    * @param title Title of object  - not significant 
    */
   AliFMDEnergyFitter(const char* title);
-  /** 
-   * Copy constructor 
-   * 
-   * @param o Object to copy from 
-   */
-  AliFMDEnergyFitter(const AliFMDEnergyFitter& o);
-  /** 
-   * Assignment operator 
-   * 
-   * @param o Object to assign from 
-   * 
-   * @return Reference to this 
-   */
-  AliFMDEnergyFitter& operator=(const AliFMDEnergyFitter& o);
 
+  // -----------------------------------------------------------------
   /** 
-   * Initialise the task
-   * 
-   * @param etaAxis The eta axis to use.  Note, that if the eta axis
-   * has already been set (using SetEtaAxis), then this parameter will be 
-   * ignored
+   * @{ 
+   * @name Setters of options and parameters 
    */
- void SetupForData(const TAxis& etaAxis);
   /** 
    * Set the eta axis to use.  This will force the code to use this
    * eta axis definition - irrespective of whatever axis is passed to
@@ -208,7 +193,7 @@ public:
    * 
    * @param n Max number of particle to try to fit 
    */
-  void SetNParticles(UShort_t n) { fNParticles = (n < 1 ? 1 : (n > 5 ? 5 : n)); }
+  void SetNParticles(UShort_t n) { fNParticles = (n<1 ? 1 : (n>5 ? 5 : n)); }
   /** 
    * Set the minimum number of entries each histogram must have 
    * before we try to fit our response function to it
@@ -299,6 +284,36 @@ public:
   {
     fRegularizationCut = cut;
   }
+  void SetSkips(UShort_t skip) { fSkips = skip; }
+  /** 
+   * Set the debug level.  The higher the value the more output 
+   * 
+   * @param dbg Debug level 
+   */
+  void SetDebug(Int_t dbg=1);
+  /* @} */
+  // -----------------------------------------------------------------
+  /** 
+   * @{ 
+   * @name Processing 
+   */
+  void Init();
+  /** 
+   * Define the output histograms.  These are put in a sub list of the
+   * passed list.   The histograms are merged before the parent task calls 
+   * AliAnalysisTaskSE::Terminate 
+   * 
+   * @param dir Directory to add to 
+   */
+  virtual void CreateOutputObjects(TList* dir);
+  /** 
+   * Initialise the task
+   * 
+   * @param etaAxis The eta axis to use.  Note, that if the eta axis
+   * has already been set (using SetEtaAxis), then this parameter will be 
+   * ignored
+   */
+ virtual void SetupForData(const TAxis& etaAxis);
   /** 
    * Fitter the input AliESDFMD object
    * 
@@ -308,36 +323,22 @@ public:
    * 
    * @return True on success, false otherwise 
    */
-  Bool_t Accumulate(const AliESDFMD& input, 
-                   Double_t         cent,
-                   Bool_t           empty);
+  virtual Bool_t Accumulate(const AliESDFMD& input, 
+                           Double_t         cent,
+                           Bool_t           empty);
   /** 
    * Scale the histograms to the total number of events 
    * 
    * @param dir Where the histograms are  
    */
-  void Fit(const TList* dir);
+  virtual void Fit(const TList* dir);
   /** 
    * Generate the corrections object 
    * 
    * @param dir List to analyse 
    */
   void MakeCorrectionsObject(TList* dir);
-  
-  /** 
-   * Define the output histograms.  These are put in a sub list of the
-   * passed list.   The histograms are merged before the parent task calls 
-   * AliAnalysisTaskSE::Terminate 
-   * 
-   * @param dir Directory to add to 
-   */
-  void CreateOutputObjects(TList* dir);
-  /** 
-   * Set the debug level.  The higher the value the more output 
-   * 
-   * @param dbg Debug level 
-   */
-  void SetDebug(Int_t dbg=1);
+  /** @} */
   /** 
    * Print information
    * 
@@ -352,9 +353,23 @@ public:
    * @return true if the parameter where read 
    */
   Bool_t ReadParameters(const TCollection* list);
-  void SetSkips(UShort_t skip) { fSkips = skip; }
 protected:
   /** 
+   * Copy constructor 
+   * 
+   * @param o Object to copy from 
+   */
+  AliFMDEnergyFitter(const AliFMDEnergyFitter& o);
+  /** 
+   * Assignment operator 
+   * 
+   * @param o Object to assign from 
+   * 
+   * @return Reference to this 
+   */
+  AliFMDEnergyFitter& operator=(const AliFMDEnergyFitter& o);
+
+  /** 
    * Internal data structure to keep track of the histograms
    */
   struct RingHistos : public AliForwardUtil::RingHistos
@@ -372,13 +387,13 @@ protected:
      */
     RingHistos(UShort_t d, Char_t r);
     /** 
-     * Copy constructor 
+     * Copy constructor - not defined
      * 
      * @param o Object to copy from 
      */
     RingHistos(const RingHistos& o);
     /** 
-     * Assignment operator 
+     * Assignment operator  - not defined
      * 
      * @param o Object to assign from 
      * 
@@ -390,11 +405,39 @@ protected:
      */
     ~RingHistos();
     /** 
+     * Make an axis with increasing bins 
+     * 
+     * @param n    Number of bins 
+     * @param min  Minimum 
+     * @param max  Maximum
+     * 
+     * @return An axis with quadratically increasing bin size 
+     */
+    virtual TArrayD MakeIncreasingAxis(Int_t    n, 
+                                      Double_t min, 
+                                      Double_t max) const;
+    /** 
+     * Make E/E_mip histogram 
+     * 
+     * @param name    Name of histogram
+     * @param title   Title of histogram
+     * @param eAxis   @f$\eta@f$ axis
+     * @param deMax   Maximum energy loss 
+     * @param nDeBins Number energy loss bins 
+     * @param incr    Whether to make bins of increasing size
+     */
+    TH2* Make(const char*  name, 
+             const char*  title, 
+             const TAxis& eAxis, 
+             Double_t     deMax=12, 
+             Int_t        nDeBins=300, 
+             Bool_t       incr=true);
+    /** 
      * Define outputs
      * 
      * @param dir 
      */
-    void CreateOutputObjects(TList* dir);
+    virtual void CreateOutputObjects(TList* dir);
     /** 
      * Initialise object 
      * 
@@ -404,11 +447,11 @@ protected:
      * @param nDEbins    Number of bins 
      * @param useIncrBin Whether to use an increasing bin size 
      */
-    void SetupForData(const TAxis& eAxis, 
-             const TAxis& cAxis,
-             Double_t     maxDE=10, 
-             Int_t        nDEbins=300, 
-             Bool_t       useIncrBin=true);
+    virtual void SetupForData(const TAxis& eAxis, 
+                             const TAxis& cAxis,
+                             Double_t     maxDE=10, 
+                             Int_t        nDEbins=300, 
+                             Bool_t       useIncrBin=true);
     /** 
      * Fill histogram 
      * 
@@ -417,12 +460,13 @@ protected:
      * @param icent  Centrality bin (1 based)
      * @param mult   Signal 
      */
-    void Fill(Bool_t empty, Int_t ieta, Int_t icent, Double_t mult);
+    virtual void Fill(Bool_t empty, Double_t eta, Int_t icent, Double_t mult);
     /** 
-     * Fit each histogram to up to @a nParticles particle responses.
+     * Get the the 2D histogram eloss name from our sub-list of @a dir
+     * and call the Fit function described below (with &fBest) as last
+     * argument.
      * 
      * @param dir         Output list 
-     * @param eta         Eta axis 
      * @param lowCut      Lower cut 
      * @param nParticles  Max number of convolved landaus to fit
      * @param minEntries  Minimum number of entries 
@@ -437,19 +481,70 @@ protected:
      * @param regCut      Regularization cut-off
      * @param residuals   Mode for residual plots
      *
-     * @return List of fits 
+     * @return List of fit parameters 
      */
-    TObjArray* Fit(TList*          dir, 
-                  const TAxis&    eta,
-                  Double_t        lowCut, 
-                  UShort_t        nParticles,
-                  UShort_t        minEntries,
-                  UShort_t        minusBins,
-                  Double_t        relErrorCut, 
-                  Double_t        chi2nuCut,
-                  Double_t        minWeight,
-                  Double_t        regCut,
-                  EResidualMethod residuals) const;
+    virtual TObjArray* Fit(TList*          dir, 
+                          Double_t        lowCut, 
+                          UShort_t        nParticles,
+                          UShort_t        minEntries,
+                          UShort_t        minusBins,
+                          Double_t        relErrorCut, 
+                          Double_t        chi2nuCut,
+                          Double_t        minWeight,
+                          Double_t        regCut,
+                          EResidualMethod residuals) const;
+    /** 
+     * Get the the 2D histogram @a name from our sub-list of @a
+     * dir. Then for each eta slice, try to fit the energu loss
+     * distribution up to @a nParticles particle responses.
+     *
+     * The fitted distributions (along with the functions fitted) are
+     * stored in a newly created sublist (<i>name</i>Dists).
+     *
+     * The fit parameters are also recorded in the newly created sub-list 
+     * <i>name</i>Results.  
+     *
+     * If @a residuals is not equal to kNoResiduals, then the
+     * residuals of the fits will be stored in the newly created
+     * sub-list <i>name</i>Residuals.
+     *
+     * A histogram named <i>name</i>Status is also generated and
+     * stored in the output list.
+     * 
+     * @param dir         Output list 
+     * @param name        Name of 2D base histogram in list
+     * @param lowCut      Lower cut 
+     * @param nParticles  Max number of convolved landaus to fit
+     * @param minEntries  Minimum number of entries 
+     * @param minusBins   Number of bins from peak to subtract to 
+     *                    get the fit range 
+     * @param relErrorCut Cut applied to relative error of parameter. 
+     *                    Note, for multi-particle weights, the cut 
+     *                    is loosend by a factor of 2 
+     * @param chi2nuCut   Cut on @f$ \chi^2/\nu@f$ - 
+     *                    the reduced @f$\chi^2@f$ 
+     * @param minWeight   Least weight ot consider
+     * @param regCut      Regularization cut-off
+     * @param residuals   Mode for residual plots
+     * @param scaleToPeak If true, scale distribution to peak value
+     * @param best        Optional array to store fits in
+     *
+     * @return List of fit parameters 
+     */
+    virtual TObjArray* FitSlices(TList*          dir, 
+                                const char*     name,
+                                Double_t        lowCut, 
+                                UShort_t        nParticles,
+                                UShort_t        minEntries,
+                                UShort_t        minusBins,
+                                Double_t        relErrorCut, 
+                                Double_t        chi2nuCut,
+                                Double_t        minWeight,
+                                Double_t        regCut,
+                                EResidualMethod residuals,
+                                Bool_t          scaleToPeak=true,
+                                TObjArray*      best=0) const;
+     
     /** 
      * Fit a signal histogram.  First, the bin @f$ b_{min}@f$ with
      * maximum bin content in the range @f$ [E_{min},\infty]@f$ is
@@ -473,55 +568,26 @@ protected:
      *                    the reduced @f$\chi^2@f$ 
      * @param minWeight   Least weight ot consider
      * @param regCut      Regularization cut-off
+     * @param scaleToPeak If true, scale distribution to peak value
+     * @param status      On return, contain the status code (0: OK, 1:
+     *                    empty, 2: low statistics, 3: fit failed)
      * 
      * @return The best fit function 
      */
-    ELossFit_t* FitHist(TH1*     dist,
-                       UShort_t bin, 
-                       Double_t lowCut, 
-                       UShort_t nParticles,
-                       UShort_t minusBins,
-                       Double_t relErrorCut, 
-                       Double_t chi2nuCut,
-                       Double_t minWeight,
-                       Double_t regCut) const;
-    /** 
-     * Calculate residuals of the fit 
-     * 
-     * @param mode   How to calculate 
-     * @param lowCut Lower cut 
-     * @param dist   Distribution 
-     * @param fit    Function fitted to distribution
-     * @param out    Output list to store residual histogram in
-     */
-    void CalculateResiduals(EResidualMethod   mode,
-                           Double_t         lowCut,
-                           TH1*             dist, 
-                           ELossFit_t*      fit, 
-                           TCollection*     out) const;
-    /** 
-     * Find the best fits 
-     * 
-     * @param d              Parent list
-     * @param obj            Object to add fits to
-     * @param eta            Eta axis 
-     * @param relErrorCut    Cut applied to relative error of parameter. 
-     *                       Note, for multi-particle weights, the cut 
-     *                       is loosend by a factor of 2 
-     * @param chi2nuCut      Cut on @f$ \chi^2/\nu@f$ - 
-     *                       the reduced @f$\chi^2@f$ 
-     * @param minWeightCut   Least valid @f$ a_i@f$ 
-     */
-    void FindBestFits(const TList*        d, 
-                     AliFMDCorrELossFit& obj,
-                     const TAxis&        eta,     
-                     Double_t            relErrorCut, 
-                     Double_t            chi2nuCut,
-                     Double_t            minWeightCut);
+    virtual ELossFit_t* FitHist(TH1*      dist,
+                               Double_t  lowCut, 
+                               UShort_t  nParticles,
+                               UShort_t  minEntries,
+                               UShort_t  minusBins,
+                               Double_t  relErrorCut, 
+                               Double_t  chi2nuCut,
+                               Double_t  minWeight,
+                               Double_t  regCut,
+                               Bool_t    scaleToPeak,
+                               UShort_t& status) const;
     /** 
      * Find the best fit 
      * 
-     * @param b              Bin number
      * @param dist           Histogram 
      * @param relErrorCut    Cut applied to relative error of parameter. 
      *                       Note, for multi-particle weights, the cut 
@@ -532,55 +598,39 @@ protected:
      * 
      * @return Best fit 
      */
-    ELossFit_t* FindBestFit(UShort_t b, 
-                           const TH1* dist,
-                           Double_t relErrorCut, 
-                           Double_t chi2nuCut,
-                           Double_t minWeightCut) const;
-#if 0
+    virtual ELossFit_t* FindBestFit(const TH1* dist,
+                                   Double_t   relErrorCut, 
+                                   Double_t   chi2nuCut,
+                                   Double_t   minWeightCut) const;
     /** 
-     * Check the result of the fit. Returns true if 
-     * - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
-     * - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters.  Note, 
-     *   for multi-particle fits, this requirement is relaxed by a 
-     *   factor of 2
-     * - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$ 
-     *   particle response 
-     * 
-     * @param r           Result to check
-     * @param relErrorCut Cut @f$ \delta_e@f$ applied to relative error 
-     *                    of parameter.  
-     * @param chi2nuCut   Cut @f$ \max{\chi^2/\nu}@f$ 
-     * 
-     * @return true if fit is good. 
-     */
-    Bool_t CheckResult(TFitResult* r,
-                      Double_t    relErrorCut, 
-                      Double_t    chi2nuCut,
-                      Double_t    minWeight) const;
-#endif
-    /** 
-     * Make an axis with increasing bins 
-     * 
-     * @param n    Number of bins 
-     * @param min  Minimum 
-     * @param max  Maximum
+     * Calculate residuals of the fit 
      * 
-     * @return An axis with quadratically increasing bin size 
+     * @param mode   How to calculate 
+     * @param lowCut Lower cut 
+     * @param dist   Distribution 
+     * @param fit    Function fitted to distribution
+     * @param out    Output list to store residual histogram in
      */
-    TArrayD MakeIncreasingAxis(Int_t n, Double_t min, Double_t max) const;
+    virtual void CalculateResiduals(EResidualMethod  mode,
+                                   Double_t         lowCut,
+                                   TH1*             dist, 
+                                   ELossFit_t*      fit, 
+                                   TCollection*     out) const;
     /** 
-     * Make E/E_mip histogram 
+     * Find the best fits.  This assumes that the array fBest has been
+     * filled with the best possible fits for each eta bin, and that
+     * the fits are placed according to the bin number of the eta bin.
+     *
+     * This is called by the parent class when generating the corretion 
+     * object. 
      * 
-     * @param ieta    Eta bin
-     * @param eMin    Least signal
-     * @param eMax    Largest signal 
-     * @param deMax   Maximum energy loss 
-     * @param nDeBins Number energy loss bins 
-     * @param incr    Whether to make bins of increasing size
+     * @param d    Parent list
+     * @param obj  Object to add fits to
+     * @param eta  Eta axis 
      */
-    void Make(Int_t ieta, Double_t eMin, Double_t eMax, 
-             Double_t deMax=12, Int_t nDeBins=300, Bool_t incr=true);
+    virtual void FindBestFits(const TList*        d, 
+                             AliFMDCorrELossFit& obj,
+                             const TAxis&        eta);
     /** 
      * Make a parameter histogram
      * 
@@ -590,9 +640,10 @@ protected:
      * 
      * @return 
      */
-    TH1D* MakePar(const char* name, const char* title, const TAxis& eta) const;
+    TH1* MakePar(const char* name, const char* title, const TAxis& eta) const;
     /** 
-     * Make a histogram that contains the results of the fit over the full ring 
+     * Make a histogram that contains the results of the fit over the
+     * full ring
      * 
      * @param name  Name 
      * @param title Title
@@ -604,22 +655,24 @@ protected:
      * 
      * @return The newly allocated histogram 
      */
-    TH1D* MakeTotal(const char* name, 
-                   const char* title, 
+    TH1* MakeTotal(const char*  name, 
+                   const char*  title, 
                    const TAxis& eta, 
-                   Int_t low, 
-                   Int_t high, 
-                   Double_t val, 
-                   Double_t err) const;
-    TH1D*        fEDist;        // Ring energy distribution 
-    TH1D*        fEmpty;        // Ring energy distribution for empty events
-    TList*       fEtaEDists;    // Energy distributions per eta bin. 
-    TList*       fList;
+                   Int_t        low, 
+                   Int_t        high, 
+                   Double_t     val, 
+                   Double_t     err) const;
+    TH1*                 fEDist;     // Ring energy distribution 
+    TH1*                 fEmpty;     // Ring energy dist for empty events
+    TH2*                 fHist;      // Two dimension Delta distribution
+    // TList*               fEtaEDists; // Energy distributions per eta bin. 
+    TList*               fList;
     mutable TObjArray    fBest;
     mutable TClonesArray fFits;
-    Int_t        fDebug;
-    ClassDef(RingHistos,3);
+    Int_t                fDebug;
+    ClassDef(RingHistos,4);
   };
+  virtual RingHistos* CreateRingHistos(UShort_t d, Char_t r) const;
   /** 
    * Get the ring histogram container 
    * 
@@ -653,27 +706,27 @@ protected:
    */
   static Bool_t CheckSkip(UShort_t d, Char_t r, UShort_t skips);
 
-  TList    fRingHistos;    // List of histogram containers
-  Double_t fLowCut;        // Low cut on energy
-  UShort_t fNParticles;    // Number of landaus to try to fit 
-  UShort_t fMinEntries;    // Minimum number of entries
-  UShort_t fFitRangeBinWidth;// Number of bins to subtract from found max
-  Bool_t   fDoFits;        // Whether to actually do the fits 
-  Bool_t   fDoMakeObject;  // Whether to make corrections object
-  TAxis    fEtaAxis;       // Eta axis 
-  TAxis    fCentralityAxis;// Centrality axis 
-  Double_t fMaxE;          // Maximum energy loss to consider 
-  Int_t    fNEbins;        // Number of energy loss bins 
-  Bool_t   fUseIncreasingBins; // Wheter to use increasing bin sizes 
-  Double_t fMaxRelParError;// Relative error cut
-  Double_t fMaxChi2PerNDF; // chi^2/nu cit
-  Double_t fMinWeight;     // Minimum weight value 
-  Int_t    fDebug;         // Debug level 
-  EResidualMethod fResidualMethod;// Whether to store residuals (debugging)
-  UShort_t        fSkips;         // Rings to skip when fitting 
+  TList           fRingHistos;        // List of histogram containers
+  Double_t        fLowCut;            // Low cut on energy
+  UShort_t        fNParticles;        // Number of landaus to try to fit 
+  UShort_t        fMinEntries;        // Minimum number of entries
+  UShort_t        fFitRangeBinWidth;  // N-bins to subtract from found max
+  Bool_t          fDoFits;            // Whether to actually do the fits 
+  Bool_t          fDoMakeObject;      // Whether to make corrections object
+  TAxis           fEtaAxis;           // Eta axis 
+  TAxis           fCentralityAxis;    // Centrality axis 
+  Double_t        fMaxE;              // Maximum energy loss to consider 
+  Int_t           fNEbins;            // Number of energy loss bins 
+  Bool_t          fUseIncreasingBins; // Wheter to use increasing bin sizes 
+  Double_t        fMaxRelParError;    // Relative error cut
+  Double_t        fMaxChi2PerNDF;     // chi^2/nu cit
+  Double_t        fMinWeight;         // Minimum weight value 
+  Int_t           fDebug;             // Debug level 
+  EResidualMethod fResidualMethod;    // Whether to store residuals (debugging)
+  UShort_t        fSkips;             // Rings to skip when fitting 
   Double_t        fRegularizationCut; // When to regularize the chi^2
 
-  ClassDef(AliFMDEnergyFitter,5); //
+  ClassDef(AliFMDEnergyFitter,6); //
 };
 
 #endif
index 89ac621..263f396 100644 (file)
@@ -91,6 +91,14 @@ AliFMDEnergyFitterTask::DefaultVertexAxis() const
 
 //____________________________________________________________________
 Bool_t
+AliFMDEnergyFitterTask::Setup()
+{
+  fEnergyFitter.Init();
+  return true;
+}
+
+//____________________________________________________________________
+Bool_t
 AliFMDEnergyFitterTask::Book()
 {
   // 
index 86d0e0e..a04aef5 100644 (file)
@@ -58,6 +58,12 @@ public:
    * @name Interface methods 
    */
   /** 
+   * Called on master when setting up the train. 
+   * 
+   * @return Always true 
+   */
+  virtual Bool_t Setup();
+  /** 
    * Book output objects. Derived class should define this to book
    * output objects on the processing output list @c fList before the
    * actual event processing.  This is called on the master and on
index d670609..6c7bc5b 100644 (file)
@@ -73,18 +73,23 @@ AliFMDEventInspector::AliFMDEventInspector()
     fVtxAxis(10,-10,10),
     fUseFirstPhysicsVertex(false),
     fUseV0AND(false),
-  fMinPileupContrib(3), 
-  fMinPileupDistance(0.8),
+    fMinPileupContrib(3), 
+    fMinPileupDistance(0.8),
     fUseDisplacedVertices(false),
-  fDisplacedVertex(),
-  fCollWords(),
-  fBgWords(),
-  fCentMethod("V0M"),
-  fMinCent(-1.0),
-  fMaxCent(-1.0),
-  fUsepA2012Vertex(false),
-  fRunNumber(0),
-  fMC(false)
+    fDisplacedVertex(),
+    fCollWords(),
+    fBgWords(),
+    fCentMethod("V0M"),
+    fMinCent(-1.0),
+    fMaxCent(-1.0),
+    fUsepA2012Vertex(false),
+    fRunNumber(0),
+    fMC(false),
+  fProdYear(-1),
+  fProdLetter('?'),
+  fProdPass(-1),
+  fProdSVN(-1),
+  fProdMC(false)
 {
   // 
   // Constructor 
@@ -130,7 +135,12 @@ AliFMDEventInspector::AliFMDEventInspector(const char* name)
   fMaxCent(-1.0),
   fUsepA2012Vertex(false),
   fRunNumber(0),
-  fMC(false)
+  fMC(false),
+  fProdYear(-1),
+  fProdLetter('?'),
+  fProdPass(-1),
+  fProdSVN(-1),
+  fProdMC(false)
 {
   // 
   // Constructor 
@@ -142,56 +152,6 @@ AliFMDEventInspector::AliFMDEventInspector(const char* name)
 }
 
 //____________________________________________________________________
-AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
-  : TNamed(o), 
-    fHEventsTr(o.fHEventsTr), 
-    fHEventsTrVtx(o.fHEventsTrVtx), 
-    fHEventsAccepted(o.fHEventsAccepted),
-    fHEventsAcceptedXY(o.fHEventsAcceptedXY),
-    fHTriggers(o.fHTriggers),
-    fHTriggerCorr(o.fHTriggerCorr),
-    fHType(o.fHType),
-    fHWords(o.fHWords),
-    fHCent(o.fHCent),
-    fHCentVsQual(o.fHCentVsQual),
-    fHStatus(o.fHStatus),
-    fHVtxStatus(o.fHVtxStatus),
-    fHTrgStatus(o.fHTrgStatus),
-    fLowFluxCut(o.fLowFluxCut),
-    fMaxVzErr(o.fMaxVzErr),
-    fList(o.fList),
-    fEnergy(o.fEnergy),
-    fField(o.fField), 
-    fCollisionSystem(o.fCollisionSystem),
-    fDebug(0),
-    fCentAxis(0),
-    fVtxAxis(o.fVtxAxis),
-    fUseFirstPhysicsVertex(o.fUseFirstPhysicsVertex),
-    fUseV0AND(o.fUseV0AND),
-    fMinPileupContrib(o.fMinPileupContrib), 
-    fMinPileupDistance(o.fMinPileupDistance),
-    fUseDisplacedVertices(o.fUseDisplacedVertices),
-    fDisplacedVertex(o.fDisplacedVertex),
-  fCollWords(),
-  fBgWords(),
-  fCentMethod(o.fCentMethod),
-  fMinCent(o.fMinCent),
-  fMaxCent(o.fMaxCent),
-  fUsepA2012Vertex(o.fUsepA2012Vertex),
-  fRunNumber(o.fRunNumber),
-  fMC(o.fMC)
-       
-{
-  // 
-  // Copy constructor 
-  // 
-  // Parameters:
-  //   o Object to copy from 
-  //
-  DGUARD(fDebug,1,"Copy CTOR of AliFMDEventInspector");
-}
-
-//____________________________________________________________________
 AliFMDEventInspector::~AliFMDEventInspector()
 {
   // 
@@ -200,74 +160,6 @@ AliFMDEventInspector::~AliFMDEventInspector()
   DGUARD(fDebug,1,"DTOR of AliFMDEventInspector");
   // if (fList)         delete fList;
 }
-//____________________________________________________________________
-AliFMDEventInspector&
-AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
-{
-  // 
-  // Assignement operator
-  // 
-  // Parameters:
-  //   o Object to assign from 
-  // 
-  // Return:
-  //    Reference to this object
-  //
-  DGUARD(fDebug,3,"Assignment of AliFMDEventInspector");
-  if (&o == this) return *this; 
-  TNamed::operator=(o);
-  fHEventsTr         = o.fHEventsTr;
-  fHEventsTrVtx      = o.fHEventsTrVtx;
-  fHEventsAccepted   = o.fHEventsAccepted;
-  fHEventsAcceptedXY = o.fHEventsAcceptedXY;
-  fHTriggers         = o.fHTriggers;
-  fHTriggerCorr      = o.fHTriggerCorr;
-  fHType             = o.fHType;
-  fHWords            = o.fHWords;
-  fHCent             = o.fHCent;
-  fHCentVsQual       = o.fHCentVsQual;
-  fHStatus           = o.fHStatus;
-  fHVtxStatus        = o.fHVtxStatus;
-  fHTrgStatus        = o.fHTrgStatus;
-  fLowFluxCut        = o.fLowFluxCut;
-  fMaxVzErr          = o.fMaxVzErr;
-  fDebug             = o.fDebug;
-  fList              = (o.fList ? new TList : 0);
-  fEnergy            = o.fEnergy;
-  fField             = o.fField;
-  fCollisionSystem   = o.fCollisionSystem;
-  fVtxAxis.Set(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(), 
-              o.fVtxAxis.GetXmax());
-  
-  fUseFirstPhysicsVertex = o.fUseFirstPhysicsVertex;
-  fUseV0AND              = o.fUseV0AND;
-  fMinPileupContrib      = o.fMinPileupContrib;
-  fMinPileupDistance     = o.fMinPileupDistance;
-  fUseDisplacedVertices  = o.fUseDisplacedVertices;
-  fDisplacedVertex       = o.fDisplacedVertex;
-  fCentMethod            = o.fCentMethod;
-  fMinCent              = o.fMinCent;
-  fMaxCent              = o.fMaxCent; 
-  fUsepA2012Vertex       = o.fUsepA2012Vertex;
-  fRunNumber             = o.fRunNumber;
-  fMC                    = o.fMC;
-
-  if (fList) { 
-    fList->SetName(GetName());
-    if (fHEventsTr)    fList->Add(fHEventsTr);
-    if (fHEventsTrVtx) fList->Add(fHEventsTrVtx);
-    if (fHTriggers)    fList->Add(fHTriggers);
-    if (fHTriggerCorr) fList->Add(fHTriggerCorr);
-    if (fHType)        fList->Add(fHType);
-    if (fHWords)       fList->Add(fHWords);
-    if (fHCent)        fList->Add(fHCent);
-    if (fHCentVsQual)  fList->Add(fHCentVsQual);
-    if (fHStatus)      fList->Add(fHStatus);
-    if (fHVtxStatus)   fList->Add(fHVtxStatus);
-    if (fHTrgStatus)   fList->Add(fHTrgStatus);
-  }
-  return *this;
-}
 
 //____________________________________________________________________
 void 
@@ -717,23 +609,23 @@ AliFMDEventInspector::StoreProduction()
 
   TString period            = p.GetLHCPeriod();
   // TString aliROOTVersion    = p.GetAlirootVersion();
-  Int_t   aliROOTSVN        = p.GetAlirootSvnVersion();
+  fProdSVN                  = p.GetAlirootSvnVersion();
   // TString rootVersion       = p.GetRootVersion();
   Int_t   rootSVN           = p.GetRootSvnVersion();
-  Int_t   pass              = p.GetRecoPass();
-  Bool_t  mc                = p.IsMC();
+  fProdPass                 = p.GetRecoPass();
+  fProdMC                   = p.IsMC();
 
   TObjArray* pp = TPRegexp("LHC([0-9]+)([a-z]+)").MatchS(period);
-  Int_t      yy = static_cast<TObjString*>(pp->At(1))->String().Atoi();
-  Char_t     ll = static_cast<TObjString*>(pp->At(2))->String()[0];
+  fProdYear     = static_cast<TObjString*>(pp->At(1))->String().Atoi();
+  fProdLetter   = static_cast<TObjString*>(pp->At(2))->String()[0];
   pp->Delete();
   
-  out->Add(AliForwardUtil::MakeParameter("year", yy));
-  out->Add(AliForwardUtil::MakeParameter("letter", Int_t(ll)));
-  out->Add(AliForwardUtil::MakeParameter("alirootSVN", aliROOTSVN));
+  out->Add(AliForwardUtil::MakeParameter("year", fProdYear));
+  out->Add(AliForwardUtil::MakeParameter("letter", Int_t(fProdLetter)));
+  out->Add(AliForwardUtil::MakeParameter("alirootSVN", fProdSVN));
   out->Add(AliForwardUtil::MakeParameter("rootSVN", rootSVN));
-  out->Add(AliForwardUtil::MakeParameter("pass", pass));
-  out->Add(AliForwardUtil::MakeParameter("mc", mc));
+  out->Add(AliForwardUtil::MakeParameter("pass", fProdPass));
+  out->Add(AliForwardUtil::MakeParameter("mc", fProdMC));
 }
 
   
@@ -1513,6 +1405,12 @@ AliFMDEventInspector::Print(Option_t*) const
   field.ReplaceAll("kG", " kG");
 
   gROOT->IncreaseDirLevel();
+  PFV("Production", "");
+  gROOT->IncreaseDirLevel();
+  PF("Period", "LHC%02d%c", fProdYear, fProdLetter);
+  PFB("MC anchor", fProdMC);
+  PFV("AliROOT Revision", fProdSVN);
+  gROOT->DecreaseDirLevel();
   PFV("Vertex bins", fVtxAxis.GetNbins());
   PF("Vertex Range", "[%+6.1f,%+6.1f]", fVtxAxis.GetXmin(), fVtxAxis.GetXmax());
   PFV("Low flux cut",          fLowFluxCut );
@@ -1531,7 +1429,8 @@ AliFMDEventInspector::Print(Option_t*) const
   if (!fCentAxis) {return; }
   Int_t nBin = fCentAxis->GetNbins();
   for (Int_t i = 0; i < nBin; i++) { 
-    if (i != 0 && (i % 10) == 0) std::cout << '\n' << "    ";
+    if (i != 0 && (i % 10) == 0) std::cout << '\n';
+    if ((i % 10) == 0)           std::cout << "    ";
     std::cout << std::setw(5) << fCentAxis->GetBinLowEdge(i+1) << '-';
   }
   std::cout << std::setw(5) << fCentAxis->GetBinUpEdge(nBin) << std::endl;
index 4996491..9d89019 100644 (file)
@@ -154,23 +154,9 @@ public:
    */
   AliFMDEventInspector(const char* name);
   /** 
-   * Copy constructor 
-   * 
-   * @param o Object to copy from 
-   */
-  AliFMDEventInspector(const AliFMDEventInspector& o);
-  /** 
    * Destructor 
    */
   virtual ~AliFMDEventInspector();
-  /** 
-   * Assignement operator
-   * 
-   * @param o Object to assign from 
-   * 
-   * @return Reference to this object
-   */
-  AliFMDEventInspector& operator=(const AliFMDEventInspector&);
 
   /** 
    * Initialize the object 
@@ -377,6 +363,39 @@ public:
    */
   ULong_t GetRunNumber() const { return fRunNumber; }
   /** 
+   * Get the production year.
+   * 
+   * - For real data, this is the year of the data taking 
+   * - For MC this is the year the production is anchored to.
+   * 
+   * @return A two-digit year (post millennium), or -1
+   */
+  Short_t GetProductionYear() const { return fProdYear; }
+  /** 
+   * Get the production period. 
+   * 
+   * - For real data, this is the period identifier of the data taking 
+   * - For MC data, this is the period identifier the production was
+   *   anchored to
+   * 
+   * @return Period identifier or null
+   */
+  Char_t  GetProductionPeriod() const { return fProdLetter; } 
+  /** 
+   * Get the AliROOT revision used for this production
+   * 
+   * @return SVN revision number or -1
+   */
+  Short_t GetProductionRevision() const { return fProdSVN; } 
+  /** 
+   * Check if the production was an MC production anchored in some
+   * real data.
+   * 
+   * @return true if this (MC) production was anchored 
+   */
+  Bool_t  IsProductionMC() const { return fProdMC; }
+  
+  /** 
    * Print information
    * 
    * @param option Not used 
@@ -415,6 +434,20 @@ public:
   static const char* CodeString(UInt_t mask);
 protected:
   /** 
+   * Copy constructor - not implemented
+   * 
+   * @param o Object to copy from 
+   */
+  AliFMDEventInspector(const AliFMDEventInspector& o);
+  /** 
+   * Assignement operator - not implemented
+   * 
+   * @param o Object to assign from 
+   * 
+   * @return Reference to this object
+   */
+  AliFMDEventInspector& operator=(const AliFMDEventInspector& o);
+  /** 
    * Cache the configure trigger classes from the physis selection.  
    * 
    * @param cache   where to cache the trigger class. 
@@ -582,45 +615,49 @@ protected:
   virtual Bool_t ReadCentrality(const AliESDEvent& esd, Double_t& cent,
                                UShort_t& qual) const;
 
-  TH1I*    fHEventsTr;    //! Histogram of events w/trigger
-  TH1I*    fHEventsTrVtx; //! Events w/trigger and vertex 
-  TH1I*    fHEventsAccepted; //! Events w/trigger and vertex in range 
-  TH2D*    fHEventsAcceptedXY; //! XY vtx with trigger and Z vertex in range 
-  TH1I*    fHTriggers;    //! Triggers
-  TH2I*    fHTriggerCorr; //! Correlation of triggers
-  TH1I*    fHType;        //! Type (low/high flux) of event
-  TH1I*    fHWords;       //! Trigger words 
-  TH1F*    fHCent;        //! Centrality 
-  TH2F*    fHCentVsQual;  //! Centrality vs quality 
-  TH1I*    fHStatus;      //! Event processing status 
-  TH1I*    fHVtxStatus;   //! Vertex processing status 
-  TH1I*    fHTrgStatus;   //! Trigger processing status 
-  Int_t    fLowFluxCut;   //  Low flux cut
-  Double_t fMaxVzErr;     //  Maximum error on v_z
-  TList*   fList;         //! Histogram container 
-  UShort_t fEnergy;       // CMS energy (per nucleon pair) [GeV]
-  Short_t  fField;        // L3 magnetic field [kG]
-  UShort_t fCollisionSystem; //  Collision system
-  Int_t    fDebug;        //  Debug level 
-  TAxis*   fCentAxis;     // Centrality axis used in histograms
-  TAxis    fVtxAxis;      //Vtx Axis 
-  Bool_t   fUseFirstPhysicsVertex; //Use the vtx code from p+p first physics
-  Bool_t   fUseV0AND;     //Use the vtx code from p+p first physics
-  UShort_t fMinPileupContrib; // Minimum number of contributors to 2nd
-                             // pile-up vertex
-  Double_t fMinPileupDistance; // Minimum distance of 2nd pile-up
-                              // vertex 
-  Bool_t   fUseDisplacedVertices; //Analyze displaced vertices?
+  TH1I*    fHEventsTr;            //! Histogram of events w/trigger
+  TH1I*    fHEventsTrVtx;         //! Events w/trigger and vertex 
+  TH1I*    fHEventsAccepted;      //! Events w/trigger and vertex in range 
+  TH2D*    fHEventsAcceptedXY;    //! XY vtx with trigger and Z vertex in range 
+  TH1I*    fHTriggers;            //! Triggers
+  TH2I*    fHTriggerCorr;         //! Correlation of triggers
+  TH1I*    fHType;                //! Type (low/high flux) of event
+  TH1I*    fHWords;               //! Trigger words 
+  TH1F*    fHCent;                //! Centrality 
+  TH2F*    fHCentVsQual;          //! Centrality vs quality 
+  TH1I*    fHStatus;              //! Event processing status 
+  TH1I*    fHVtxStatus;           //! Vertex processing status 
+  TH1I*    fHTrgStatus;           //! Trigger processing status 
+  Int_t    fLowFluxCut;           //  Low flux cut
+  Double_t fMaxVzErr;             //  Maximum error on v_z
+  TList*   fList;                 //! Histogram container 
+  UShort_t fEnergy;               // CMS energy (per nucleon pair) [GeV]
+  Short_t  fField;                // L3 magnetic field [kG]
+  UShort_t fCollisionSystem;      //  Collision system
+  Int_t    fDebug;                //  Debug level 
+  TAxis*   fCentAxis;             // Centrality axis used in histograms
+  TAxis    fVtxAxis;              // IP_z Axis 
+  Bool_t   fUseFirstPhysicsVertex;//Use the vtx code from p+p first physics
+  Bool_t   fUseV0AND;             // Use the vtx code from p+p first physics
+  UShort_t fMinPileupContrib;     // Min contributors to 2nd pile-up IP
+  Double_t fMinPileupDistance;    // Min distance of 2nd pile-up IP
+  Bool_t   fUseDisplacedVertices; // Analyze displaced vertices?
   AliDisplacedVertexSelection fDisplacedVertex; //Displaced vertex selector
-  TList    fCollWords;     //! Configured collision words 
-  TList    fBgWords;       //! Configured background words 
-  TString  fCentMethod;
-  Double_t fMinCent;  //min centrality
-  Double_t fMaxCent;  //max centrailty
-  Bool_t   fUsepA2012Vertex;// flag to use pA2012 Veretx selection
-  ULong_t  fRunNumber; // Current run number 
-  Bool_t   fMC;        // Is this MC input
-  ClassDef(AliFMDEventInspector,11); // Inspect the event 
+  TList    fCollWords;            //! Configured collision words 
+  TList    fBgWords;              //! Configured background words 
+  TString  fCentMethod;           // Centrality method
+  Double_t fMinCent;              // min centrality
+  Double_t fMaxCent;              // max centrailty
+  Bool_t   fUsepA2012Vertex;      // flag to use pA2012 Veretx selection
+  ULong_t  fRunNumber;            // Current run number 
+  Bool_t   fMC;                   // Is this MC input
+  Short_t  fProdYear;             // Production year 
+  Char_t   fProdLetter;           // Production letter 
+  Short_t  fProdPass;             // Pass number 
+  Int_t    fProdSVN;              // AliROOT revision used in production
+  Bool_t   fProdMC;               // True if anchor production
+
+  ClassDef(AliFMDEventInspector,12); // Inspect the event 
 };
 
 #endif
diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitter.cxx b/PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitter.cxx
new file mode 100644 (file)
index 0000000..7e1b111
--- /dev/null
@@ -0,0 +1,405 @@
+#include "AliFMDMCHitEnergyFitter.h"
+#include "AliMCEvent.h"
+#include "AliESDEvent.h"
+#include "AliStack.h"
+#include "AliFMDHit.h"
+#include "AliESDFMD.h"
+#include "AliMCAuxHandler.h"
+#include "AliGenEventHeader.h"
+#include "AliLog.h"
+#include "AliFMDEncodedEdx.h"
+#include <TH1.h>
+#include <TH2.h>
+#include <TNtuple.h>
+#include <TMath.h>
+
+
+//____________________________________________________________________
+AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter()
+  : AliFMDEnergyFitter(),
+    fSumPrimary(0),
+    fSumSecondary(0),
+    fIp(),
+    fNTrack(0),
+    fNPrimary(0),
+    fTuple(0)
+{
+}
+
+//____________________________________________________________________
+AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter(const char* title,
+                                                Bool_t useTuple)
+  : AliFMDEnergyFitter(title),
+    fSumPrimary(0),
+    fSumSecondary(0),
+    fIp(),
+    fNTrack(0),
+    fNPrimary(0),
+    fTuple(0)
+{
+  // Some defaults 
+  fUseIncreasingBins = true;
+  fDoFits            = true;
+  fDoMakeObject      = false;
+  fResidualMethod    = kNoResiduals;
+
+  if (!useTuple) return;
+
+  fTuple = new TNtuple("hits", "Hits", 
+                      "primary:eta:phi:edep:length:"
+                      "detector:ring:sector:strip:ipz");
+
+}
+//____________________________________________________________________
+AliFMDMCHitEnergyFitter::~AliFMDMCHitEnergyFitter()
+{}
+
+//____________________________________________________________________
+void
+AliFMDMCHitEnergyFitter::CreateOutputObjects(TList* dir)
+{
+  DGUARD(fDebug, 2, "Create output objects w/MC hits (%p)", dir);
+  AliFMDEnergyFitter::CreateOutputObjects(dir);
+  // TList* d = static_cast<TList*>(l->FindObject(GetName()));
+}
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitEnergyFitter::PreEvent(const AliMCEvent& mcInput)
+{
+  DGUARD(fDebug, 5, "Reset for event w/MC hits (%p)", &mcInput);
+  fSumPrimary.Reset(0);
+  fSumSecondary.Reset(0);
+
+  AliMCEvent&        mc        = const_cast<AliMCEvent&>(mcInput);
+  AliHeader*         header    = mc.Header();
+  AliStack*          stack     = mc.Stack();
+  AliGenEventHeader* genHeader = header->GenEventHeader();
+  
+  genHeader->PrimaryVertex(fIp);
+  fNTrack   = stack->GetNtrack();
+  fNPrimary = stack->GetNprimary();
+
+  return true;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitEnergyFitter::Event(const AliESDEvent& esdInput,
+                              const AliMCEvent&  mcInput,
+                              AliMCAuxHandler&   handler)
+{
+  DGUARD(fDebug, 5, "Process an event for MC hits (%p,%p,%p)",
+        &esdInput, &mcInput, &handler);
+  PreEvent(mcInput);
+
+  Int_t  nEntries = handler.GetNEntry();
+  DMSG(fDebug,5, "We got a total of %d particles", nEntries);
+
+  for (Int_t i = 0; i < nEntries; i++) {
+    TClonesArray* hits = handler.GetEntryArray(i);
+    if (!hits) continue;
+
+    AccumulateHits(mcInput, *hits);
+  }
+  return PostEvent(esdInput);
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDMCHitEnergyFitter::AccumulateHits(const AliMCEvent& mcInput,
+                                       const TClonesArray& hits)
+{
+  DGUARD(fDebug, 5, "Accumulate MC hit energy loss (%p,%p,%d)",
+        &mcInput, &hits, hits.GetEntries());
+
+  AliStack*  stack = const_cast<AliMCEvent&>(mcInput).Stack();
+  Float_t    cache[10];
+  Int_t      nHit  = hits.GetEntries();
+  Int_t      oldTr = -1;
+  UShort_t   oldD  = 0;
+  Char_t     oldR  = '\0';
+  for (Int_t j = 0; j < nHit; j++) { 
+       
+    // Check the hit 
+    AliFMDHit* hit  = static_cast<AliFMDHit*>(hits.At(j));
+    if (!hit) continue;
+
+    // Zero fill array 
+    if (fTuple) for (Int_t k = 0; k < 10; k++) cache[k] = 0;
+
+    // Get the track number 
+    Int_t iTr = hit->Track();
+    if (iTr < 0 || iTr >= fNTrack) 
+      AliWarningF("Track # %d seems out of bounds [0,%d]", 
+                 iTr, fNTrack-1);
+
+    // Get the track 
+    AliMCParticle* particle = 
+      static_cast<AliMCParticle*>(mcInput.GetTrack(iTr));
+       
+    // Particle type 0: unknown, 1: primary, 2: secondary
+    // Int_t flag = 0; - not used at the moment 
+
+    // Check if this charged and a primary 
+    if (particle) {
+      // Only charged tracks 
+      if (!particle->Charge() != 0) continue;
+    }
+    Bool_t   isPrimary = stack->IsPhysicalPrimary(iTr) && iTr < fNPrimary;
+    Double_t eloss     = hit->Edep();
+    Double_t length    = hit->Length();
+    // Double_t eta  = particle->Eta();
+    Double_t x         = hit->X() - fIp[0];
+    Double_t y         = hit->Y() - fIp[1];
+    Double_t z         = hit->Z() - fIp[2];
+    Double_t r         = TMath::Sqrt(x*x + y*y);
+    Double_t phi       = TMath::ATan2(y, x);
+    Double_t theta     = TMath::ATan2(r, z);
+    Double_t eta       = -TMath::Log(TMath::Tan(theta/2));
+    UShort_t det       = hit->Detector();
+    Char_t   rng       = hit->Ring();
+    UShort_t sec       = hit->Sector();
+    UShort_t str       = hit->Strip();
+
+    Bool_t isDistinct = (iTr != oldTr) || (det != oldD) || (rng != oldR);
+    if (isDistinct) {
+      RingHistos* h = static_cast<RingHistos*>(GetRingHistos(det, rng));
+      h->fKind->Fill(eta, isPrimary ? .5 : 1.5);
+    }
+    oldTr = iTr;
+    oldD  = det;
+    oldR  = rng;
+
+    AliFMDFloatMap& sum     =  (isPrimary ? fSumPrimary : fSumSecondary);
+    // Float_t         old     =  sum(det, rng, sec, str);
+    sum(det, rng, sec, str) += eloss;
+#if 0
+    Info("", "FMD%d%c[%02d,%03d]-%s: Adding %10.7f to %10.7f -> %10.7f",
+        det, rng, sec, str, (isPrimary ? "1st" : "2nd"), 
+        eloss, old, sum(det, rng, sec, str));
+#endif
+    if (!fTuple) continue;
+
+    // Leaves primary:eta:phi:edep:length:detector:ring:sector:strip:ipz
+    cache[0] = isPrimary ? 1 : 0;
+    cache[1] = eta;
+    cache[2] = phi;
+    cache[3] = eloss;
+    cache[4] = length;
+    cache[5] = hit->Detector();
+    cache[6] = Int_t(hit->Ring());
+    cache[7] = hit->Sector();
+    cache[8] = hit->Strip();
+    cache[9] = fIp[2];
+    fTuple->Fill(cache);
+  }
+  return true;
+}
+
+//____________________________________________________________________
+Bool_t AliFMDMCHitEnergyFitter::PostEvent(const AliESDEvent& input)
+{
+  DGUARD(fDebug, 5, "Convert MC hit energy loss (%p)", &input);
+
+  AliESDFMD*  esdFMD    = input.GetFMDData();
+  for (UShort_t d = 1; d <= 3; d++) { 
+    UShort_t nQ = (d == 1 ? 1 : 2);
+    for (UShort_t q = 0; q < nQ; q++) { 
+      UShort_t    nS = (q == 0 ?  20 :  40);
+      UShort_t    nT = (q == 0 ? 512 : 256);
+      Char_t      r  = (q == 0 ? 'I' : 'O');
+      RingHistos* h  = 
+       static_cast<AliFMDMCHitEnergyFitter::RingHistos*>(GetRingHistos(d,r));
+      if (!h) continue;
+
+      for (UShort_t s = 0; s < nS; s++) { 
+       for (UShort_t t = 0; t < nT; t++) {
+         Float_t  totPrim = fSumPrimary(d, r, s, t);
+         Float_t  totSec  = fSumSecondary(d, r, s, t);
+         Float_t  totAll  = totPrim + totSec;
+         Double_t esdEta  = esdFMD->Eta(d, r, s, t);
+         if (totAll  > 0) h->FillMC(0,  esdEta, totAll);
+         if (totPrim > 0) h->FillMC(1,  esdEta, totPrim);
+         if (totSec  > 0) h->FillMC(2,  esdEta, totSec);
+       }
+      }
+    }
+  }
+  return true;
+}
+
+//____________________________________________________________________
+AliFMDEnergyFitter::RingHistos*
+AliFMDMCHitEnergyFitter::CreateRingHistos(UShort_t d, Char_t r) const
+{ 
+  // DGUARD(fDebug, 1, "Create Ring cache of MC hit energy loss (FMD%d%c)",d,r);
+  return new AliFMDMCHitEnergyFitter::RingHistos(d,r);
+}
+
+//====================================================================
+AliFMDMCHitEnergyFitter::RingHistos::RingHistos()
+  : AliFMDEnergyFitter::RingHistos(), 
+    fPrimary(0), 
+    fSecondary(0),
+    fKind(0)
+{}
+
+//____________________________________________________________________
+AliFMDMCHitEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
+  : AliFMDEnergyFitter::RingHistos(d,r), 
+    fPrimary(0),
+    fSecondary(0),
+    fKind(0)
+{}
+//____________________________________________________________________
+TArrayD
+AliFMDMCHitEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t, 
+                                                       Double_t,
+                                                       Double_t) const
+{
+  // Make an increasing axis for ELoss distributions. 
+  // 
+  // We use the service function of AliFMDEncodedEdx to do this 
+  TArrayD ret;
+  const AliFMDEncodedEdx::Spec& s = AliFMDEncodedEdx::GetdEAxis();
+  s.FillBinArray(ret);
+
+  return ret;
+}
+//____________________________________________________________________
+void
+AliFMDMCHitEnergyFitter::RingHistos::SetupForData(const TAxis& eAxis, 
+                                                 const TAxis& /*cAxis*/, 
+                                                 Double_t maxDE, 
+                                                 Int_t    nDEbins, 
+                                                 Bool_t   useIncrBin)
+{
+  // AliFMDEnergyFitter::RingHistos::SetupForData(eAxis, cAxis, maxDE, nDEbins, 
+  // useIncrBin);
+
+  fHist      = Make("eloss", "#sum#Delta_{true} of all", 
+                   eAxis, maxDE, nDEbins, useIncrBin);
+  fPrimary   = Make("primary", "#sum#Delta_{true} of primaries", 
+                   eAxis, maxDE, nDEbins, useIncrBin);
+  fSecondary = Make("secondary","#sum#Delta_{true} of secondaries",
+                   eAxis, maxDE, nDEbins, useIncrBin);
+  fHist->SetXTitle("#Delta_{true}");
+  fSecondary->SetXTitle("#Delta_{true}");
+  fPrimary->SetXTitle("#Delta_{true}");
+
+  fKind       = 0;
+  if (eAxis.GetXbins()->GetArray()) 
+    fKind = new TH2I("kind", "Particle types", 
+                    eAxis.GetNbins(), eAxis.GetXbins()->GetArray(), 
+                    2, 0, 2);
+  else 
+    fKind = new TH2I("kind", "Particle types", 
+                    eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
+                    2, 0, 2);
+  fKind->SetXTitle("#eta");
+  fKind->GetYaxis()->SetBinLabel(1, "Primary");
+  fKind->GetYaxis()->SetBinLabel(2, "Secondary");
+  fKind->SetDirectory(0);
+  fKind->Sumw2();
+
+  fList->Add(fHist);
+  fList->Add(fPrimary);
+  fList->Add(fSecondary);
+  fList->Add(fKind);
+}
+
+//____________________________________________________________________
+void
+AliFMDMCHitEnergyFitter::RingHistos::FillMC(UShort_t flag, Double_t eta, 
+                                           Double_t mult)
+{
+  switch (flag) { 
+    // case 0: AliFMDEnergyFitter::RingHistos::Fill(false, eta, 0, mult); break;
+  case 0: fHist->Fill(eta, mult); break;
+  case 1: fPrimary->Fill(eta, mult); break;
+  case 2: fSecondary->Fill(eta, mult); break;
+  }
+}
+
+//____________________________________________________________________
+TObjArray*
+AliFMDMCHitEnergyFitter::RingHistos::Fit(TList*           dir, 
+                                        Double_t         lowCut, 
+                                        UShort_t         nParticles,
+                                        UShort_t         minEntries,
+                                        UShort_t         minusBins, 
+                                        Double_t         relErrorCut, 
+                                        Double_t         chi2nuCut,
+                                        Double_t         minWeight,
+                                        Double_t         regCut,
+                                        EResidualMethod  residuals) const
+{
+  TObjArray* all  = FitSlices(dir, "eloss", lowCut, nParticles, minEntries, 
+                             minusBins, relErrorCut, chi2nuCut, minWeight, 
+                             regCut, residuals, false, 0);
+  TObjArray* prim = FitSlices(dir, "primary", lowCut, nParticles, minEntries, 
+                             minusBins, relErrorCut, chi2nuCut, minWeight, 
+                             regCut, residuals, false, 0);
+  TObjArray* sec  = FitSlices(dir, "secondary", lowCut, nParticles, 
+                             minEntries, minusBins, relErrorCut, chi2nuCut, 
+                             minWeight, regCut, residuals, false, 0);
+  if (!all || !prim || !sec) {
+    AliWarningF("Failed to get results for all(%p), primary(%p), and/or "
+               "secondary(%p)", all, prim, sec);
+    return 0;
+  }
+  // Already added to the sub-folders
+  // dir->Add(all);
+  // dir->Add(prim);
+  // dir->Add(sec);
+
+  Int_t      nPrim = prim->GetEntriesFast();
+  TObjArray* out   = new TObjArray;
+  for (Int_t i = 0; i < nPrim; i++) { 
+    TH1* h = static_cast<TH1*>(prim->At(i));
+    if (!h) continue;
+
+    TAxis* xAxis  = h->GetXaxis();
+    TH2*   hh     = 0;
+    if (xAxis->GetXbins()->GetArray()) 
+      hh = new TH2D(h->GetName(), h->GetTitle(), 
+                   xAxis->GetNbins(), xAxis->GetXbins()->GetArray(),
+                   3, 0, 3);
+    else 
+      hh = new TH2D(h->GetName(), h->GetTitle(), 
+                   xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
+                   3, 0, 3);
+    hh->GetYaxis()->SetBinLabel(1, "Primary");
+    hh->GetYaxis()->SetBinLabel(2, "Secondary");
+    hh->GetYaxis()->SetBinLabel(3, "All");
+    hh->GetXaxis()->SetTitle("#eta");
+    out->Add(hh);
+  }
+  for (Int_t i = 0; i < nPrim; i++) { 
+    TH2* h = static_cast<TH2*>(out->At(i));
+    if (!h) continue;
+    
+    TH1* hp = static_cast<TH1*>(prim->At(i));
+    TH1* hs = static_cast<TH1*>(sec->At(i));
+    TH1* ha = static_cast<TH1*>(all->At(i));
+    TH1* hh[] = { hp, hs, ha };
+    for (Int_t j = 0; j < 3; j++) { 
+      TH1* ph = hh[j];
+      if (!ph) continue;
+
+      for (Int_t k = 1; k <= ph->GetNbinsX(); k++) { 
+       Double_t c = ph->GetBinContent(k);
+       Double_t e = ph->GetBinError(k);
+       h->SetBinContent(k, j+1, c);
+       h->SetBinError(k, j+1, e);
+      }
+    }
+  }
+  TList* l = GetOutputList(dir);
+  if (!l) return 0; 
+
+  out->SetName("compartive");
+  out->SetOwner();
+  l->Add(out);
+  return out;
+}
diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitter.h b/PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitter.h
new file mode 100644 (file)
index 0000000..e879cc0
--- /dev/null
@@ -0,0 +1,219 @@
+#ifndef ALIFMDMCHITENERGYFITTER_H
+#define ALIFMDMCHITENERGYFITTER_H
+#include "AliFMDEnergyFitter.h"
+#include "AliFMDFloatMap.h"
+#include "TArrayF.h"
+class AliMCAuxHandler;
+class AliMCEvent;
+class AliESDEvent;
+class TNtuple;
+class TClonesArray;
+
+class AliFMDMCHitEnergyFitter : public AliFMDEnergyFitter
+{
+public:
+  /** 
+   * Constructor - do not use
+   */
+  AliFMDMCHitEnergyFitter();
+  /** 
+   * Constructor 
+   * 
+   * @param title    Title of object - not significant 
+   * @param useTuple If true, also make an NTuple of hits on output 3
+   */
+  AliFMDMCHitEnergyFitter(const char* title, Bool_t useTuple=false);
+  /**
+   * Destructor
+   */
+  virtual ~AliFMDMCHitEnergyFitter();
+
+  /** 
+   * Define the output histograms.  These are put in a sub list of the
+   * passed list.   The histograms are merged before the parent task calls 
+   * AliAnalysisTaskSE::Terminate 
+   * 
+   * @param dir Directory to add to 
+   */
+  virtual void CreateOutputObjects(TList* dir);
+  /** 
+   * Set-up before processing an event 
+   * 
+   * @param mcInput MC input 
+   * 
+   * @return true
+   */
+  virtual Bool_t PreEvent(const AliMCEvent& mcInput);
+  /** 
+   * Process a single event 
+   * 
+   * @param esdInput ESD input 
+   * @param mcInput  MC input
+   * 
+   * @return true
+   */
+  virtual Bool_t Event(const AliESDEvent&  esdInput,
+                      const AliMCEvent&   mcInput,
+                      AliMCAuxHandler&    handler);
+  /** 
+   * Accumulate signals from MC hits
+   * 
+   * @param mcInput  MC input event
+   * @param hits     Clones array of hits 
+   * 
+   * @return 
+   */
+  virtual Bool_t AccumulateHits(const AliMCEvent&   mcInput,
+                               const TClonesArray& hits);
+  /** 
+   * Post-process accumulated signals 
+   * 
+   * @param esdInput ESD event 
+   * 
+   * @return true 
+   */
+  virtual Bool_t PostEvent(const AliESDEvent& esdInput);
+  
+  TNtuple* GetTuple() { return fTuple; }
+protected:
+  /**
+   * copy constructor - not implemented
+   */
+  AliFMDMCHitEnergyFitter(const AliFMDMCHitEnergyFitter&);
+  /**
+   * Assignment operator - not implemented
+   * 
+   */
+  AliFMDMCHitEnergyFitter& operator=(const AliFMDMCHitEnergyFitter&);
+
+  /**
+   * Container of ring histograms 
+   * 
+   */
+  struct RingHistos : public AliFMDEnergyFitter::RingHistos
+  {
+    /** 
+     * Default CTOR - do not use 
+     */
+    RingHistos();
+    /** 
+     * User CTOR 
+     * 
+     * @param d Detector number 
+     * @param r Ring identifier 
+     */
+    RingHistos(UShort_t d, Char_t r);
+    /** 
+     * DTOR
+     */
+    ~RingHistos() {}
+    /** 
+     * Copy constructor - not defined
+     * 
+     * @param o Object to copy from 
+     */
+    RingHistos(const RingHistos& o);
+    /** 
+     * Assignment operator  - not defined
+     * 
+     * @param o Object to assign from 
+     * 
+     * @return Reference to this 
+     */
+    RingHistos& operator=(const RingHistos& o);
+    /** 
+     * Create a bin array of increasing bins. This overload uses the
+     * service AliFMDEncodedEdx::Spec::FillBinArray. 
+     * 
+     * @param nBins Number of bins - ignored
+     * @param low   Low cut - ignored
+     * @param high  High cut - ignored
+     * 
+     * @return Array of bin boundaries 
+     */
+    TArrayD MakeIncreasingAxis(Int_t nBins, 
+                              Double_t low,
+                              Double_t high) const;
+    /** 
+     * Initialise object 
+     * 
+     * @param eAxis      Eta axis
+     * @param cAxis      Centrality axis 
+     * @param maxDE      Max energy loss to consider 
+     * @param nDEbins    Number of bins 
+     * @param useIncrBin Whether to use an increasing bin size 
+     */
+    virtual void SetupForData(const TAxis& eAxis, 
+                             const TAxis& cAxis,
+                             Double_t     maxDE=10, 
+                             Int_t        nDEbins=300, 
+                             Bool_t       useIncrBin=true);
+    /** 
+     * Fill in observation 
+     * 
+     * @param flag 0 - fill all, 1 - primary, 2 - secondary
+     * @param eta  Eta of particle observations
+     * @param mult Scaled energy loss 
+     */
+    virtual void FillMC(UShort_t flag, Double_t eta, Double_t mult);
+    /** 
+     * Fit the final distributions - called via Terminate 
+     * 
+     * @param dir           Containing directory
+     * @param lowCut        Lower cut on @f$\Delta/\Delta_{mip}@f$ 
+     * @param nParticles    Max. number of particle peaks to fit
+     * @param minEntries    Least number of entries required before fitting
+     * @param minusBins     Number of bins below the 1st peak we start fitting
+     * @param relErrorCut   Largest relative error on paramters
+     * @param chi2nuCut     Largest value of the @f$\chi^2/\nu@f$ 
+     * @param minWeight     Least weight to consider 
+     * @param regCut        When to regalurize 
+     * @param residuals     How to do residuals - if at all 
+     * 
+     * @return List of histograms of parameters 
+     */
+    TObjArray* Fit(TList*           dir, 
+                  Double_t         lowCut, 
+                  UShort_t         nParticles,
+                  UShort_t         minEntries,
+                  UShort_t         minusBins, 
+                  Double_t         relErrorCut, 
+                  Double_t         chi2nuCut,
+                  Double_t         minWeight,
+                  Double_t         regCut,
+                  EResidualMethod  residuals) const;
+    
+    TH2* fPrimary;   // @f$\Delta@f$ vs @f$\eta@f$ for primaries
+    TH2* fSecondary; // @f$\Delta@f$ vs @f$\eta@f$ for second.
+    TH2* fKind;      // Particle kind
+    ClassDef(RingHistos,1); // Cache of histograms per ring 
+  };
+  /** 
+   * Create a container of histograms for a single ring
+   * 
+   * @param d Detector 
+   * @param r Ring 
+   * 
+   * @return Newly allocated container 
+   */
+  AliFMDEnergyFitter::RingHistos* CreateRingHistos(UShort_t d, Char_t r) const;
+  /** Cache of per-strip energy loss of primaries */
+  AliFMDFloatMap   fSumPrimary;
+  /** Cache of per-strip energy loss of secondaries */
+  AliFMDFloatMap   fSumSecondary;
+  /** Cache of current MC IP */
+  TArrayF          fIp;
+  /** Cache of number of tracks */
+  Int_t            fNTrack;
+  /** Cache of numbr of primaries */
+  Int_t            fNPrimary;
+  /** Output nTuple of per-hit information */
+  TNtuple*         fTuple;
+
+  ClassDef(AliFMDMCHitEnergyFitter,1);
+};
+
+#endif
+// Local Variables:
+//  mode: C++
+// End:
diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitterTask.cxx b/PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitterTask.cxx
new file mode 100644 (file)
index 0000000..a10e847
--- /dev/null
@@ -0,0 +1,134 @@
+#include "AliFMDMCHitEnergyFitterTask.h"
+#include "AliForwardCorrectionManager.h"
+#include "AliAODForwardMult.h"
+#include "AliMCAuxHandler.h"
+#include <TROOT.h>
+#include <TNtuple.h>
+
+//____________________________________________________________________
+AliFMDMCHitEnergyFitterTask::AliFMDMCHitEnergyFitterTask(const char* name, 
+                                                        Bool_t      useTuple)
+  : AliBaseESDTask(name, "", &(AliForwardCorrectionManager::Instance())),
+    fEventInspector("event"), 
+    fEnergyFitter("fitter", useTuple), 
+    fHitHandler(0)
+{
+  fCloneList = true;
+  if (useTuple) DefineOutput(3, TNtuple::Class());
+}
+
+//____________________________________________________________________
+TAxis*
+AliFMDMCHitEnergyFitterTask::DefaultEtaAxis() const
+{
+  static TAxis* a = new TAxis(0, 0, 0);
+  return a;
+}
+//____________________________________________________________________
+TAxis*
+AliFMDMCHitEnergyFitterTask::DefaultVertexAxis() const
+{
+  static TAxis* a = new TAxis(10, -10, 10);
+  return a;
+}
+
+//____________________________________________________________________
+Bool_t 
+AliFMDMCHitEnergyFitterTask::Setup()
+{
+  DGUARD(fDebug,1,"Setting up the MC hit energy loss task");
+  fEnergyFitter.Init(); 
+  return true; 
+}
+
+//____________________________________________________________________
+Bool_t 
+AliFMDMCHitEnergyFitterTask::Book()
+{
+  DGUARD(fDebug,1,"Booking histograms for the MC hit energy loss task");
+  fNeededCorrections = 0;
+  fExtraCorrections  = 0;
+
+  fEnergyFitter.CreateOutputObjects(fList);  
+  // We have to add the handler here, since the subsiduary handlers of
+  // AliMCEventHandler are not streamed with the parent object - sigh!
+  // Hence, we add it at init time of the slaves. 
+  fHitHandler = AliMCAuxHandler::Create("FMD", "AliFMDHit");
+  return true;
+}
+
+//____________________________________________________________________
+Bool_t 
+AliFMDMCHitEnergyFitterTask::PreData(const TAxis& /*ipz*/, const TAxis& eta)
+{
+  DGUARD(fDebug,2,"Final setup of the MC hit energy loss task");
+  fEnergyFitter.SetupForData(eta);
+  if (fEnergyFitter.GetTuple()) 
+    PostData(3, fEnergyFitter.GetTuple());
+  return true;
+}
+//____________________________________________________________________
+Bool_t 
+AliFMDMCHitEnergyFitterTask::Event(AliESDEvent& esd)
+{
+  // --- Read in the data --------------------------------------------
+  LoadBranches();
+
+  Bool_t   lowFlux   = kFALSE;
+  UInt_t   triggers  = 0;
+  UShort_t ivz       = 0;
+  TVector3 ip;
+  Double_t cent      = 0;
+  UShort_t nClusters = 0;
+  UInt_t   found     = fEventInspector.Process(&esd, triggers, lowFlux, 
+                                              ivz, ip, cent, nClusters);
+  if (found & AliFMDEventInspector::kNoEvent)    return false;
+  if (found & AliFMDEventInspector::kNoTriggers) return false;
+  if (found & AliFMDEventInspector::kNoSPD)      return false;
+  if (found & AliFMDEventInspector::kNoFMD)      return false;
+  if (found & AliFMDEventInspector::kNoVertex)   return false;
+  if (found & AliFMDEventInspector::kBadVertex)  return false;
+
+  // do not process pile-up, A, C, and E events 
+  if (triggers & AliAODForwardMult::kPileUp)     return false;
+  if (triggers & AliAODForwardMult::kA)          return false;
+  if (triggers & AliAODForwardMult::kC)          return false;
+  if (triggers & AliAODForwardMult::kE)          return false;
+  
+  // We want only the events found by off-line 
+  if (!(triggers & AliAODForwardMult::kOffline)) return false;
+
+  // Perhaps we should also insist on MB only 
+  // if (fOnlyMB && (!(triggers & AliAODForwardMult::kInel))) return false;
+
+  AliMCEvent* mc = MCEvent();
+  if (!mc) return false;
+  
+  Bool_t ret = fEnergyFitter.Event(esd, *mc, *fHitHandler);
+
+  if (fEnergyFitter.GetTuple()) 
+    PostData(3, fEnergyFitter.GetTuple());
+
+  return ret;
+}
+#define PFB(N,FLAG)                            \
+  do {                                                                 \
+    AliForwardUtil::PrintName(N);                                      \
+    std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
+  } while(false)
+
+//____________________________________________________________________
+void   
+AliFMDMCHitEnergyFitterTask::Print(Option_t* option) const
+{
+  AliBaseESDTask::Print(option);
+  gROOT->IncreaseDirLevel();
+  fEnergyFitter.Print(option);
+  gROOT->DecreaseDirLevel();
+}
+
+//
+// EOF
+// 
+
+
diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitterTask.h b/PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitterTask.h
new file mode 100644 (file)
index 0000000..6a1f653
--- /dev/null
@@ -0,0 +1,191 @@
+#ifndef ALIFMDMCHITENERGYFITTERTASK_H
+#define ALIFMDMCHITENERGYFITTERTASK_H
+#include <AliBaseESDTask.h>
+#include <AliFMDMCHitEnergyFitter.h>
+#include <AliFMDEventInspector.h>
+class AliMCAuxHandler;
+
+/** 
+ * This task is designed to read in MC truth information about the
+ * energy loss in each strip, and then fit the distributions from
+ * secondaries and primaries separately.
+ * 
+ * Then (this is not implemented yet) it tries to deconvolve the
+ * contributions from secondaries and primaries separately from the
+ * total sum energy loss distribution, and in that way estimate the
+ * secondary to primary charge particle ratios per @f$\eta@f$ bin.  If
+ * the same procedure is applied to real data, then we can have an
+ * estimate - from data - of the ratio of secondaries to primaries and
+ * thus perhaps get a clearer picture of the secondary particle
+ * contamination.
+ *
+ * The the function fitted to the (scaled) energy loss
+ * (@f$\Delta/\Delta_{mip}@f$ distributions is 
+ *
+ * @f[
+ *    F(\Delta;\Delta_p,\xi,\sigma,\mathbf{a}) = 
+ *    \sum_{i=1}^{N}a_i f(\Delta;\Delta_i,\xi_i,\sigma_i)
+ * @f]
+ * where @f$ a@f$ is of length @f$ N@f$ and 
+ * @f[
+ *   f(\Delta;\Delta,\xi,\sigma) = 
+ *     \int_{-\infty}^{+\infty}d\Delta' L(\Delta,\Delta',\xi)
+ *     \frac{1}{\sqrt{2\pi\sigma^2}}
+ *     e^{-\frac{(\Delta'-\Delta_mp^2}{2\sigma^2}}
+ * @f]
+ * and 
+ * @f[ 
+ *     \Delta_i = i(\Delta_1+\xi\log(i))\\
+ *     \xi_i    = i\xi_1\\
+ *     \sigma_i = \sqrt{i}\sigma_1\\
+ *     a_1      = 1\\
+ *     \Delta_p = \Delta_1\\
+ *     \xi      = \xi_1\\
+ *     \sigma   = \sigma_1
+ * @f]
+ *
+ * See also AliForwardUtil::ELossFitter
+ */
+class AliFMDMCHitEnergyFitterTask : public AliBaseESDTask
+{
+public:
+  /** 
+   * Default CTOR - do not use
+   */
+  AliFMDMCHitEnergyFitterTask()
+    : AliBaseESDTask(),
+      fEventInspector(), 
+      fEnergyFitter(), 
+      fHitHandler(0)
+  {}
+  /** 
+   * CTOR
+   * 
+   * @param name     Name - not used
+   * @param useTuple Whether to use store an NTuple
+   */
+  AliFMDMCHitEnergyFitterTask(const char* name, 
+                             Bool_t      useTuple=false);
+  /** 
+   * DTOR
+   */
+  ~AliFMDMCHitEnergyFitterTask() {}
+
+  /** 
+   * Called when setting up the train on the client - i.e., called
+   * once before the job hits the worker nodes
+   * 
+   * @return true
+   */  
+  Bool_t Setup();
+  /** 
+   * Called at start-up on the clients - i.e., called once per worker
+   * node before the event processing.
+   * 
+   * @return true
+   */
+  Bool_t Book();
+  /** 
+   * Called after the first event was seen - i.e., called once per
+   * worker node just after the first event was seen, but before event
+   * processing.
+   * 
+   * @param ipz Interaction point Z--coordinate axis to use 
+   * @param eta @f$\eta=-\log[\tan^{-1}(\theta/2)]@f$ axis to use
+   * 
+   * @return true
+   */
+  Bool_t PreData(const TAxis& ipz, const TAxis& eta);
+  /** 
+   * Process a single event 
+   * 
+   * @param esd ESD input event 
+   * 
+   * @return true on success
+   */
+  Bool_t Event(AliESDEvent& esd);
+  /** 
+   * Finalize the task.  This is called once after all event
+   * processing and after all outputs have been merged.  This is
+   * called on the master/client once.
+   * 
+   * @return true on success
+   */
+  Bool_t Finalize() { fEnergyFitter.Fit(fResults); return true; }
+  /** 
+   * Print information to standard out
+   * 
+   * @param option Passed to sub objects as-is
+   */
+  void   Print(Option_t* option="") const;
+  /** 
+   * Get the event inspector 
+   * 
+   * @return Reference to the event inspector 
+   */
+  AliFMDEventInspector& GetEventInspector() { return fEventInspector; }
+  /** 
+   * Get the event inspector 
+   * 
+   * @return Constant reference to the event inspector 
+   */
+  const AliFMDEventInspector& GetEventInspector() const 
+  { 
+    return fEventInspector; 
+  }
+  /** 
+   * Get the energy fitter 
+   * 
+   * @return Reference to the energy fitter 
+   */
+  AliFMDMCHitEnergyFitter& GetEnergyFitter() { return fEnergyFitter; }
+  /** 
+   * Get the energy fitter 
+   * 
+   * @return Constant reference to the energy fitter 
+   */
+  const AliFMDMCHitEnergyFitter& GetEnergyFitter() const 
+  { 
+    return fEnergyFitter; 
+  }
+  void SetDebug(Int_t dbg) { 
+    AliBaseESDTask::SetDebug(dbg); GetEnergyFitter().SetDebug(dbg); }
+protected:
+  /** 
+   * Copy CTOR - not Implemented
+   * 
+   * @param o Object to copy from 
+   */  
+  AliFMDMCHitEnergyFitterTask(const AliFMDMCHitEnergyFitterTask& o);
+  /** 
+   * Assignment operator - not implemented
+   * 
+   * @param o Object to assign from 
+   * 
+   * @return Reference to this 
+   */
+  AliFMDMCHitEnergyFitterTask& operator=(const AliFMDMCHitEnergyFitterTask&o);
+  /** 
+   * Get the default @f$\eta=-\log[\tan^{-1}(\theta/2)]@f$ axis to use 
+   * 
+   * @return Pointer to static object 
+   */  
+  TAxis* DefaultEtaAxis() const;
+  /** 
+   * Get the default interaction point Z--coordinate axis to use 
+   * 
+   * @return Pointer to static object 
+   */  
+  TAxis* DefaultVertexAxis() const;
+  
+  AliFMDEventInspector    fEventInspector; // The event inspector
+  AliFMDMCHitEnergyFitter fEnergyFitter;   // The energy loss fitter 
+  AliMCAuxHandler*        fHitHandler;     // Handler for reading hits in
+
+  ClassDef(AliFMDMCHitEnergyFitterTask,1); // Task to fit Delta from MC hits
+};
+
+#endif
+// Local Variables:
+//  mode: C++
+// End:
index d207320..bba2552 100644 (file)
@@ -19,7 +19,7 @@
 #include "AliFMDMCDensityCalculator.h"
 #include "AliFMDMCCorrector.h"
 #include "AliFMDHistCollector.h"
-#include "AliFMDEnergyFitter.h"
+// #include "AliFMDEnergyFitter.h"
 #include "AliFMDEventPlaneFinder.h"
 #include <AliESDFMD.h>
 class AliESDEvent;
index eed8db4..327b550 100644 (file)
@@ -97,6 +97,9 @@ AliForwardMultiplicityBase::Book()
   GetHistCollector()   .CreateOutputObjects(fList);
   GetEventPlaneFinder()        .CreateOutputObjects(fList);
 
+  TAxis tmp(1, 0, 1);
+  fHistos.Init(tmp);
+
   if (fDebug > 1) fDoTiming = true;
   if (fDoTiming) { 
     fHTiming = new TProfile("timing", "Timing of task", 
@@ -135,9 +138,6 @@ AliForwardMultiplicityBase::CreateBranches(AliAODHandler* ah)
   TObject* obj   = &fAODFMD; ah->AddBranch("AliAODForwardMult", &obj);
   TObject* epobj = &fAODEP;  ah->AddBranch("AliAODForwardEP", &epobj);
 
-  TAxis tmp(1, 0, 1);
-  fHistos.Init(tmp);
-
   if (!fStorePerRing) return;
   
   AliWarning("Per-ring histograms in AOD\n"
index 75190d8..a3b92ac 100644 (file)
@@ -18,7 +18,7 @@
 #include "AliForwardUtil.h"
 #include "AliAODForwardMult.h"
 #include "AliAODForwardEP.h"
-class AliFMDEnergyFitter;
+// class AliFMDEnergyFitter;
 class AliFMDSharingFilter;
 class AliFMDDensityCalculator;
 class AliFMDCorrector;
index c7c0906..a52cca5 100644 (file)
@@ -20,7 +20,7 @@
 #include "AliFMDDensityCalculator.h"
 #include "AliFMDCorrector.h"
 #include "AliFMDHistCollector.h"
-#include "AliFMDEnergyFitter.h"
+// #include "AliFMDEnergyFitter.h"
 #include "AliFMDEventPlaneFinder.h"
 #include <AliESDFMD.h>
 class AliESDEvent;
index 0f03ec1..ee1aa45 100644 (file)
@@ -105,6 +105,14 @@ AliForwardQATask::DefaultVertexAxis() const
 
 //____________________________________________________________________
 Bool_t
+AliForwardQATask::Setup()
+{
+  fEnergyFitter.Init();
+  return true;
+}
+
+//____________________________________________________________________
+Bool_t
 AliForwardQATask::Book()
 {
   // 
index 55f7ed1..1731f46 100644 (file)
@@ -57,6 +57,12 @@ public:
    * @name Interface methods 
    */
   /** 
+   * Called when initialising the train. 
+   * 
+   * @return Always true
+   */
+  virtual Bool_t Setup();
+  /** 
    * Book output objects. Derived class should define this to book
    * output objects on the processing output list @c fList before the
    * actual event processing.  This is called on the master and on
index 28d2687..61fc1ea 100644 (file)
@@ -5,6 +5,7 @@
 #include <TFile.h>
 #include <TClonesArray.h>
 #include <TROOT.h>
+#include <TSystem.h>
 #include <AliStack.h>
 #include <AliMCEvent.h>
 
@@ -58,9 +59,13 @@ AliMCAuxHandler::Init(Option_t* opt)
   // 
   // @return true on success 
   AliDebugF(10,"AliMCAuxHandler::Init(\"%s\")", opt);
+  // Info("","AliMCAuxHandler::Init(\"%s\") - detector: %s, class: %s", 
+  //      opt, GetName(), GetTitle());
 
   TString option(opt);
-  if (option.EqualTo("proof") || option.EqualTo("local")) return true;
+  if (option.EqualTo("proof") || 
+      option.EqualTo("local") || 
+      option.EqualTo("lite")) return true;
 
   TString t = "Tree";
   TString b = "";
@@ -81,6 +86,14 @@ AliMCAuxHandler::Init(Option_t* opt)
     else 
       t = "";
   }
+  else {
+    AliWarningF("Couldn't get the class description of %s", GetTitle());
+    AliWarning("The list of loaded libraries is");
+    gSystem->ListLibraries("");
+    return false;
+  }
+      
+  // Info("Init", "Tree-name: %s, file-base: %s", t.Data(), b.Data());
   if (!t.IsNull()) fTreeName = t;
   if (!b.IsNull()) fFileBase = b;
 
@@ -113,6 +126,7 @@ AliMCAuxHandler::BeginEvent(Long64_t entry)
   // 
   // @return true on success
   AliDebugF(10,"AliMCAuxHandler::BeginEvent(%lld)", entry);
+  // Info("","AliMCAuxHandler::BeginEvent(%lld)", entry);
 
   if (entry == -1) 
     fEvent++;
@@ -151,6 +165,7 @@ AliMCAuxHandler::Notify(const char* path)
   //
   // @return true on success
   AliDebugF(10,"AliMCAuxHandler::Notify(\"%s\")", path);
+  // Info("","AliMCAuxHandler::Notify(\"%s\")", path);
   return true;
 }
 //____________________________________________________________________
@@ -160,7 +175,7 @@ AliMCAuxHandler::FinishEvent()
   // Called at the end of an event 
   // 
   // @return true on success
-  AliDebug(10,"AliMCAuxHandler::FinishEvent()");
+  // AliDebug(10,"AliMCAuxHandler::FinishEvent()");
   return true;
 }
 //____________________________________________________________________
@@ -195,6 +210,7 @@ AliMCAuxHandler::ResetIO()
 
   TString* path = GetParentPath();
   AliDebugF(10,"Got parent path %s", path ? path->Data() : "null");
+  // Info("ResetIO","Got parent path %s", path ? path->Data() : "null");
 
   if (fFile) { 
     delete fFile;
@@ -207,6 +223,7 @@ AliMCAuxHandler::OpenFile(Int_t fileNo)
 {
   TString* path = GetParentPath();
   AliDebugF(10,"Got parent path %s", path ? path->Data() : "null");
+  // Info("OpenFile","Got parent path %s", path ? path->Data() : "null");
   if (!path) return false;
 
   TString ext("");
@@ -218,7 +235,8 @@ AliMCAuxHandler::OpenFile(Int_t fileNo)
   TString fn = TString::Format("%s%s.%s%s.root", 
                               path->Data(), GetName(), 
                               fFileBase.Data(), ext.Data());
-  Info("Init", "Opening %s", fn.Data());
+  AliDebugF(10,"Opening %s", fn.Data());
+  // Info("OpenFile", "Opening %s", fn.Data());
   fFile = TFile::Open(fn, "READ");
   if (!fFile) { 
     AliErrorF("Failed to open %s", fn.Data());
@@ -238,6 +256,7 @@ AliMCAuxHandler::LoadEvent(Int_t iev)
   // 
   // @return true on success 
   AliDebugF(10,"AliMCAuxHandler::LoadEvent(%d)", iev);
+  // Info("","AliMCAuxHandler::LoadEvent(%d)", iev);
 
   Int_t iNew = iev / fNEventsPerFile;
   if (iNew != fFileNumber) { 
index 6187682..45f8b3f 100644 (file)
@@ -184,7 +184,14 @@ ForwardAODConfig(AliForwardMultiplicityBase* task)
   //   kDistance Only bins that are more than half the size of it neighbors
   task->GetHistCollector().SetFiducialMethod(AliFMDHistCollector::kByCut);
   // Additional diagnostics output - off by default
+  // 
+  // If this option is enabled, then the summed per-vertex, per-ring
+  // d2N/detadphi histograms will be stored in the output, as well as
+  // copies of the secondary maps
   // task->GetHistCollector().SetMakeBGHitMaps(true);
+  //
+  // If this option is enabled, then a 3D histogram will be made for
+  // each ring, summing dN/deta for each centrality bin.
   // task->GetHistCollector().SetMakeCentralitySums(true);
 
   // --- Eventplane Finder -------------------------------------------
index c6549fb..9f041f4 100644 (file)
@@ -4,7 +4,7 @@
  * @ingroup pwglf_forward_scripts
  */
 void
-LoadLibs()
+LoadLibs(bool alsoHit=false)
 {
 
   gROOT->LoadClass("TVirtualMC",              "libVMC");
@@ -18,6 +18,18 @@ LoadLibs()
   gROOT->LoadClass("AliAnalysisTaskSE",       "libANALYSISalice");
   gROOT->LoadClass("AliOADBPhysicsSelection"  "libOADB");
   gROOT->LoadClass("AliAODForwardMult",       "libPWGLFforward2");
+
+  if (!alsoHit) return;
+  
+  gROOT->LoadClass("TProof",                  "libProof");
+  gROOT->LoadClass("TGFrame",                 "libGui");
+  gROOT->LoadClass("AliCDBManager",           "libCDB");
+  gROOT->LoadClass("AliRawVEvent",            "libRAWDatabase");
+  gROOT->LoadClass("AliHit",                  "libSTEER");
+  gROOT->LoadClass("AliFMDDigit"              "libFMDbase");
+  gROOT->LoadClass("AliFMDHit",               "libFMDsim");
+  gROOT->LoadClass("AliFMDMCHitEnergyFitter", "libPWGLFforwardhit");
+  
 }
 //
 // EOF
index 3db5002..a104a3b 100644 (file)
@@ -75,7 +75,7 @@ protected:
     if (fOptions.Has("residuals")) resi = fOptions.Get("residuals"); 
 
     // --- Add the task ----------------------------------------------
-    AddTask("AddTaskFMDELoss.C", Form("%d,%d,%d,%d,\"%s\")", 
+    AddTask("AddTaskFMDELoss.C", Form("%d,%d,%d,%d,\"%s\"", 
                                      mc, cent, onlyMB, verb, resi.Data()));
   }
   /** 
diff --git a/PWGLF/FORWARD/analysis2/trains/MakeFMDMCHitTrain.C b/PWGLF/FORWARD/analysis2/trains/MakeFMDMCHitTrain.C
new file mode 100644 (file)
index 0000000..7a7e10a
--- /dev/null
@@ -0,0 +1,116 @@
+/**
+ * @file   MakeFMDMCHitTrain.C
+ * @author Christian Holm Christensen <cholm@nbi.dk>
+ * @date   Fri Jun  1 13:54:47 2012
+ * 
+ * @brief  
+ * 
+ * @ingroup pwglf_forward_trains_specific
+ */
+#include "TrainSetup.C"
+// #include "ParUtilities.C"
+
+//====================================================================
+/**
+ * Analysis train to make Forward and Central MC corrections
+ * 
+ *
+ * @ingroup pwglf_forward_mc
+ * @ingroup pwglf_forward_trains_specific
+ */
+class MakeFMDMCHitTrain : 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)
+   */
+  MakeFMDMCHitTrain(const  char* name) 
+    : TrainSetup(name)
+  {
+    fOptions.Add("use-tuple", "Whether to make an NTuple of hits");
+    fOptions.Set("type", "ESD");
+  }
+protected:
+  /** 
+   * Create the tasks 
+   * 
+   * @param mgr  Analysis manager 
+   */
+  void CreateTasks(AliAnalysisManager* mgr)
+  {
+    // --- Output file name ------------------------------------------
+    AliAnalysisManager::SetCommonFileName("forward_mchits.root");
+
+
+    // --- Load libraries/pars ---------------------------------------
+    fHelper->LoadLibrary("PWGLFforward2");
+    fHelper->LoadLibrary("Proof");
+    fHelper->LoadLibrary("Gui"); // Sigh! CDB depends on GUI!
+    fHelper->LoadLibrary("CDB");
+    fHelper->LoadLibrary("RAWDatabase");
+    fHelper->LoadLibrary("STEER");
+    fHelper->LoadLibrary("FMDbase");
+    fHelper->LoadLibrary("FMDsim");
+    fHelper->LoadLibrary("PWGLFforwardhit");
+    
+    // --- Set load path ---------------------------------------------
+    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2",
+                            gROOT->GetMacroPath()));
+
+    // --- Check if this is MC ---------------------------------------
+    if (!mgr->GetMCtruthEventHandler()) 
+      Fatal("CreateTasks", "No MC truth handler");
+    
+    TString args = TString::Format("%d,%d", 
+                                  fOptions.AsBool("use-tuple"),
+                                  fOptions.AsInt("verbose"));
+    if (!AddTask("AddTaskFMDMCHit.C", args))
+      Fatal("CreateTasks", "Couldn't add our task");
+  }
+  //__________________________________________________________________
+  /** 
+   * 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);
+  }
+  //__________________________________________________________________
+  /** 
+   * @return 0 - AOD disabled 
+   */
+  virtual AliVEventHandler* CreateOutputHandler(UShort_t) { return 0; }
+  /** 
+   * Do not the centrality selection
+   */
+  // void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {}
+  //__________________________________________________________________
+  const char* ClassName() const { return "MakeFMDMCHitTrain"; }
+};
+
+//
+// EOF
+//