]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/FORWARD/analysis2/AliForwardUtil.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardUtil.cxx
index c950fe58cae7bf76a5b518b4d817989eba93be0c..4d8b87e7fe2f6f184564b84296cdfffe99e42fbc 100644 (file)
 #include <TMath.h>
 #include <TError.h>
 #include <TROOT.h>
+#define FIT_OPTIONS "RNQS"
 
 //====================================================================
 ULong_t AliForwardUtil::AliROOTRevision()
 {
 #ifdef ALIROOT_SVN_REVISION
   return ALIROOT_SVN_REVISION;
-#else 
-  return 0;
+#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;
 #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[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));
@@ -61,9 +99,6 @@ ULong_t AliForwardUtil::AliROOTBranch()
     (pAn ? 0xAA : 0));
   
   return ret;
-#else 
-  return 0;
-#endif
 }
 
 //====================================================================
@@ -159,6 +194,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 +232,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 +307,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;
   }
 
@@ -332,10 +378,10 @@ TObject* AliForwardUtil::MakeParameter(const Char_t* name, ULong_t value)
 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);
+  // Float_t v = value;
+  // UInt_t* tmp = reinterpret_cast<UInt_t*>(&v);
   ret->SetMergeMode('f');
-  ret->SetUniqueID(*tmp);
+  // ret->SetUniqueID(*tmp);
   return ret;
 }
 //_____________________________________________________________________
@@ -351,54 +397,134 @@ TObject* AliForwardUtil::MakeParameter(const Char_t* name, Bool_t value)
 void AliForwardUtil::GetParameter(TObject* o, UShort_t& value)
 {
   if (!o) return;
-  value = o->GetUniqueID();
+  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;
-  value = o->GetUniqueID();
+  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;
-  value = o->GetUniqueID();
+  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;
-  UInt_t  i = o->GetUniqueID();
-  Float_t v = *reinterpret_cast<Float_t*>(&i);
-  value = v;
+  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;
-  value = o->GetUniqueID();
+  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
-  const Double_t iMinR = 4.5213;
-  const Double_t iMaxR = 17.2;
-  const Double_t oMinR = 15.4;
-  const Double_t oMaxR = 28.0;
-  
-  Double_t   minR    = (ring == 'I' || ring == 'i') ? iMinR : oMinR;
-  Double_t   maxR    = (ring == 'I' || ring == 'i') ? iMaxR : oMaxR;
-  Double_t   nStrips = (ring == 'I' || ring == 'i') ? 512   : 256;
-  Double_t   rad     =  maxR - minR;
-  Double_t   segment = rad / nStrips;
-  Double_t   r       =  minR + segment*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 
+  
+  // 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, 
@@ -412,6 +538,8 @@ Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring,
   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;
@@ -425,6 +553,7 @@ Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring,
   
   return eta;
 }
+#endif
 
 //_____________________________________________________________________
 Double_t AliForwardUtil::GetPhiFromStrip(Char_t ring, UShort_t strip, 
@@ -447,6 +576,111 @@ 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);
+}
+
+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;
+}
 
 //====================================================================
 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
@@ -484,9 +718,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));
@@ -889,7 +1123,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 
@@ -983,15 +1217,34 @@ AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
   landau1->SetParameters(1,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
+  if (peakE >= minE && peakE <= fMaxRange) {
+    // printf("Fit1: Set par limits on Delta: %f, %f\n", minE, fMaxRange);
+    landau1->SetParLimits(kDelta, minE, fMaxRange);
+  }
+  if (peakE/10 >= 0 && peakE <= 0.1) {
+    // printf("Fit1: Set par limits on xi: %f, %f\n", 0., 0.1);
+    landau1->SetParLimits(kXi,    0.00, 0.1); // Was fMaxRange - too wide
+  }
+  if (peakE/5 >= 0 && peakE/5 <= 0.1) {
+    // printf("Fit1: Set par limits on sigma: %f, %f\n", 0., 0.1);
+    landau1->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
+  }
   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
-  else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
+  else {
+    // printf("Fit1: Set par limits on sigmaN: %f, %f\n", 0., fMaxRange);
+    landau1->SetParLimits(kSigmaN, 0, fMaxRange);
+  }
 
   // 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, FIT_OPTIONS, "", 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);
@@ -1058,21 +1311,39 @@ 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
+  if (minE      <= r->Parameter(kDelta) &&
+      fMaxRange >= r->Parameter(kDelta)) {
+    // Protect against warning from ParameterSettings
+    // printf("FitN: Set par limits on Delta: %f, %f\n", minE, fMaxRange);
+    landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
+  }
+  if (r->Parameter(kXi) >= 0 && r->Parameter(kXi) <= 0.1) {
+    // printf("FitN: Set par limits on xi: %f, %f\n", 0., 0.1);
+    landaun->SetParLimits(kXi,     0.00, 0.1);  // was fMaxRange - too wide
+  }
+  if (r->Parameter(kSigma) >= 1e-5 && r->Parameter(kSigma) <= 0.1) {
+    // printf("FitN: Set par limits on sigma: %f, %f\n", 1e-5, 0.1);
+    landaun->SetParLimits(kSigma,  1e-5, 0.1);  // was fMaxRange - too wide
+  }
   // Check if we're using the noise sigma 
   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
-  else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
+  else {
+    // printf("FitN: Set par limits on sigmaN: %f, %f\n", 0., fMaxRange);
+    landaun->SetParLimits(kSigmaN, 0, fMaxRange);
+  }
 
   // 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);
+    if (a[i-2] >= 0 && a[i-2] <= 1) {
+      // printf("FitN: Set par limits on a_%d: %f, %f\n", i, 0., 1.);
+      landaun->SetParLimits(kA+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);
+  if (fDebug) 
+    ::Info("FitNParticle", "Fitting in the range %f,%f (%d)", minE, maxEi, n);
+  TFitResultPtr tr = dist->Fit(landaun, FIT_OPTIONS, "", minE, maxEi);
   
   // landaun->SetRange(minE, fMaxRange);
   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
@@ -1142,10 +1413,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",
@@ -1158,7 +1448,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 
@@ -1167,12 +1456,14 @@ 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);
+  if (fDebug) 
+    ::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
+  /* TFitResultPtr r = */ dist->Fit(comp, FIT_OPTIONS, "", minE, maxE);
 
 #if 0
   TF1* part1 = static_cast<TF1*>(seed->Clone("part1"));
@@ -1301,11 +1592,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);
 }
 
 //____________________________________________________________________
@@ -1318,11 +1609,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); }
 }
 
 //____________________________________________________________________