]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/FORWARD/analysis2/AliForwardUtil.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardUtil.cxx
index b8ef2d162191f66475813c69f472eaef73c3431e..ac7176862be852557305db5c1342a3f59aa0c6d2 100644 (file)
@@ -3,14 +3,20 @@
 // 
 //
 #include "AliForwardUtil.h"
+//#include <ARVersion.h>
 #include <AliAnalysisManager.h>
 #include "AliAODForwardMult.h"
 #include <AliLog.h>
 #include <AliInputEventHandler.h>
+#include <AliAODInputHandler.h>
+#include <AliAODHandler.h>
+#include <AliAODEvent.h>
 #include <AliESDEvent.h>
+#include <AliAnalysisTaskSE.h>
 #include <AliPhysicsSelection.h>
 #include <AliTriggerAnalysis.h>
 #include <AliMultiplicity.h>
+#include <TParameter.h>
 #include <TH2D.h>
 #include <TH1I.h>
 #include <TF1.h>
 #include <TMath.h>
 #include <TError.h>
 #include <TROOT.h>
+#define FIT_OPTIONS "RNS"
+
+//====================================================================
+ULong_t AliForwardUtil::AliROOTRevision()
+{
+#ifdef ALIROOT_SVN_REVISION
+  return ALIROOT_SVN_REVISION;
+#elif defined(ALIROOT_REVISION)
+  static ULong_t ret = 0;
+  if (ret != 0) return ret;
+
+  // Select first 32bits of the 40byte long check-sum
+  TString rev(ALIROOT_REVISION, 8);
+  for (ULong_t i = 0; i < 8; i++) {
+    ULong_t p = 0;
+    switch (rev[i]) { 
+    case '0': p = 0; break;
+    case '1': p = 1; break;
+    case '2': p = 2; break;
+    case '3': p = 3; break;
+    case '4': p = 4; break;
+    case '5': p = 5; break;
+    case '6': p = 6; break;
+    case '7': p = 7; break;
+    case '8': p = 8; break;
+    case '9': p = 9; break;
+    case 'a': case 'A': p = 10; break;
+    case 'b': case 'B': p = 11; break;
+    case 'c': case 'C': p = 12; break;
+    case 'd': case 'D': p = 13; break;
+    case 'e': case 'E': p = 14; break;
+    case 'f': case 'F': p = 15; break;
+    }
+    ret |= (p << (32-4*(i+1)));
+  }
+  return ret;
+#else
+  return 0;
+#endif
+}
+//____________________________________________________________________
+ULong_t AliForwardUtil::AliROOTBranch()
+{
+  // Do something here when we switch to git - sigh!
+#if !defined(ALIROOT_SVN_BRANCH) && !defined(ALIROOT_BRANCH) 
+  return 0;
+#endif
+  static ULong_t ret = 0;
+  if (ret != 0) return ret;
+  TString str;
+  TString top;
+#ifdef ALIROOT_SVN_BRANCH
+  str = ALIROOT_SVN_BRANCH;
+  top = "trunk";
+#elif defined(ALIROOT_BRANCH)
+  str = ALIROOT_BRANCH;
+  top = "master";
+#endif
+  if (str.IsNull()) return 0xFFFFFFFF;
+  if (str[0] == 'v') str.Remove(0,1);
+  if (str.EqualTo(top)) return ret = 0xFFFFFFFF;
+
+  TObjArray*   tokens = str.Tokenize("-");
+  TObjString*  pMajor = tokens->GetEntries()>0 ? 
+    (static_cast<TObjString*>(tokens->At(0))) : 0;
+  TObjString*  pMinor = tokens->GetEntries()>1 ? 
+    (static_cast<TObjString*>(tokens->At(1))) : 0;
+  TObjString*  pRelea = tokens->GetEntries() > 2 ? 
+    static_cast<TObjString*>(tokens->At(2)) : 0;
+  TObjString* pAn     = tokens->GetEntries() > 3 ? 
+    static_cast<TObjString*>(tokens->At(3)) : 0;
+  TString sMajor,sMinor,sRelea;
+  if (pMajor) sMajor = pMajor->String().Strip(TString::kLeading, '0'); 
+  if (pMinor) sMinor = pMinor->String().Strip(TString::kLeading, '0');
+  if (pRelea) sRelea = pRelea->String().Strip(TString::kLeading, '0');
+  //
+  ret = (((sMajor.Atoi() & 0xFF) << 12) |
+    ((sMinor.Atoi() & 0xFF) <<  8) |
+    ((sRelea.Atoi() & 0xFF) <<  4) |
+    (pAn ? 0xAA : 0));
+  
+  return ret;
+}
 
 //====================================================================
 UShort_t
@@ -42,6 +131,7 @@ AliForwardUtil::ParseCollisionSystem(const char* sys)
   // we do pA first to avoid pp catch on ppb string (AH)
   if (s.Contains("p-pb")  || s.Contains("ppb"))   return AliForwardUtil::kPPb;
   if (s.Contains("p-a")   || s.Contains("pa"))    return AliForwardUtil::kPPb;
+  if (s.Contains("a-p")   || s.Contains("ap"))    return AliForwardUtil::kPPb;
   if (s.Contains("p-p")   || s.Contains("pp"))    return AliForwardUtil::kPP; 
   if (s.Contains("pb-pb") || s.Contains("pbpb"))  return AliForwardUtil::kPbPb;
   if (s.Contains("a-a")   || s.Contains("aa"))    return AliForwardUtil::kPbPb;
@@ -71,33 +161,90 @@ AliForwardUtil::CollisionSystemString(UShort_t sys)
   return "unknown";
 }
 //____________________________________________________________________
+Float_t
+AliForwardUtil::BeamRapidity(Float_t beam, UShort_t z, UShort_t a)
+{
+  const Double_t pMass = 9.38271999999999995e-01;
+  const Double_t nMass = 9.39564999999999984e-01;
+  Double_t       beamE = z * beam / 2;
+  Double_t       beamM = z * pMass + (a - z) * nMass;
+  Double_t       beamP = TMath::Sqrt(beamE * beamE - beamM * beamM);
+  Double_t       beamY = .5* TMath::Log((beamE+beamP) / (beamE-beamP));
+  return beamY;
+}
+//____________________________________________________________________
+Float_t
+AliForwardUtil::CenterOfMassEnergy(Float_t beam, 
+                                  UShort_t z1, 
+                                  UShort_t a1, 
+                                  Short_t z2, 
+                                  Short_t a2) 
+{
+  // Calculate the center of mass energy given target/projectile 
+  // mass and charge numbers
+  if (z2 < 0) z2 = z1;
+  if (a2 < 0) a2 = a1;
+  return TMath::Sqrt(Float_t(z1*z2)/a1/a2) * beam;
+}
+//____________________________________________________________________
+Float_t
+AliForwardUtil::CenterOfMassRapidity(UShort_t z1, 
+                                    UShort_t a1, 
+                                    Short_t z2, 
+                                    Short_t a2) 
+{
+  // Calculate the center of mass rapidity (shift) given target/projectile 
+  // mass and charge numbers
+  if (z2 < 0) z2 = z1;
+  if (a2 < 0) a2 = a1;
+  if (z2 == z1 && a2 == a1) return 0;
+  return .5 * TMath::Log(Float_t(z1*a2)/z2/a1);
+}
+
+namespace {
+  UShort_t CheckSNN(Float_t energy)
+  {
+    if (TMath::Abs(energy - 900.)   < 10)  return 900;
+    if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
+    if (TMath::Abs(energy - 2760.)  < 20)  return 2760;
+    if (TMath::Abs(energy - 4400.)  < 10)  return 4400;
+    if (TMath::Abs(energy - 5022.)  < 10)  return 5023;
+    if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
+    if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
+    if (TMath::Abs(energy - 8000.)  < 10)  return 8000;
+    if (TMath::Abs(energy - 10000.) < 10)  return 10000;
+    if (TMath::Abs(energy - 14000.) < 10)  return 14000;
+    return 0;
+  }
+}
+//____________________________________________________________________
 UShort_t
-AliForwardUtil::ParseCenterOfMassEnergy(UShort_t /* sys */, Float_t v)
+AliForwardUtil::ParseCenterOfMassEnergy(UShort_t sys, Float_t beam)
 {
   // 
   // Parse the center of mass energy given as a float and return known 
   // values as a unsigned integer
   // 
   // Parameters:
-  //    sys  Collision system (needed for AA)
-  //    cms  Center of mass energy * total charge 
+  //    sys   Collision system (needed for AA)
+  //    beam  Center of mass energy * total charge 
   // 
   // Return:
   //    Center of mass energy per nucleon
   //
-  Float_t energy = v
+  Float_t energy = beam
   // Below no longer needed apparently
   // if (sys == AliForwardUtil::kPbPb) energy = energy / 208 * 82;
-  if (TMath::Abs(energy - 900.)   < 10)  return 900;
-  if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
-  if (TMath::Abs(energy - 2750.)  < 20)  return 2750;
-  if (TMath::Abs(energy - 4400.)  < 10)  return 4400;
-  if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
-  if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
-  if (TMath::Abs(energy - 8000.)  < 10)  return 8000;
-  if (TMath::Abs(energy - 10000.) < 10)  return 10000;
-  if (TMath::Abs(energy - 14000.) < 10)  return 14000;
-  return 0;
+  if (sys == AliForwardUtil::kPPb) 
+    energy = CenterOfMassEnergy(beam, 82, 208, 1, 1);
+  else if (sys == AliForwardUtil::kPbPb) 
+    energy = CenterOfMassEnergy(beam, 82, 208, 82, 208);
+  UShort_t ret = CheckSNN(energy);
+  if (ret > 1) return ret;
+  if (sys == AliForwardUtil::kPbPb || sys == AliForwardUtil::kPPb) {
+    ret = CheckSNN(beam);
+  }
+  return ret;
 }
 //____________________________________________________________________
 const char* 
@@ -148,31 +295,261 @@ AliForwardUtil::MagneticFieldString(Short_t f)
   return Form("%01dkG", f);
 }
 //_____________________________________________________________________
-Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Double_t zvtx)
+AliAODEvent* AliForwardUtil::GetAODEvent(AliAnalysisTaskSE* task)
 {
-  //Calculate eta from strip with vertex (redundant with AliESDFMD::Eta but support displaced vertices)
+  // Check if AOD is the output event
+  if (!task) ::Fatal("GetAODEvent", "Null task given, cannot do that");
+
+  AliAODEvent* ret = task->AODEvent();
+  if (ret) return ret; 
   
-  //Get max R of ring
-  Double_t maxR = 0;
-  Double_t minR = 0;
-  Bool_t inner = false;
-  switch (ring) { 
-  case 'i': case 'I': maxR = 17.2; minR =  4.5213; inner = true;  break;
-  case 'o': case 'O': maxR = 28.0; minR = 15.4;    inner = false; break;
-  default: 
-    return -99999;
-  }
-
-  Double_t   rad       =  maxR- minR;
-  Double_t   nStrips   = (ring == 'I' ? 512 : 256);
-  Double_t   segment   = rad / nStrips;
-  Double_t   r         =  minR + segment*strip;
-  Int_t hybrid = sec / 2;
+  // Check if AOD is the input event 
+  ret = dynamic_cast<AliAODEvent*>(task->InputEvent());
+  if (!ret) ::Warning("GetAODEvent", "No AOD event found");
+  
+  return ret; 
+}
+//_____________________________________________________________________
+UShort_t AliForwardUtil::CheckForAOD()
+{
+  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
+  if (dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler())) {
+    // ::Info("CheckForAOD", "Found AOD Input handler");
+    return 1;
+  }
+  if (dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler())) {
+    // ::Info("CheckForAOD", "Found AOD Output handler");
+    return 2;
+  }
+
+  ::Warning("CheckForAOD", 
+           "Neither and input nor output AOD handler is specified");
+  return 0;
+}
+//_____________________________________________________________________
+Bool_t AliForwardUtil::CheckForTask(const char* clsOrName, Bool_t cls)
+{
+  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
+  if (!cls) { 
+    AliAnalysisTask* t = am->GetTask(clsOrName);
+    if (!t) { 
+      ::Warning("CheckForTask", "Task %s not found in manager", clsOrName);
+      return false;
+    }
+    ::Info("CheckForTask", "Found task %s", clsOrName);
+    return true;
+  }
+  TClass* dep = gROOT->GetClass(clsOrName);
+  if (!dep) { 
+    ::Warning("CheckForTask", "Unknown class %s for needed task", clsOrName);
+    return false;
+  }
+  TIter next(am->GetTasks());
+  TObject* o = 0;
+  while ((o = next())) { 
+    if (o->IsA()->InheritsFrom(dep)) {
+      ::Info("CheckForTask", "Found task of class %s: %s", 
+            clsOrName, o->GetName());
+      return true;
+    }
+  }
+  ::Warning("CheckForTask", "No task of class %s was found", clsOrName);
+  return false;
+}
+
+//_____________________________________________________________________
+TObject* AliForwardUtil::MakeParameter(const Char_t* name, UShort_t value)
+{
+  TParameter<int>* ret = new TParameter<int>(name, value);
+  ret->SetMergeMode('f');
+  ret->SetUniqueID(value);
+  return ret;
+}
+//_____________________________________________________________________
+TObject* AliForwardUtil::MakeParameter(const Char_t* name, Int_t value)
+{
+  TParameter<int>* ret = new TParameter<int>(name, value);
+  ret->SetMergeMode('f');
+  ret->SetUniqueID(value);
+  return ret;
+}
+//_____________________________________________________________________
+TObject* AliForwardUtil::MakeParameter(const Char_t* name, ULong_t value)
+{
+  TParameter<Long_t>* ret = new TParameter<Long_t>(name, value);
+  ret->SetMergeMode('f');
+  ret->SetUniqueID(value);
+  return ret;
+}
+//_____________________________________________________________________
+TObject* AliForwardUtil::MakeParameter(const Char_t* name, Double_t value)
+{
+  TParameter<double>* ret = new TParameter<double>(name, value);
+  // Float_t v = value;
+  // UInt_t* tmp = reinterpret_cast<UInt_t*>(&v);
+  ret->SetMergeMode('f');
+  // ret->SetUniqueID(*tmp);
+  return ret;
+}
+//_____________________________________________________________________
+TObject* AliForwardUtil::MakeParameter(const Char_t* name, Bool_t value)
+{
+  TParameter<bool>* ret = new TParameter<bool>(name, value);
+  ret->SetMergeMode('f');
+  ret->SetUniqueID(value);
+  return ret;
+}
+
+//_____________________________________________________________________
+void AliForwardUtil::GetParameter(TObject* o, UShort_t& value)
+{
+  if (!o) return;
+  TParameter<int>* p = static_cast<TParameter<int>*>(o);
+  if (p->TestBit(BIT(19)))
+    value = p->GetVal(); 
+  else
+    value = o->GetUniqueID();
+}
+//_____________________________________________________________________
+void AliForwardUtil::GetParameter(TObject* o, Int_t& value)
+{
+  if (!o) return;
+  TParameter<int>* p = static_cast<TParameter<int>*>(o);
+  if (p->TestBit(BIT(19)))
+    value = p->GetVal(); 
+  else
+    value = o->GetUniqueID();
+}
+//_____________________________________________________________________
+void AliForwardUtil::GetParameter(TObject* o, ULong_t& value)
+{
+  if (!o) return;
+  TParameter<Long_t>* p = static_cast<TParameter<Long_t>*>(o);
+  if (p->TestBit(BIT(19)))
+    value = p->GetVal(); 
+  else
+    value = o->GetUniqueID();
+}
+//_____________________________________________________________________
+void AliForwardUtil::GetParameter(TObject* o, Double_t& value)
+{
+  if (!o) return;
+  TParameter<double>* p = static_cast<TParameter<double>*>(o);
+  if (p->TestBit(BIT(19)))
+    value = p->GetVal(); // o->GetUniqueID();
+  else {
+    UInt_t  i = o->GetUniqueID();
+    Float_t v = *reinterpret_cast<Float_t*>(&i);
+    value = v;
+  }
+}
+//_____________________________________________________________________
+void AliForwardUtil::GetParameter(TObject* o, Bool_t& value)
+{
+  if (!o) return;
+  TParameter<bool>* p = static_cast<TParameter<bool>*>(o);
+  if (p->TestBit(BIT(19)))
+    value = p->GetVal(); // o->GetUniqueID();
+  else
+    value = o->GetUniqueID();
+}
+  
+#if 0
+//_____________________________________________________________________
+Double_t AliForwardUtil::GetStripR(Char_t ring, UShort_t strip)
+{
+  // Get max R of ring
+  // 
+  // Optimized version that has a cache 
+  static TArrayD inner;
+  static TArrayD outer; 
+  if (inner.GetSize() <= 0 || outer.GetSize() <= 0) {
+    const Double_t minR[] = {  4.5213, 15.4 };
+    const Double_t maxR[] = { 17.2,    28.0 };
+    const Int_t    nStr[] = { 512,     256  };
+    for (Int_t q = 0; q < 2; q++) { 
+      TArrayD& a = (q == 0 ? inner : outer);
+      a.Set(nStr[q]);
+
+      for (Int_t it = 0; it < nStr[q]; it++) {
+       Double_t   rad     = maxR[q] - minR[q];
+       Double_t   segment = rad / nStr[q];
+       Double_t   r       = minR[q] + segment*strip;
+       a[it]              = r;
+      }
+    }
+  }
+  if (ring == 'I' || ring == 'i') return inner.At(strip);
+  return outer.At(strip);
+}
+#else
+//_____________________________________________________________________
+Double_t AliForwardUtil::GetStripR(Char_t ring, UShort_t strip)
+{
+  // Get max R of ring
+  // 
+  // New implementation has only one branch
+  const Double_t minR[] = {  4.5213, 15.4 };
+  const Double_t maxR[] = { 17.2,    28.0 };
+  const Int_t    nStr[] = { 512,     256  };
+
+  Int_t      q       = (ring == 'I' || ring == 'i') ? 0 : 1;  
+  Double_t   rad     = maxR[q] - minR[q];
+  Double_t   segment = rad / nStr[q];
+  Double_t   r       = minR[q] + segment*strip;
+
+  return r;
+}
+#endif
+
+#if 1
+//_____________________________________________________________________
+Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, 
+                                        UShort_t sec, UShort_t strip, 
+                                        Double_t zvtx)
+{
+  // Calculate eta from strip with vertex (redundant with
+  // AliESDFMD::Eta but support displaced vertices)
+  //
+  // Slightly more optimized version that uses less branching 
   
-  Double_t z = 0;
+  // Get R of the strip
+  Double_t   r         = GetStripR(ring, strip);
+  Int_t      hybrid    = sec / 2;
+  Int_t      q        = (ring == 'I' || ring == 'i') ? 0 : 1;
+
+  const Double_t zs[][2] = { { 320.266, -999999 }, 
+                           {  83.666,  74.966 },
+                           { -63.066, -74.966 } };
+  if (det > 3 || zs[det-1][q] == -999999) return -999999;
+
+  Double_t z = zs[det-1][q];
+  if ((hybrid % 2) == 0) z -= .5;
+  
+  Double_t   theta = TMath::ATan2(r,z-zvtx);
+  Double_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
+  
+  return eta;
+}
+#else
+//_____________________________________________________________________
+Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, 
+                                        UShort_t sec, UShort_t strip, 
+                                        Double_t zvtx)
+{
+  // Calculate eta from strip with vertex (redundant with
+  // AliESDFMD::Eta but support displaced vertices)
+  
+  //Get max R of ring
+  Double_t   r         = GetStripR(ring, strip);
+  Int_t      hybrid    = sec / 2;
+  Bool_t     inner     = (ring == 'I' || ring == 'i');
+  Double_t   z         = 0;
+
+
   switch (det) { 
-  case 1: z = 320.266; break;
-  case 2: z = (inner ? 83.666 : 74.966); break;
+  case 1: z = 320.266;                     break;
+  case 2: z = (inner ?  83.666 :  74.966); break;
   case 3: z = (inner ? -63.066 : -74.966); break; 
   default: return -999999;
   }
@@ -183,9 +560,142 @@ Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec
   
   return eta;
 }
+#endif
 
+//_____________________________________________________________________
+Double_t AliForwardUtil::GetPhiFromStrip(Char_t ring, UShort_t strip, 
+                                        Double_t phi,
+                                        Double_t xvtx, Double_t yvtx)
+{
+  // Calculate eta from strip with vertex (redundant with
+  // AliESDFMD::Eta but support displaced vertices)
 
+  // Unknown x,y -> no change
+  if (yvtx > 999 || xvtx > 999) return phi;
+  
+  //Get max R of ring
+  Double_t r   = GetStripR(ring, strip);
+  Double_t amp = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx) / r;
+  Double_t pha = (TMath::Abs(yvtx) < 1e-12  ? 0 : TMath::ATan2(xvtx, yvtx));
+  Double_t cha = amp * TMath::Cos(phi+pha);
+  phi += cha;
+  if (phi < 0)              phi += TMath::TwoPi();
+  if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
+  return phi;
+}
 //====================================================================
+TAxis*
+AliForwardUtil::MakeFullIpZAxis(Int_t nCenter)
+{
+  TArrayD bins;
+  MakeFullIpZAxis(nCenter, bins);
+  TAxis* a = new TAxis(bins.GetSize()-1,bins.GetArray());
+  return a;
+}
+//____________________________________________________________________
+void
+AliForwardUtil::MakeFullIpZAxis(Int_t nCenter, TArrayD& bins)
+{
+  // Custom vertex axis that will include satellite vertices 
+  // Satellite vertices are at k*37.5 where k=-10,-9,...,9,10 
+  // Nominal vertices are usually in -10 to 10 and we should have 
+  // 10 bins in that range.  That gives us a total of 
+  //
+  //   10+10+10=30 bins 
+  // 
+  // or 31 bin boundaries 
+  if (nCenter % 2 == 1) 
+    // Number of central bins is odd - make it even
+    nCenter--;
+  const Double_t mCenter = 20;
+  const Int_t    nSat    = 10;
+  const Int_t    nBins   = 2*nSat + nCenter;
+  const Int_t    mBin    = nBins / 2;
+  Double_t       dCenter = 2*mCenter / nCenter;
+  bins.Set(nBins+1);
+  bins[mBin] = 0;
+  for (Int_t i = 1; i <= nCenter/2; i++) { 
+    // Assign from the middle out 
+    Double_t  v  = i * dCenter;
+    // Printf("Assigning +/-%7.2f to %3d/%3d", v,mBin-i,mBin+i);
+    bins[mBin-i] = -v;
+    bins[mBin+i] = +v;
+  }
+  for (Int_t i = 1; i <= nSat; i++) { 
+    Double_t v = (i+.5) * 37.5;
+    Int_t    o = nCenter/2+i;
+    // Printf("Assigning +/-%7.2f to %3d/%3d", v,mBin-o,mBin+o);
+    bins[mBin-o] = -v;
+    bins[mBin+o] = +v;
+  }
+}
+//____________________________________________________________________
+void 
+AliForwardUtil::MakeLogScale(Int_t    nBins, 
+                            Int_t    minOrder, 
+                            Int_t    maxOrder, 
+                            TArrayD& bins)
+{
+  Double_t dO = Double_t(maxOrder-minOrder) / nBins; 
+  bins.Set(nBins+1);
+  for (Int_t i = 0; i <= nBins; i++) bins[i] = TMath::Power(10, i * dO+minOrder);
+}
+
+//____________________________________________________________________
+void 
+AliForwardUtil::PrintTask(const TObject& o)
+{
+  Int_t ind = gROOT->GetDirLevel();
+  if (ind > 0) 
+    // Print indention 
+    std::cout << std::setfill(' ') << std::setw(ind) << " " << std::flush;
+
+  TString t = TString::Format("%s %s", o.GetName(), o.ClassName());
+  const Int_t maxN = 75;
+  std::cout << "--- " << t << " " << std::setfill('-') 
+           << std::setw(maxN-ind-5-t.Length()) << "-" << std::endl;
+}
+//____________________________________________________________________
+void
+AliForwardUtil::PrintName(const char* name)
+{
+  Int_t ind = gROOT->GetDirLevel();
+  if (ind > 0) 
+    // Print indention 
+    std::cout << std::setfill(' ') << std::setw(ind) << " " << std::flush;
+    
+  // Now print field name 
+  const Int_t maxN  = 29;
+  Int_t       width = maxN - ind;
+  TString     n(name);
+  if (n.Length() > width-1) {
+    // Truncate the string, and put in "..."
+    n.Remove(width-4);
+    n.Append("...");
+  }
+  n.Append(":");
+  std::cout << std::setfill(' ') << std::left << std::setw(width) 
+           << n << std::right << std::flush;
+}
+//____________________________________________________________________
+void
+AliForwardUtil::PrintField(const char* name, const char* value, ...)
+{
+  PrintName(name);
+
+  // Now format the field value 
+  va_list ap;
+  va_start(ap, value);
+  static char buf[512];
+  vsnprintf(buf, 511, value, ap);
+  buf[511] = '\0';
+  va_end(ap);
+
+  std::cout << buf << std::endl;
+}
+
+//====================================================================
+#if 0 // Moved to separate classes 
 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
 namespace {
@@ -213,6 +723,22 @@ namespace {
     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigmaN);
   }
 
+  Double_t landauGausComposite(Double_t* xp, Double_t* pp)
+  {
+    Double_t x           = xp[0];
+    Double_t cP          = pp[AliForwardUtil::ELossFitter::kC];
+    Double_t deltaP      = pp[AliForwardUtil::ELossFitter::kDelta];
+    Double_t xiP         = pp[AliForwardUtil::ELossFitter::kXi];
+    Double_t sigmaP      = pp[AliForwardUtil::ELossFitter::kSigma];
+    Double_t cS          = pp[AliForwardUtil::ELossFitter::kSigma+1];
+    Double_t deltaS      = deltaP; // pp[AliForwardUtil::ELossFitter::kSigma+2];
+    Double_t xiS         = pp[AliForwardUtil::ELossFitter::kSigma+2/*3*/];
+    Double_t sigmaS      = sigmaP; // pp[AliForwardUtil::ELossFitter::kSigma+4];
+
+    return (cP * AliForwardUtil::LandauGaus(x,deltaP,xiP,sigmaP,0) + 
+           cS * AliForwardUtil::LandauGaus(x,deltaS,xiS,sigmaS,0));
+  }
+    
   // 
   // Utility function to use in TF1 defintition 
   //
@@ -270,7 +796,8 @@ AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
   // Return:
   //    @f$ f'_{L}(x;\Delta,\xi) @f$
   //
-  return TMath::Landau(x, delta - xi * mpshift, xi);
+  Double_t deltaP = delta - xi * mpshift;
+  return TMath::Landau(x, deltaP, xi, true);
 }
 //____________________________________________________________________
 Double_t 
@@ -310,7 +837,9 @@ AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
   // Return:
   //    @f$ f@f$ evaluated at @f$ x@f$.  
   //
-  Double_t deltap = delta - xi * mpshift;
+  if (xi <= 0) return 0;
+
+  Double_t deltaP = delta; // - sigma * sigmaShift; // + sigma * mpshift;
   Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
   Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
@@ -322,12 +851,35 @@ AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
     Double_t x1 = xlow  + (i - .5) * step;
     Double_t x2 = xhigh - (i - .5) * step;
     
-    sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
-    sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
+    //sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
+    //sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
+    sum += Landau(x1, deltaP, xi) * TMath::Gaus(x, x1, sigma1);
+    sum += Landau(x2, deltaP, xi) * TMath::Gaus(x, x2, sigma1);
   }
   return step * sum * invSq2pi / sigma1;
 }
 
+namespace { 
+  const Double_t sigmaShift = 0.36390; // TMath::Log(TMath::Sqrt(2.));
+  double deltaSigmaShift(Int_t i, Double_t sigma)
+  {
+    return 0; // - sigma * sigmaShift;
+  }
+  void getIPars(Int_t i, Double_t& delta, Double_t& xi, Double_t& sigma)
+  {
+    Double_t dsig = deltaSigmaShift(i, sigma);
+    if (i == 1) { 
+      delta += dsig;
+      return; // { delta = delta + xi*mpshift; return; } // Do nothing 
+    }
+    
+    delta = i * (delta + xi * TMath::Log(i)) + dsig;
+    xi    = i * xi;
+    sigma = TMath::Sqrt(Double_t(i)) * sigma;
+  }
+}
+    
+    
 //____________________________________________________________________
 Double_t 
 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi, 
@@ -357,9 +909,10 @@ AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi,
   // Return:
   //    @f$ f_i @f$ evaluated
   //  
-  Double_t deltaI =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
-  Double_t xiI    =  i * xi;
-  Double_t sigmaI =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
+  Double_t deltaI = delta;
+  Double_t xiI    = xi;
+  Double_t sigmaI = sigma;
+  getIPars(i, deltaI, xiI, sigmaI);
   if (sigmaI < 1e-10) { 
     // Fall back to landau 
     return AliForwardUtil::Landau(x, deltaI, xiI);
@@ -610,7 +1163,7 @@ AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut,
                                         Double_t maxRange, 
                                         UShort_t minusBins) 
   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins), 
-    fFitResults(0), fFunctions(0)
+    fFitResults(0), fFunctions(0), fDebug(false)
 {
   // 
   // Constructor 
@@ -644,6 +1197,19 @@ AliForwardUtil::ELossFitter::Clear()
   fFitResults.Clear();
   fFunctions.Clear();
 }
+namespace { 
+  void setParLimit(TF1* f, Int_t iPar, Bool_t debug,
+                  Double_t test, Double_t low, Double_t high)
+  {
+    if (test >= low && test <= high) { 
+      if (debug) 
+       printf("Fit: Set par limits on %s: %f, %f\n", 
+              f->GetParName(iPar), low, high);
+      f->SetParLimits(iPar, low, high);
+    }
+  }
+}
+
 //____________________________________________________________________
 TF1*
 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
@@ -667,37 +1233,64 @@ AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
   Clear();
   
   // Find the fit range 
-  dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
+  // Find the fit range 
+  Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
+  Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
+                               dist->GetNbinsX());
+  dist->GetXaxis()->SetRange(cutBin, maxBin);
+  // dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
   
   // Get the bin with maximum 
-  Int_t    maxBin = dist->GetMaximumBin();
-  Double_t maxE   = dist->GetBinLowEdge(maxBin);
+  Int_t    peakBin = dist->GetMaximumBin();
+  Double_t peakE   = dist->GetBinLowEdge(peakBin);
+  Double_t rmsE    = dist->GetRMS();
   
   // Get the low edge 
-  dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
-  Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
+  // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
+  Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
-  Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
+  Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
 
+  Int_t    minEb = dist->GetXaxis()->FindBin(minE);
+  Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
+  Double_t intg  = dist->Integral(minEb, maxEb);
+  if (intg <= 0) {
+    ::Warning("Fit1Particle", 
+             "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0", 
+             dist->GetName(), minE, maxE, minEb, maxEb, intg);
+    return 0;
+  }
+    
   // Restore the range 
-  dist->GetXaxis()->SetRangeUser(0, fMaxRange);
+  dist->GetXaxis()->SetRange(1, maxBin);
   
   // Define the function to fit 
-  TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
+  TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxE,kSigmaN+1);
 
   // Set initial guesses, parameter names, and limits  
-  landau1->SetParameters(1,0.5,0.07,0.1,sigman);
+  landau1->SetParameters(intg,peakE,peakE/10,peakE/5,sigman);
   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
   landau1->SetNpx(500);
-  landau1->SetParLimits(kDelta, minE, fMaxRange);
-  landau1->SetParLimits(kXi,    0.00, fMaxRange);
-  landau1->SetParLimits(kSigma, 0.01, 0.1);
+  setParLimit(landau1, kDelta, fDebug, peakE,   minE, fMaxRange);
+  setParLimit(landau1, kXi,    fDebug, peakE,   0,    rmsE); // 0.1
+  setParLimit(landau1, kSigma, fDebug, peakE/5, 1e-5, rmsE); // 0.1
   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
-  else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
+  else 
+    setParLimit(landau1, kSigmaN, fDebug, peakE, 0, rmsE);
+  
 
+  TString opts(Form("%s%s", FIT_OPTIONS, fDebug ? "" : "Q"));
   // Do the fit, getting the result object 
-  TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
-  landau1->SetRange(minE, fMaxRange);
+  if (fDebug) 
+    ::Info("Fit1Particle", "Fitting in the range %f,%f", minE, maxE);
+  TFitResultPtr r = dist->Fit(landau1, opts, "", minE, maxE);
+  if (!r.Get()) { 
+    ::Warning("Fit1Particle", 
+             "No fit returned when processing %s in the range [%f,%f] "
+             "options %s", dist->GetName(), minE, maxE, FIT_OPTIONS);
+    return 0;
+  }
+  // landau1->SetRange(minE, fMaxRange);
   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
   fFunctions.AddAtAndExpand(landau1, 0);
 
@@ -742,38 +1335,194 @@ AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
   Double_t minE    = f->GetXmin();
 
+  Int_t    minEb = dist->GetXaxis()->FindBin(minE);
+  Int_t    maxEb = dist->GetXaxis()->FindBin(maxEi);
+  Double_t rmsE  = dist->GetRMS();
+  Double_t intg  = dist->Integral(minEb, maxEb);
+  if (intg <= 0) {
+    ::Warning("FitNParticle",
+             "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0", 
+             dist->GetName(), minE, maxEi, minEb, maxEb, intg);
+    return 0;
+  }
+
   // Array of weights 
   TArrayD a(n-1);
   for (UShort_t i = 2; i <= n; i++) 
     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
   // Make the fit function 
-  TF1* landaun     = MakeNLandauGaus(r->Parameter(kC),
-                                    r->Parameter(kDelta),
-                                    r->Parameter(kXi),
-                                    r->Parameter(kSigma),
-                                    r->Parameter(kSigmaN),
-                                    n,a.fArray,minE,maxEi);
-  landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
-  landaun->SetParLimits(kXi,     0.00, fMaxRange);       // xi
-  landaun->SetParLimits(kSigma,  0.01, 1);            // sigma
-  // Check if we're using the noise sigma 
+  TF1* landaun = MakeNLandauGaus(r->Parameter(kC),
+                                r->Parameter(kDelta),
+                                r->Parameter(kXi),
+                                r->Parameter(kSigma),
+                                r->Parameter(kSigmaN),
+                                n, a.fArray, minE, maxEi);
+  setParLimit(landaun, kDelta, fDebug, r->Parameter(kDelta), minE, fMaxRange);
+  setParLimit(landaun, kXi,    fDebug, r->Parameter(kXi),    0,    rmsE); // 0.1
+  setParLimit(landaun, kSigma, fDebug, r->Parameter(kSigma), 1e-5, rmsE); // 0.1
   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
-  else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
+  else 
+    setParLimit(landaun, kSigmaN, fDebug, r->Parameter(kSigmaN), 0, rmsE);
 
   // Set the range and name of the scale parameters 
   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
-    landaun->SetParLimits(kA+i-2, 0,1);
+    setParLimit(landaun, kA+i-2, fDebug, a[i-2], 0, 1);
   }
 
   // Do the fit 
-  TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
+  TString opts(Form("%s%s", FIT_OPTIONS, fDebug ? "" : "Q"));
+  if (fDebug) 
+    ::Info("FitNParticle", "Fitting in the range %f,%f (%d)", minE, maxEi, n);
+  TFitResultPtr tr = dist->Fit(landaun, opts, "", minE, maxEi);
   
-  landaun->SetRange(minE, fMaxRange);
+  // landaun->SetRange(minE, fMaxRange);
   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
   fFunctions.AddAtAndExpand(landaun, n-1);
   
   return landaun;
 }  
+//____________________________________________________________________
+TF1*
+AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
+{
+  // 
+  // Fit a composite particle signal to the passed energy loss
+  // distribution  // 
+  // Parameters:
+  //    dist    Data to fit the function to 
+  //    sigman If larger than zero, the initial guess of the
+  //               detector induced noise. If zero or less, then this 
+  //               parameter is ignored in the fit (fixed at 0)
+  // 
+  // Return:
+  //    The function fitted to the data 
+  //
+
+  // Find the fit range 
+  Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
+  Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
+                               dist->GetNbinsX());
+  dist->GetXaxis()->SetRange(cutBin, maxBin);
+  
+  // Get the bin with maximum 
+  Int_t    peakBin = dist->GetMaximumBin();
+  Double_t peakE   = dist->GetBinLowEdge(peakBin);
+  
+  // Get the low edge 
+  // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
+  Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
+  Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
+  Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
+
+  // Get the range in bins and the integral of that range 
+  Int_t    minEb = dist->GetXaxis()->FindBin(minE);
+  Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
+  Double_t intg  = dist->Integral(minEb, maxEb);
+  if (intg <= 0) {
+    ::Warning("Fit1Particle", 
+             "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0", 
+             dist->GetName(), minE, maxE, minEb, maxEb, intg);
+    return 0;
+  }
+    
+  // Restore the range 
+  dist->GetXaxis()->SetRange(1, maxBin);
+  
+  // Define the function to fit 
+  TF1* seed = new TF1("landauSeed", landauGaus1, minE,maxE,kSigmaN+1);
+
+  // Set initial guesses, parameter names, and limits  
+  seed->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
+  seed->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
+  seed->SetNpx(500);
+  seed->SetParLimits(kDelta, minE, fMaxRange);
+  seed->SetParLimits(kXi,    0.00, 0.1); // Was fMaxRange - too wide
+  seed->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
+  if (sigman <= 0)  seed->FixParameter(kSigmaN, 0);
+  else              seed->SetParLimits(kSigmaN, 0, fMaxRange);
+
+  // Do the fit, getting the result object 
+  if (fDebug) 
+    ::Info("FitComposite", "Fitting seed in the range %f,%f", minE, maxE);
+  /* TFitResultPtr r = */ dist->Fit(seed, FIT_OPTIONS, "", minE, maxE);
+
+  maxE = dist->GetXaxis()->GetXmax();
+#if 1
+  TF1* comp = new TF1("composite", landauGausComposite, 
+                     minE, maxE, kSigma+1+2);
+  comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
+                   "C#prime", "#xi#prime");
+  comp->SetParameters(0.8 * seed->GetParameter(kC),  // 0 Primary weight 
+                     seed->GetParameter(kDelta),    // 1 Primary Delta
+                     seed->GetParameter(kDelta)/10, // 2 primary Xi
+                     seed->GetParameter(kDelta)/5,  // 3 primary sigma
+                     1.20 * seed->GetParameter(kC), // 5 Secondary weight
+                     seed->GetParameter(kXi));      // 7 secondary Xi
+  // comp->SetParLimits(kC,       minE, fMaxRange); // C
+  comp->SetParLimits(kDelta,      minE, fMaxRange); // Delta
+  comp->SetParLimits(kXi,         0.00, fMaxRange); // Xi 
+  comp->SetParLimits(kSigma,      1e-5, fMaxRange); // Sigma
+  // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
+  comp->SetParLimits(kSigma+2,    0.00, fMaxRange); // Xi'
+#else
+  TF1* comp = new TF1("composite", landauGausComposite, 
+                     minE, maxE, kSigma+1+4);
+  comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
+                   "C#prime", "#Delta_{p}#prime", "#xi#prime", "#sigma#prim");
+  comp->SetParameters(0.8 * seed->GetParameter(kC),  // 0 Primary weight 
+                     seed->GetParameter(kDelta),    // 1 Primary Delta
+                     seed->GetParameter(kDelta)/10, // 2 primary Xi
+                     seed->GetParameter(kDelta)/5,  // 3 primary sigma
+                     1.20 * seed->GetParameter(kC), // 5 Secondary weight
+                     seed->GetParameter(kDelta),    // 6 secondary Delta
+                     seed->GetParameter(kXi),       // 7 secondary Xi
+                     seed->GetParameter(kSigma));   // 8 secondary sigma
+  // comp->SetParLimits(kC,       minE, fMaxRange); // C
+  comp->SetParLimits(kDelta,      minE, fMaxRange); // Delta
+  comp->SetParLimits(kXi,         0.00, fMaxRange); // Xi 
+  comp->SetParLimits(kSigma,      1e-5, fMaxRange); // Sigma
+  // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
+  comp->SetParLimits(kSigma+2,    minE/10, fMaxRange); // Delta
+  comp->SetParLimits(kSigma+3,    0.00,    fMaxRange); // Xi 
+  comp->SetParLimits(kSigma+4,    1e-6,    fMaxRange); // Sigma
+#endif               
+  comp->SetLineColor(kRed+1);
+  comp->SetLineWidth(3);
+  
+  // Do the fit, getting the result object 
+  TString opts(Form("%s%s", FIT_OPTIONS, fDebug ? "" : "Q"));
+  if (fDebug) 
+    ::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
+  /* TFitResultPtr r = */ dist->Fit(comp, opts, "", minE, maxE);
+
+#if 0
+  TF1* part1 = static_cast<TF1*>(seed->Clone("part1"));
+  part1->SetLineColor(kGreen+1);
+  part1->SetLineWidth(4);
+  part1->SetRange(minE, maxE);
+  part1->SetParameters(comp->GetParameter(0), // C 
+                      comp->GetParameter(1), // Delta
+                      comp->GetParameter(2), // Xi
+                      comp->GetParameter(3), // sigma
+                      0);
+  part1->Save(minE,maxE,0,0,0,0);
+  dist->GetListOfFunctions()->Add(part1);
+
+  TF1* part2 = static_cast<TF1*>(seed->Clone("part2"));
+  part2->SetLineColor(kBlue+1);
+  part2->SetLineWidth(4);
+  part2->SetRange(minE, maxE);
+  part2->SetParameters(comp->GetParameter(4), // C 
+                      comp->GetParameter(5), // Delta
+                      comp->GetParameter(6), // Xi
+                      comp->GetParameter(7), // sigma
+                      0);
+  part2->Save(minE,maxE,0,0,0,0);
+  dist->GetListOfFunctions()->Add(part2);
+#endif
+  return comp;
+}
+#endif
 
 //====================================================================
 AliForwardUtil::Histos::~Histos()
@@ -781,17 +1530,28 @@ AliForwardUtil::Histos::~Histos()
   // 
   // Destructor
   //
+}
+
+//____________________________________________________________________
+void
+AliForwardUtil::Histos::Delete(Option_t* opt)
+{
   if (fFMD1i) delete fFMD1i;
   if (fFMD2i) delete fFMD2i;
   if (fFMD2o) delete fFMD2o;
   if (fFMD3i) delete fFMD3i;
   if (fFMD3o) delete fFMD3o;
+  fFMD1i = 0;
+  fFMD2i = 0;
+  fFMD2o = 0;
+  fFMD3i = 0;
+  fFMD3o = 0;
+  TObject::Delete(opt);
 }
 
 //____________________________________________________________________
 TH2D*
-AliForwardUtil::Histos::Make(UShort_t d, Char_t r, 
-                            const TAxis& etaAxis) const
+AliForwardUtil::Histos::Make(UShort_t d, Char_t r, const TAxis& etaAxis)
 {
   // 
   // Make a histogram 
@@ -805,10 +1565,17 @@ AliForwardUtil::Histos::Make(UShort_t d, Char_t r,
   //    Newly allocated histogram 
   //
   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
-  TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r), 
-                       Form("FMD%d%c cache", d, r),
-                       etaAxis.GetNbins(), etaAxis.GetXmin(), 
-                       etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
+  TH2D* hist = 0;
+  if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
+    hist = new TH2D(Form("FMD%d%c_cache", d, r), 
+                   Form("FMD%d%c cache", d, r),
+                   etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray(), 
+                   ns, 0, TMath::TwoPi());
+  else
+    hist = new TH2D(Form("FMD%d%c_cache", d, r), 
+                   Form("FMD%d%c cache", d, r),
+                   etaAxis.GetNbins(), etaAxis.GetXmin(), 
+                   etaAxis.GetXmax(), ns, 0, TMath::TwoPi());
   hist->SetXTitle("#eta");
   hist->SetYTitle("#phi [radians]");
   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
@@ -817,6 +1584,19 @@ AliForwardUtil::Histos::Make(UShort_t d, Char_t r,
 
   return hist;
 }
+//____________________________________________________________________
+void
+AliForwardUtil::Histos::RebinEta(TH2D* hist, const TAxis& etaAxis)
+{
+  TAxis* xAxis = hist->GetXaxis();
+  if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
+    xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray());
+  else
+    xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax());
+  hist->Rebuild();
+}
+  
+
 //____________________________________________________________________
 void
 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
@@ -833,6 +1613,23 @@ AliForwardUtil::Histos::Init(const TAxis& etaAxis)
   fFMD3i = Make(3, 'I', etaAxis);
   fFMD3o = Make(3, 'O', etaAxis);
 }
+//____________________________________________________________________
+void
+AliForwardUtil::Histos::ReInit(const TAxis& etaAxis)
+{
+  // 
+  // Initialize the object 
+  // 
+  // Parameters:
+  //    etaAxis Eta axis to use 
+  //
+  if (!fFMD1i) fFMD1i = Make(1, 'i', etaAxis); else RebinEta(fFMD1i, etaAxis);
+  if (!fFMD2i) fFMD2i = Make(2, 'i', etaAxis); else RebinEta(fFMD2i, etaAxis);
+  if (!fFMD2o) fFMD2o = Make(2, 'o', etaAxis); else RebinEta(fFMD2o, etaAxis);
+  if (!fFMD3i) fFMD3i = Make(3, 'i', etaAxis); else RebinEta(fFMD3i, etaAxis);
+  if (!fFMD3o) fFMD3o = Make(3, 'o', etaAxis); else RebinEta(fFMD3o, etaAxis);
+}
+
 //____________________________________________________________________
 void
 AliForwardUtil::Histos::Clear(Option_t* option)
@@ -843,11 +1640,11 @@ AliForwardUtil::Histos::Clear(Option_t* option)
   // Parameters:
   //    option Not used 
   //
-  if (fFMD1i) fFMD1i->Reset(option);
-  if (fFMD2i) fFMD2i->Reset(option);
-  if (fFMD2o) fFMD2o->Reset(option);
-  if (fFMD3i) fFMD3i->Reset(option);
-  if (fFMD3o) fFMD3o->Reset(option);
+  if (fFMD1i) { fFMD1i->Reset(option); fFMD1i->ResetBit(kSkipRing); }
+  if (fFMD2i) { fFMD2i->Reset(option); fFMD2i->ResetBit(kSkipRing); }
+  if (fFMD2o) { fFMD2o->Reset(option); fFMD2o->ResetBit(kSkipRing); }
+  if (fFMD3i) { fFMD3i->Reset(option); fFMD3i->ResetBit(kSkipRing); }
+  if (fFMD3o) { fFMD3o->Reset(option); fFMD3o->ResetBit(kSkipRing); }
 }
 
 //____________________________________________________________________
@@ -934,29 +1731,49 @@ AliForwardUtil::DebugGuard::DebugGuard(Int_t lvl, Int_t msgLvl,
   if (lvl < msgLvl) return; 
   va_list ap;
   va_start(ap, format);
-  static char buf[512];
-  Int_t n = gROOT->GetDirLevel() + 1;
-  for (Int_t i = 0; i < n; i++) buf[i] = ' ';
-  vsnprintf(&(buf[n]), 511-n, format, ap);
-  buf[511] = '\0';
-  fMsg = buf;
+  Format(fMsg, format, ap);
   va_end(ap);
-  Format(true);
+  Output(+1, fMsg);
 }
 //____________________________________________________________________
 AliForwardUtil::DebugGuard::~DebugGuard()
 {
   if (fMsg.IsNull()) return;
-  Format(false);
+  Output(-1, fMsg);
+}
+//____________________________________________________________________
+void
+AliForwardUtil::DebugGuard::Message(Int_t lvl, Int_t msgLvl, 
+                                   const char* format, ...)
+{
+  if (lvl < msgLvl) return; 
+  TString msg;
+  va_list ap;
+  va_start(ap, format);
+  Format(msg, format, ap);
+  va_end(ap);
+  Output(0, msg);
+}
+
+//____________________________________________________________________
+void
+AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
+{
+  static char buf[512];
+  Int_t n = gROOT->GetDirLevel() + 2;
+  for (Int_t i = 0; i < n; i++) buf[i] = ' ';
+  vsnprintf(&(buf[n]), 511-n, format, ap);
+  buf[511] = '\0';
+  out = buf;  
 }
 //____________________________________________________________________
 void
-AliForwardUtil::DebugGuard::Format(bool in)
+AliForwardUtil::DebugGuard::Output(int in, TString& msg)
 {
-  fMsg[0] = in ? '>' : '<';
-  AliLog::Debug(0, fMsg, "PWGLF-forward", "", "", "", 0);
-  if (in) gROOT->IncreaseDirLevel();
-  else    gROOT->DecreaseDirLevel();
+  msg[0] = (in > 0 ? '>' :  in < 0 ? '<' : '=');
+  AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
+  if      (in > 0) gROOT->IncreaseDirLevel();
+  else if (in < 0) gROOT->DecreaseDirLevel();
 }