PeakFitter.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004-2006 by Olivier Stezowski & Christophe Theisen     *
 *   stezow(AT)ipnl.in2p3.fr, christophe.theisen(AT)cea.fr                       *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
 
#ifndef GW_PEAKFITTER_H
#include "PeakFitter.h"
#endif
#ifndef GW_CLICKCOLLECTOR_H
#include "ClickCollector.h"
#endif
#ifndef ROOT_TMath_H
#include "TMath.h"
#endif
#ifndef ROOT_TAxis_H
#include "TAxis.h"
#endif
#ifndef ROOT_TLine_H
#include "TLine.h"
#endif
#ifndef ROOT_TF1_H
#include "TF1.h"
#endif
#ifndef ROOT_TCanvas_H
#include "TCanvas.h"
#endif
#ifndef ROOT_TVirtualFitter_H
#include "TVirtualFitter.h"
#endif
#ifndef ROOT_TObject_H
#include "TObject.h"
#endif
#ifndef ROOT_TQObject_H
#include "TQObject.h"
#endif


#include <iostream>
#include <fstream>

using namespace std;
using namespace GwD;
 
// ROOT dictionnary
ClassImp(PeakFitter);

Double_t SkewedGaussian(Double_t *x, Double_t *par){  //SkewedGaussian function (See Radford WebSite for details)
        return(par[0]*TMath::Exp(-0.5*TMath::Power(((x[0]-par[1])/par[2]),2)) +
               par[3]*TMath::Exp((x[0]-par[1])/par[4])*TMath::Erfc(((x[0]-par[1])/(TMath::Sqrt(2.0)*par[2])) +
               (par[2]/(TMath::Sqrt(2.0)*par[4]))) + par[5]*TMath::Erfc((x[0]-par[1])/(TMath::Sqrt(2.0)*par[2])));
}

PeakFitter::PeakFitter() : fPeaks(0), fType(kGauss)//Defaut Constructor
{
        //fGaussian = new TF1("Gaussian",Gaussian);
        //fGaussian->SetParNames("Constant","Mean","Sigma");
        //fGaussian->SetLineColor(2);
        //fSkewed = new TF1("Skewed",SkewedGaussian);
        //fSkewed->SetParNames("ConstantGaus","Mean","Sigma","ConstantSkewed","Beta","ConstantStep");
        //fSkewed->SetLineColor(3);

        
         fMainFrame = new TGMainFrame(0,10,10,kMainFrame | kVerticalFrame);
   fMainFrame->SetLayoutBroken(kTRUE);
   
   fGroupFrame1 = new TGGroupFrame(fMainFrame,"PeakFitter");
   fGroupFrame1->SetLayoutBroken(kTRUE);

         // button "CollectOn"
   fTextButtonCollectOn = new TGTextButton(fGroupFrame1,"CollectOn");
   fTextButtonCollectOn->SetTextJustify(36);
   fGroupFrame1->AddFrame(fTextButtonCollectOn, new TGLayoutHints(kLHintsLeft | kLHintsTop,2,2,2,2));
   fTextButtonCollectOn->MoveResize(10,30,100,40);

         // button "CollectOff"
   fTextButtonCollectOff = new TGTextButton(fGroupFrame1,"CollectOff");
   fTextButtonCollectOff->SetTextJustify(36);
   fGroupFrame1->AddFrame(fTextButtonCollectOff, new TGLayoutHints(kLHintsLeft | kLHintsTop,2,2,2,2));
   fTextButtonCollectOff->MoveResize(120,30,100,40);

         // button "DoFit"
         fTextButtonDoFit = new TGTextButton(fGroupFrame1,"DoFit");
   fTextButtonDoFit->SetTextJustify(36);
   fGroupFrame1->AddFrame(fTextButtonDoFit, new TGLayoutHints(kLHintsLeft | kLHintsTop,2,2,2,2));
   fTextButtonDoFit->MoveResize(10,80,100,40);

         // Type combo box
   fComboBoxFunctionType = new TGComboBox(fGroupFrame1,"Function Type");
   fComboBoxFunctionType->AddEntry("Gaussian",0);
   fComboBoxFunctionType->AddEntry("Skewed",1);
   fComboBoxFunctionType->Select(0);
   fGroupFrame1->AddFrame(fComboBoxFunctionType, new TGLayoutHints(kLHintsLeft | kLHintsTop,2,2,2,2));
   fComboBoxFunctionType->MoveResize(120,80,100,40);

         // button "ClearMarkers"
   fTextButtonClearMarkers = new TGTextButton(fGroupFrame1,"ClearMarkers");
   fTextButtonClearMarkers->SetTextJustify(36);
   fGroupFrame1->AddFrame(fTextButtonClearMarkers, new TGLayoutHints(kLHintsLeft | kLHintsTop,2,2,2,2));
   fTextButtonClearMarkers->MoveResize(10,130,100,40);

         // button "ClearLastMarker"
   fTextButtonClearLastMarker = new TGTextButton(fGroupFrame1,"ClearLastMarker");
   fTextButtonClearLastMarker->SetTextJustify(36);
   fGroupFrame1->AddFrame(fTextButtonClearLastMarker, new TGLayoutHints(kLHintsLeft | kLHintsTop,2,2,2,2));
   fTextButtonClearLastMarker->MoveResize(120,130,100,40);


   //Gestion des boutons
   fTextButtonCollectOn->Connect("Clicked()","Gw::PeakFitter",this,"CollectOn()");
         fTextButtonCollectOff->Connect("Clicked()","Gw::PeakFitter",this,"CollectOff()");
   fTextButtonDoFit->Connect("Clicked()","Gw::PeakFitter",this,"DoFit1(=1.0)");
   fTextButtonClearMarkers->Connect("Clicked()","Gw::PeakFitter",this,"ClearMarkers()");
   fTextButtonClearLastMarker->Connect("Clicked()","Gw::PeakFitter",this,"ClearLastMarker()");
         fComboBoxFunctionType->Connect("Selected(Int_t)","Gw::PeakFitter",this,"SetFunctionType(Int_t)");
   

   fGroupFrame1->SetLayoutManager(new TGVerticalLayout(fGroupFrame1));
   fMainFrame->AddFrame(fGroupFrame1, new TGLayoutHints(kLHintsLeft | kLHintsTop,2,2,2,2));
   fGroupFrame1->MoveResize(10,10,250,180);

   fMainFrame->SetMWMHints(kMWMDecorAll,kMWMFuncAll,kMWMInputModeless);
   fMainFrame->MapSubwindows();
   fMainFrame->MapWindow();
   fMainFrame->Resize(270,200);
}

PeakFitter::~PeakFitter(){  //Destructor
        delete fMainFrame;
        delete fGroupFrame1;
        delete fTextButtonCollectOn;
        delete fTextButtonDoFit;
        delete fTextButtonClearMarkers;
        delete fTextButtonClearLastMarker;
        delete fComboBoxFunctionType;
}

void PeakFitter::CollectOn()
{
        if ( TVirtualPad::Pad() ) {
                TCanvas *c = TVirtualPad::Pad()->GetCanvas(); ClickCollector *thecollector = ClickCollector::TheCollector();
                c->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)","Gw::ClickCollector",thecollector,"DoCollect(Int_t,Int_t,Int_t,TObject*)");
                cout << " - Start collecting clicks in " << c->GetName() << endl ;      
                
                thecollector->SetMode();
        }
        else 
                cout << " - Cannot collect clicks, there is no current pad !!! " << endl ;       

}

void PeakFitter::CollectOff()
{
        if ( TVirtualPad::Pad() ) {
          TCanvas *c = TVirtualPad::Pad()->GetCanvas(); ClickCollector *thecollector = ClickCollector::TheCollector();
          thecollector->SetMode(false);
          cout << " - Stop collecting clicks in " << c->GetName() << endl ;
        }
}


void PeakFitter::ClearMarkers()
{
        TObject *obj; if ( TVirtualPad::Pad() == NULL ) return;
   
        TIter next(TVirtualPad::Pad()->GetListOfPrimitives());
        while ( (obj = next()) ) { // look for markers
                if ( obj->InheritsFrom("TMarker") ) { delete obj; }
        }
        TVirtualPad::Pad()->Update();
}

void PeakFitter::ClearLastMarker()
{
        TObject *obj; if ( TVirtualPad::Pad() == NULL ) return;
   
        TIter next(TVirtualPad::Pad()->GetListOfPrimitives(),kIterBackward);
        while ( (obj = next()) ) { // look for markers
                if ( obj->InheritsFrom("TMarker") ) { delete obj; break; }
        }
        TVirtualPad::Pad()->Update();
}


void PeakFitter::DoFit1(Float_t sigma)
{
   Float_t beta = 5.;
   const Int_t MAXMARKERS = 50 ; TVirtualFitter::Fitter(0,3*MAXMARKERS); //to get a maximum of 150 parameters
  
   Int_t i, nbmarkers, nbgauss, ind[MAXMARKERS]; 
   Float_t x[MAXMARKERS], y[MAXMARKERS], /*xgauss[MAXMARKERS], ygauss[MAXMARKERS],*/ xlow, xhigh; Double_t par[3*MAXMARKERS];

         TObject *obj;
   
   if ( TVirtualPad::Pad() == NULL ) return;

   TH1 *hist;
   TIter next(TVirtualPad::Pad()->GetListOfPrimitives());

   hist = NULL; nbmarkers = 0;
   while ( (obj = next()) ) { // look for an histogram and for markers
         if ( obj->InheritsFrom("TMarker") ) {
                x[nbmarkers]   = ((TMarker *)obj)->GetX();
                y[nbmarkers++] = ((TMarker *)obj)->GetY();
         }
         if ( obj->InheritsFrom("TH1") && !obj->InheritsFrom("TH2") ) { hist = (TH1 *)obj; }
   }
   if ( hist == NULL || nbmarkers < 3 ) {
        printf("There is no histograms or enough clicks in the current pad to realise a graphical fit \n"); return; }
        

   // first one and last one are for background
   TMath::Sort(nbmarkers,x,ind,false); xlow = x[ind[0]]; xhigh = x[ind[nbmarkers-1]];

   // histograms used to fit
   TH1 *lochist = (TH1 *)hist->Clone("fitter"); lochist->Reset("ICE");

   // to get the backgound
   Int_t first_bin, last_bin; Stat_t first_bin_c, last_bin_c, bin_c;

   first_bin = hist->FindBin(  xlow ); 
   last_bin  = hist->FindBin( xhigh );
   first_bin_c = hist->GetBinContent(first_bin); last_bin_c = hist->GetBinContent(last_bin);

   // determine the linear function for the background
   TF1 *fbg;
   if ( ( fbg = (TF1 *)gROOT->GetListOfFunctions()->FindObject("Background") ) ) { delete fbg; }

   fbg = new TF1("Background","pol1",xlow,xhigh);
   fbg->SetParameters(first_bin_c-xlow*(last_bin_c-first_bin_c)/(xhigh-xlow),(last_bin_c-first_bin_c)/(xhigh-xlow));

   fbg->SetLineColor(30); fbg->SetLineWidth(3); fbg->SetLineStyle(1); fbg->Draw("same");

   // fill the histogram that is used to fit
   for ( i = first_bin; i < last_bin ; i++ ) {
            bin_c = hist->GetBinContent(i) -            
                    ( Float_t(i-first_bin)*(last_bin_c-first_bin_c)/(last_bin-first_bin) + first_bin_c ) ;
            lochist->SetBinContent(i,bin_c);
   }
   // now the function used to fit is created
   char fname[500]; TString ffit, fdisplay = "pol1(0)";

   nbgauss = 0;
   for ( i = 1; i < nbmarkers - 1; i++ ) {

     if ( y[ind[i]] < fbg->Eval(x[ind[i]]) ) { // background too high cannot add a gaussian at this position
       printf("background too high at position %.1f \n",x[ind[i]]); continue;
     }

                 if (fType == kGauss) {
        // add a new gaussian and init parameters
         par[3*nbgauss]   = y[ind[i]] - fbg->Eval(x[ind[i]]);
     par[3*nbgauss+1] = x[ind[i]];
     par[3*nbgauss+2] = sigma;
     sprintf(fname,"gaus(%d)",nbgauss*3);   ffit += "+"; ffit += fname;
     sprintf(fname,"gaus(%d)",2+nbgauss*3); fdisplay += "+"; fdisplay += fname;
                 }

                 if (fType == kSkewed) {
        // add a new SkewedGaussian and init parameters
         par[6*nbgauss]   = y[ind[i]] - fbg->Eval(x[ind[i]]);
     par[6*nbgauss+1] = x[ind[i]];
     par[6*nbgauss+2] = sigma;
                 par[6*nbgauss+3] = (y[ind[i]] - fbg->Eval(x[ind[i]]))/10.;
                 par[6*nbgauss+4] = beta;
                 par[6*nbgauss+5] = 0.;
                 
                        sprintf(fname, "[%d]*TMath::Exp(-0.5*TMath::Power(((x-[%d])/[%d]),2)) + [%d]*TMath::Exp((x-[%d])/[%d])*TMath::Erfc(((x-[%d])/(TMath::Sqrt(2.0)*[%d])) + ([%d]/(TMath::Sqrt(2.0)*[%d]))) + [%d]*TMath::Erfc((x-[%d])/(TMath::Sqrt(2.0)*[%d]))",6*nbgauss,6*nbgauss+1,6*nbgauss+2,6*nbgauss+3,6*nbgauss+1,6*nbgauss+4,6*nbgauss+1,6*nbgauss+2,6*nbgauss+2,6*nbgauss+4,6*nbgauss+5,6*nbgauss+1,6*nbgauss+2);

                        ffit += "+"; ffit += fname;

                        sprintf(fname, "[%d]*TMath::Exp(-0.5*TMath::Power(((x-[%d])/[%d]),2)) + [%d]*TMath::Exp((x-[%d])/[%d])*TMath::Erfc(((x-[%d])/(TMath::Sqrt(2.0)*[%d])) + ([%d]/(TMath::Sqrt(2.0)*[%d]))) + [%d]*TMath::Erfc((x-[%d])/(TMath::Sqrt(2.0)*[%d]))",2+6*nbgauss,2+6*nbgauss+1,2+6*nbgauss+2,2+6*nbgauss+3,2+6*nbgauss+1,2+6*nbgauss+4,2+6*nbgauss+1,2+6*nbgauss+2,2+6*nbgauss+2,2+6*nbgauss+4,2+6*nbgauss+5,2+6*nbgauss+1,2+6*nbgauss+2);
      fdisplay += "+"; fdisplay += fname;
                 }
                 
     nbgauss++;
   }
   
   // perform the final fit and display the result
   TF1 *f; if ( ( f = (TF1 *)gROOT->GetListOfFunctions()->FindObject("ResultFit") ) ) { delete f; }
   
   ffit.Remove(0,1);
   printf("%s \n",ffit.Data());
   f = new TF1("ResultFit",ffit.Data(),xlow,xhigh);

   // init parameters and limits 
   for ( i = 0; i < nbgauss; i++ ) {
           if (fType == kGauss) {f->SetParameter(3*i,par[3*i]); f->SetParameter(3*i+1,par[3*i+1]); f->SetParameter(3*i+2,par[3*i+2]);}
                 if (fType == kSkewed) {f->SetParameter(6*i,par[6*i]); f->SetParameter(6*i+1,par[6*i+1]); f->SetParameter(6*i+2,par[6*i+2]);  f->SetParameter(6*i+3,par[6*i+3]); f->SetParameter(6*i+4,par[6*i+4]); f->FixParameter(6*i+5,par[6*i+5]);}
        // futur developments in case we would like to add limits
        // f->SetParLimits(3*i,);
        // f->SetParLimits(3*i+1,);
        // f->SetParLimits(3*i+2,0.5*sigma,1.5*sigma);  
        //f->ReleaseParameter(3*i);f->ReleaseParameter(3*i+1);f->ReleaseParameter(3*i+2);
   }
   lochist->Fit(f,"VNR");

   Stat_t intensity, intensity_err;
   printf("\n");
   printf("******************* DOFIT1 RESULTS ******************* \n");
   printf("  Chi2 %.2f / %d = %.2f \n",f->GetChisquare(),f->GetNDF(),f->GetChisquare()/f->GetNDF());
  
//   printf("  Background %f %f \n",par[0],par[1]);

   for ( i = 0; i < nbgauss; i++ ) {
   
     intensity = f->GetParameter(2+3*i)*f->GetParameter(3*i)*TMath::Sqrt(2*TMath::Pi());
     intensity_err = intensity * (f->GetParError(1+3*i)/f->GetParameter(1+3*i) + f->GetParError(2+3*i)/f->GetParameter(2+3*i));
     
     printf(" Peak  %d : Position %.2f +- %.3f, width %.3f +- %.3f, Height %.1f +- %.1f, \t Intensity %.1f +- %.1f\n",
            i+1,f->GetParameter(1+3*i),f->GetParError(1+3*i),2*f->GetParameter(2+3*i),2*f->GetParError(2+3*i),
            f->GetParameter(3*i),f->GetParError(3*i),intensity,intensity_err);
   
     sprintf(fname,"peak%d",i);
     TF1 *fpeak;
     if ( ( fpeak = (TF1 *)gROOT->GetListOfFunctions()->FindObject(fname) ) ) { delete fpeak; }

                 if (fType == kGauss) {
        fpeak = new TF1(fname,"gaus",f->GetParameter(1+3*i)-10,f->GetParameter(1+3*i)+10); 
        fpeak->SetLineStyle(3); fpeak->SetLineColor(16);
    
        fpeak->SetParameters(f->GetParameter(3*i),f->GetParameter(1+3*i),f->GetParameter(2+3*i));
                        fpeak->Draw("same");
                 }
                 if (fType == kSkewed) {
        fpeak = new TF1(fname,SkewedGaussian,f->GetParameter(1+6*i)-10,f->GetParameter(1+6*i)+10,6);
        fpeak->SetLineStyle(3); fpeak->SetLineColor(16);
    
        fpeak->SetParameters(f->GetParameter(6*i),f->GetParameter(1+6*i),f->GetParameter(2+6*i),f->GetParameter(3+6*i),f->GetParameter(4+6*i),f->GetParameter(5+6*i));
                        fpeak->Draw("same");
                 }
   }
   printf("****************************************************** \n");

   //  to display the final result
   par[0] = fbg->GetParameter(0);
   par[1] = fbg->GetParameter(1); f->GetParameters(&par[2]);

   if ( ( f = (TF1 *)gROOT->GetListOfFunctions()->FindObject("DisplayFit") ) ) { delete f; }

   f = new TF1("DisplayFit",fdisplay.Data(),xlow,xhigh); f->SetParameters(par);
   f->SetLineColor(6); f->SetLineWidth(3);
   
   delete lochist; // no longer useful

   f->Draw("same");  TVirtualPad::Pad()->Update(); 
}

void PeakFitter::SetFunctionType(const FunctionType type){
        fType = type;
}