mcr/analysis/PeakFinder.C

//                                             //
// This macro is an example for using the      //
//  most important facilities provided by the  //
//  FIN_AN namespace:                          //
//    (see $FROOTSYS/ANALYSIS/FIN_AN.h)        //
//                                             //
// Peak fitting is also supported by the       //
//  TFndFitter class                           //
//                                             //

#include "Riostream.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TRandom.h"
#include "TThread.h"

#include "FIN_AN.h"

#include "HistogramCreators.C"
using namespace FIN_AN;

Int_t fNFoundPeaks = 0;// evaluated number of peaks (TSpectrum search is performed)

Double_t fResolutionFactor = 1;
// fResolutionFactor see the "FIN_AN" namespace
//
Float_t fPeakPos[FinAn_Nmaxpeaks];
Float_t fPeakHeight[FinAn_Nmaxpeaks];


//_____________________
void FitPeak(TH1F *his,Int_t peak_id,Double_t sigma_init,Double_t sigma_maxerr){
  // this methods uses peaks already evaluated (by FIN_AN::FindPeaks)

  Printf("Now fitting peaks (id:%d)",peak_id);
  Double_t FitRange[2] = {0};
  FitRange[0] = fPeakPos[peak_id] - (3./fResolutionFactor)*sigma_init;
  FitRange[1] = fPeakPos[peak_id] + (3./fResolutionFactor)*sigma_init;

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

  SetFitFCN_ParNames_GausYield(cur_peak);

  cur_peak->SetParameter(0,fPeakHeight[peak_id]*(sigma_init/fPeakPos[peak_id]));
  cur_peak->SetParameter(1,fPeakPos[peak_id]);
  cur_peak->SetParameter(2,sigma_init);

  cur_peak->SetParLimits(1,fPeakPos[peak_id]-sigma_init,fPeakPos[peak_id]+sigma_init);
  cur_peak->SetParLimits(2,
                         (sigma_init > sigma_maxerr) ? (sigma_init-sigma_maxerr) : 1e-10,
                         sigma_init+sigma_maxerr);
  
  Int_t fit_res = his->Fit(curpknam,"RQ+");
  if( fit_res != 0 ) cout << "Fit error: " << fit_res << endl;
  
  delete cur_peak;
  cur_peak = 0;

}

//_____________________
void FitAllPeaks(TH1F *his,Double_t sigma_start,Double_t sigma_maxerr){
  
  for(Int_t pi=0;pi<fNFoundPeaks;pi++){
    FitPeak(his,pi,sigma_start,sigma_maxerr);
  }
  
}

//_____________________
void PrintFittedPeaks(TH1F *his){
  
  TIter next(his->GetListOfFunctions());
  TObject *obj;
  
  TString strchk = "TestPeak";
  while ((obj = next())) {
    TString objnam = obj->GetName();
    objnam.Resize(8);
    if( strchk.CompareTo(objnam) ) continue;
  
    TF1 *cur_peak = (TF1 *)obj;
    // ---
    TString resmsg = "\nFit Results:\n";
    TString 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) / his->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());
    // ---
  }
  
}

//_____________________
void LoadPeaks(TSpectrum *TmpSpectrum){
  
  fNFoundPeaks = TmpSpectrum->GetNPeaks();

  Printf("Found %d candidate peaks to fit\n",fNFoundPeaks);
  for(Int_t pi=0;pi<fNFoundPeaks;pi++){
    // note: peaks are evaluated from right to left by TSpectrum.
    Int_t peakid = fNFoundPeaks - pi -1;
    fPeakPos[pi] = (*(TmpSpectrum->GetPositionX() + peakid));
    fPeakHeight[pi] = (*(TmpSpectrum->GetPositionY() + peakid));
  }
}

//_____________________
void PrintFoundPeaks(){
  TString msg = "\nPeaks position: \n";
  TString linestr ="";
  for(Int_t pi=0;pi<fNFoundPeaks;pi++){
    linestr.Form("Peak number %d: position: %.4f ; height: %.4f \n",
                 pi+1,
                 fPeakPos[pi],
                 fPeakHeight[pi]
                 );
    msg+=linestr;
  }
  Printf(msg.Data());
}


void PeakFinder(){
  gStyle->SetOptStat(1001111);

  // --- peaks simulation
  Int_t n_simul_peaks=7;
  Double_t simul_sigma=0.01;
  Double_t sig_noise_ratio=2.;
     
  fResolutionFactor = 1; // must be >= 1 (1 means 3sigma)
  Bool_t eval_background = kTRUE;
  Int_t n_iter_background = 50;

  // --- peaks evaluation
  Int_t n_max_expected_peaks = 10;
  Double_t expected_sigma = 0.015;
  Double_t sigma_maxerr=0.005;
  Double_t peaks_threshold = 0.3;

  TH1F *start_his = BuildExampleHistogram(n_simul_peaks,simul_sigma,sig_noise_ratio,kTRUE);
  
  //  TH1F *start_his = new TH1F("starthis","starthis",3000,0,1);
  //   for(Int_t i=0;i<20000;i++)  start_his->Fill(gRandom->Gaus(.5,.02) );
  
  
  // --- evaluate background
  TH1F *fBckHis = 0;
  if(eval_background) fBckHis = EvalBackground(start_his,n_max_expected_peaks,n_iter_background);

  // --- check peaks before subtracting background
  LoadPeaks( FindPeaks(start_his,n_max_expected_peaks,expected_sigma,peaks_threshold,fResolutionFactor,"") );
  UpdateDisplay(start_his,fBckHis);
  PrintFoundPeaks();
  //  

  // --- subtract background (if required)
  TH1F *clean_his = (TH1F*)(start_his->Clone("CloneHis") );
  clean_his->SetTitle("Clean histogram (simple clone)");
  if(fBckHis){
    Printf("Now evaluating clean histogram (subtracting background)...");
    clean_his->SetName("CleanHis");
    clean_his->SetTitle("Clean histogram");
    clean_his->Add(fBckHis,-1);
    Printf("...clean histogram evaluated");
  }
  //  

  // --- check peaks after subtracting background
  LoadPeaks( FindPeaks(clean_his,n_max_expected_peaks,expected_sigma,peaks_threshold,fResolutionFactor,"") );
  UpdateDisplay(start_his,fBckHis,clean_his);
  
  // --- perform fit for every peak
  FitAllPeaks(clean_his,expected_sigma,sigma_maxerr);
  Printf("...peak fit completed... :)");
  PrintFittedPeaks(clean_his);
}

---