]> 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 c7cfccefa621ffd49e29344727feae8050822f04..ac7176862be852557305db5c1342a3f59aa0c6d2 100644 (file)
 #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;
-#else 
+#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()
 {
-#ifdef ALIROOT_SVN_BRANCH
+  // 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(ALIROOT_SVN_BRANCH);
+  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("trunk")) return ret = 0xFFFFFFFF;
+  if (str.EqualTo(top)) return ret = 0xFFFFFFFF;
 
   TObjArray*   tokens = str.Tokenize("-");
-  TObjString*  pMajor = static_cast<TObjString*>(tokens->At(0));
-  TObjString*  pMinor = static_cast<TObjString*>(tokens->At(1));
-  TObjString*  pRelea = static_cast<TObjString*>(tokens->At(2));
-  TObjString* pAn     = (tokens->GetEntries() > 3 ? 
-    static_cast<TObjString*>(tokens->At(3)) : 0);
-  TString sMajor = pMajor->String().Strip(TString::kLeading, '0');
-  TString sMinor = pMinor->String().Strip(TString::kLeading, '0');
-  TString sRelea = pRelea->String().Strip(TString::kLeading, '0');
-  
+  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;
-#else 
-  return 0;
-#endif
 }
 
 //====================================================================
@@ -159,6 +201,22 @@ AliForwardUtil::CenterOfMassRapidity(UShort_t z1,
   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 beam)
@@ -181,17 +239,12 @@ AliForwardUtil::ParseCenterOfMassEnergy(UShort_t sys, Float_t beam)
     energy = CenterOfMassEnergy(beam, 82, 208, 1, 1);
   else if (sys == AliForwardUtil::kPbPb) 
     energy = CenterOfMassEnergy(beam, 82, 208, 82, 208);
-  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 ret = CheckSNN(energy);
+  if (ret > 1) return ret;
+  if (sys == AliForwardUtil::kPbPb || sys == AliForwardUtil::kPPb) {
+    ret = CheckSNN(beam);
+  }
+  return ret;
 }
 //____________________________________________________________________
 const char* 
@@ -261,11 +314,11 @@ UShort_t AliForwardUtil::CheckForAOD()
 {
   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
   if (dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler())) {
-    ::Info("CheckForAOD", "Found AOD Input handler");
+    // ::Info("CheckForAOD", "Found AOD Input handler");
     return 1;
   }
   if (dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler())) {
-    ::Info("CheckForAOD", "Found AOD Output handler");
+    // ::Info("CheckForAOD", "Found AOD Output handler");
     return 2;
   }
 
@@ -352,7 +405,7 @@ void AliForwardUtil::GetParameter(TObject* o, UShort_t& value)
 {
   if (!o) return;
   TParameter<int>* p = static_cast<TParameter<int>*>(o);
-  if (p->TestBit(BIT(17)))
+  if (p->TestBit(BIT(19)))
     value = p->GetVal(); 
   else
     value = o->GetUniqueID();
@@ -362,7 +415,7 @@ void AliForwardUtil::GetParameter(TObject* o, Int_t& value)
 {
   if (!o) return;
   TParameter<int>* p = static_cast<TParameter<int>*>(o);
-  if (p->TestBit(BIT(17)))
+  if (p->TestBit(BIT(19)))
     value = p->GetVal(); 
   else
     value = o->GetUniqueID();
@@ -372,7 +425,7 @@ 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(17)))
+  if (p->TestBit(BIT(19)))
     value = p->GetVal(); 
   else
     value = o->GetUniqueID();
@@ -382,7 +435,7 @@ void AliForwardUtil::GetParameter(TObject* o, Double_t& value)
 {
   if (!o) return;
   TParameter<double>* p = static_cast<TParameter<double>*>(o);
-  if (p->TestBit(BIT(17)))
+  if (p->TestBit(BIT(19)))
     value = p->GetVal(); // o->GetUniqueID();
   else {
     UInt_t  i = o->GetUniqueID();
@@ -395,7 +448,7 @@ void AliForwardUtil::GetParameter(TObject* o, Bool_t& value)
 {
   if (!o) return;
   TParameter<bool>* p = static_cast<TParameter<bool>*>(o);
-  if (p->TestBit(BIT(17)))
+  if (p->TestBit(BIT(19)))
     value = p->GetVal(); // o->GetUniqueID();
   else
     value = o->GetUniqueID();
@@ -530,8 +583,119 @@ Double_t AliForwardUtil::GetPhiFromStrip(Char_t ring, UShort_t strip,
   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 {
@@ -567,9 +731,9 @@ namespace {
     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      = pp[AliForwardUtil::ELossFitter::kSigma+2];
-    Double_t xiS         = pp[AliForwardUtil::ELossFitter::kSigma+3];
-    Double_t sigmaS      = pp[AliForwardUtil::ELossFitter::kSigma+4];
+    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));
@@ -632,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 
@@ -672,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;
@@ -684,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, 
@@ -719,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);
@@ -972,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 
@@ -1006,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)
@@ -1039,6 +1243,7 @@ AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
   // Get the bin with maximum 
   Int_t    peakBin = dist->GetMaximumBin();
   Double_t peakE   = dist->GetBinLowEdge(peakBin);
+  Double_t rmsE    = dist->GetRMS();
   
   // Get the low edge 
   // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
@@ -1063,18 +1268,28 @@ AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxE,kSigmaN+1);
 
   // Set initial guesses, parameter names, and limits  
-  landau1->SetParameters(1,peakE,peakE/10,peakE/5,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, 0.1); // Was fMaxRange - too wide
-  landau1->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
+  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 
-  ::Info("Fit1Particle", "Fitting in the range %f,%f", minE, maxE);
-  TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxE);
+  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);
@@ -1122,6 +1337,7 @@ AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
 
   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",
@@ -1141,21 +1357,23 @@ AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
                                 r->Parameter(kSigma),
                                 r->Parameter(kSigmaN),
                                 n, a.fArray, minE, maxEi);
-  landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
-  landaun->SetParLimits(kXi,     0.00, 0.1);  // was fMaxRange - too wide
-  landaun->SetParLimits(kSigma,  1e-5, 0.1);  // was fMaxRange - too wide
-  // Check if we're using the noise sigma 
+  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 
-  ::Info("Fit1Particle", "Fitting in the range %f,%f", minE, maxEi);
-  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);
   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
@@ -1169,8 +1387,7 @@ AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
 {
   // 
   // Fit a composite particle signal to the passed energy loss
-  // distribution
-  // 
+  // distribution  // 
   // Parameters:
   //    dist    Data to fit the function to 
   //    sigman If larger than zero, the initial guess of the
@@ -1225,10 +1442,29 @@ AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
   else              seed->SetParLimits(kSigmaN, 0, fMaxRange);
 
   // Do the fit, getting the result object 
-  ::Info("FitComposite", "Fitting seed in the range %f,%f", minE, maxE);
-  /* TFitResultPtr r = */ dist->Fit(seed, "RNQS", "", minE, maxE);
+  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",
@@ -1241,7 +1477,6 @@ AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
                      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 
@@ -1250,12 +1485,15 @@ AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
   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 
-  ::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
-  /* TFitResultPtr r = */ dist->Fit(comp, "RNQS", "", minE, maxE);
+  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"));
@@ -1284,6 +1522,7 @@ AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
 #endif
   return comp;
 }
+#endif
 
 //====================================================================
 AliForwardUtil::Histos::~Histos()
@@ -1384,11 +1623,11 @@ AliForwardUtil::Histos::ReInit(const TAxis& etaAxis)
   // Parameters:
   //    etaAxis Eta axis to use 
   //
-  RebinEta(fFMD1i, etaAxis);
-  RebinEta(fFMD2i, etaAxis);
-  RebinEta(fFMD2o, etaAxis);
-  RebinEta(fFMD3i, etaAxis);
-  RebinEta(fFMD3o, etaAxis);
+  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);
 }
 
 //____________________________________________________________________
@@ -1401,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); }
 }
 
 //____________________________________________________________________