DoFit1.C

Go to the documentation of this file.

void DoFit1(Float_t sigma = 1.0, Char_t *logname = NULL)
{
   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];
   
   FILE *flog; if ( TVirtualPad::Pad() == NULL ) return;

   TH1 *hist;
   TIter next(gPad->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; }
        
   if ( logname ) { // results are written in a file
        flog = fopen(logname,"a"); if ( flog == NULL ) 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;
     }
     // 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;

     nbgauss++;
   }
   
   // perform the final fit and display the result
   TF1 *f; if ( ( f = (TF1 *)gROOT->GetListOfFunctions()->FindObject("ResultFit") ) ) { delete f; }
   
   ffit.Remove(0,1); 
   f = new TF1("ResultFit",ffit.Data(),xlow,xhigh); 

   // init parameters and limits 
   for ( i = 0; i < nbgauss; i++ ) {
        f->SetParameter(i,par[i]);
        // 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,"NR");

   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; }
     
     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 ( logname ) { // results are kept in a file
        fprintf(flog,"#");
        fprintf(flog,"# Position %.2f +- %.3f, width %.3f +- %.3f, Height %.1f +- %.1f, \t Intensity %.1f +- %.1f");
        fprintf(flog,"#");      
        fprintf(flog,"%.2f\t %.3f\t%.3f\t%.3f\t%.1f\t%.1f\t%.1f\t%.1f\n",
            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);     
     }   
   }
   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"); if ( logname ) { fclose(flog); }  TVirtualPad::Pad()->Update(); 
}