ANALYSIS/TFndFitter.cxx

// @(#)fROOT/ANALYSYS:$Name:  $:$Id: TFndFitter.cxx,v 1.11 2007/09/24 07:32:41 Diego_Faso Exp $
// Author: Diego Faso <mailto:faso@to.infn.it>, 11/05/07

//                                                  //
//                    TFndFitter                    //
//                                                  //
//  This class provides a general-purposes          //
//  fit-environment, based on the FIN_AN namespace  //
//                                                  //

#include "TFndFitter.h"
#include "TROOT.h"

ClassImp(TFndFitter)

//_______________________
TFndFitter::TFndFitter():
  fCurHist(),fCurSpec(),
  fResolutionFactor(1.)
{
  
}

//_______________________
TFndFitter::~TFndFitter(){
  //   delete fCurHist;
  //   delete fCurSpec;
}

//_____________________
Int_t TFndFitter::FitAllPeaks_GausYield(Double_t sigma_start,Double_t sigma_maxerr){
  // return value
  //             0: ok
  //            -1: at least one error
  
  Int_t retres = 0;

  Int_t n_found_peaks = fCurSpec->GetNPeaks();
  Printf("Found %d candidate peaks to fit\n",n_found_peaks);

  for(Int_t peak_num=1;peak_num<=n_found_peaks;peak_num++){
    Int_t fres = FitSinglePeak_GausYield(peak_num,sigma_start,sigma_maxerr);
    if (fres != 0) retres = -1;
  }
  
  return retres;
}

//_____________________
Int_t TFndFitter::FitSinglePeak_GausYield(Int_t peak_num,Double_t sigma_init,Double_t sigma_maxerr){
  // this methods uses peaks already evaluated (by FIN_AN::FindPeaks)
  // return value (fit_result)
  //             0: ok
  //           !=0: see root-documentation  
  //          -999: error within this method

  if(!fCurSpec){
    Warning("FitSinglePeak","No spectrum selected: nothing will happen...");
    return -999;
  }
  
  Int_t n_found_peaks = fCurSpec->GetNPeaks();
  
  if(n_found_peaks==0){
    Info("FitSinglePeak","No peak found in current spectrum.");
    return -999;
  }

  if(peak_num > n_found_peaks){
    gROOT->Warning("FitSinglePeak","Peak number %d not available.",peak_num);
    return -999;
  }
  
  // note: peaks are evaluated from right to left by TSpectrum.
  Float_t peak_pos = (*(fCurSpec->GetPositionX() + peak_num - 1));
  Float_t peak_height = (*(fCurSpec->GetPositionY() + peak_num - 1 ));
  

  //  Printf("Now fitting peak number %d",peak_num);
  Double_t FitRange[2] = {0};
  FitRange[0] = peak_pos - (3./fResolutionFactor)*sigma_init;
  FitRange[1] = peak_pos + (3./fResolutionFactor)*sigma_init;

  TString curpknam = "FndPeak_GausYield"; curpknam+=peak_num;
  TF1 *cur_peak = new TF1(curpknam,FitFCN_GausYield,FitRange[0],FitRange[1],3);

  SetFitFCN_ParNames_GausYield(cur_peak);

  cur_peak->SetParameter(0,peak_height * (sigma_init/peak_pos));
  cur_peak->SetParameter(1,peak_pos);
  cur_peak->SetParameter(2,sigma_init);

  cur_peak->SetParLimits(1,peak_pos - sigma_init,peak_pos + sigma_init);
  cur_peak->SetParLimits(2,
                         (sigma_init > sigma_maxerr) ? (sigma_init-sigma_maxerr) : 1e-10,
                         sigma_init+sigma_maxerr);
  
  Int_t fit_res = fCurHist->Fit(curpknam,"RQ+");
  if( fit_res != 0 ) Printf("Fit error: %d",fit_res);
  
  delete cur_peak;
  cur_peak = 0;

  return fit_res;
}

//_____________________
void TFndFitter::PrintFittedPeaks_GausYield(){
  
  if(!fCurHist){
    Warning("PrintFittedPeaks_GausYield","No histogram selected: nothing will happen...");
    return;
  }

  TIter next(fCurHist->GetListOfFunctions());
  TObject *obj;
  
  while ((obj = next())) {
    TString objnam = obj->GetName();
    if( objnam.BeginsWith("FndPeak_GausYield") == kFALSE) continue;
  
    TF1 *cur_peak = (TF1 *)obj;
    // ---
    TString resmsg = "";
    TString resline = "";
    resline.Form("\nFit Results (peak name: \"%s\"):\n",cur_peak->GetName());
    resmsg+=resline;
    
    for(Int_t i=0;i<3;i++){
      resline.Form("\t-> %s:\t%.4f\n",cur_peak->GetParName(i),cur_peak->GetParameter(i));
      resmsg+=resline;
    }
    Double_t Area = cur_peak->GetParameter(0) / fCurHist->GetBinWidth(1);
    Double_t Area_68 = Area * 0.683;
    Double_t Area_95 = Area * 0.9545;
    Double_t Area_99 = Area * 0.997;
    resline.Form("\t-> area:\t%.4f  (total)\n",Area);  resmsg+=resline;
    resline.Form("\t-> area:\t%.4f  (C.L: 68% [+- 1 sigma])\n",Area_68);  resmsg+=resline;
    resline.Form("\t-> area:\t%.4f  (C.L: 95% [+- 2 sigma])\n",Area_95);  resmsg+=resline;
    resline.Form("\t-> area:\t%.4f  (C.L: 99% [+- 3 sigma])\n",Area_99);  resmsg+=resline;

    Printf(resmsg.Data());
    // ---
  }
  
}

//_____________________
TF1 *TFndFitter::LocateFittedPeak_GausYield(Double_t low_lim,Double_t high_lim){
  // 

  if(!fCurHist){
    Warning("LocateFittedPeaks_GausYield","No histogram selected: nothing will happen...");
    return 0;
  }

  TIter next(fCurHist->GetListOfFunctions());
  TObject *obj;
 
  TF1 *cur_peak = 0;
  while ((obj = next())) {
    TString objnam = obj->GetName();
    if( objnam.BeginsWith("FndPeak_GausYield") == kFALSE) continue;
    
    cur_peak = (TF1 *)obj;
    Double_t cur_avg = cur_peak->GetParameter(1);
    if(cur_avg >= low_lim && cur_avg <= high_lim){
      //Printf("LocateFittedPeak_GausYield: peak found at pos: %.4f ",cur_avg);
      return cur_peak;
    }    
  }

  Printf("LocateFittedPeak_GausYield: peak not found in range [%.4f ; %.4f]",low_lim,high_lim);
  return 0;
}

---