test/macro/effSigma.C

0001 //******************************************//

0002 // author: Chris Seez

0003 //

0004 //  This function evaluates the effective

0005 //  sigma of a histogram.

0006 //

0007 //

0008 //  TODO list:

0009 //

0010 //

0011 //

0012 //  Modified by:

0013 //

0014 //******************************************//

0015
0016
0017 Double_t effSigma(TH1 * hist)
0018 {
0019
0020   TAxis *xaxis = hist->GetXaxis();
0021   Int_t nb = xaxis->GetNbins();
0022   if(nb < 10) {
0023     cout << "effsigma: Not a valid histo. nbins = " << nb << endl;
0024     return 0.;
0025   }
0026
0027   Double_t bwid = xaxis->GetBinWidth(1);
0028   if(bwid == 0) {
0029     cout << "effsigma: Not a valid histo. bwid = " << bwid << endl;
0030     return 0.;
0031   }
0032   Double_t xmax = xaxis->GetXmax();
0033   Double_t xmin = xaxis->GetXmin();
0034   Double_t ave = hist->GetMean();
0035   Double_t rms = hist->GetRMS();
0036
0037   Double_t total=0.;
0038   for(Int_t i=0; i<nb+2; i++) {
0039     total+=hist->GetBinContent(i);
0040   }
0041   if(total < 100.) {
0042     cout << "effsigma: Too few entries " << total << endl;
0043     return 0.;
0044   }
0045   Int_t ierr=0;
0046   Int_t ismin=999;
0047
0048   Double_t rlim=0.683*total;
0049   Int_t nrms=rms/(bwid);    // Set scan size to +/- rms

0050   if(nrms > nb/10) nrms=nb/10; // Could be tuned...

0051
0052   Double_t widmin=9999999.;
0053   for(Int_t iscan=-nrms;iscan<nrms+1;iscan++) { // Scan window centre

0054     Int_t ibm=(ave-xmin)/bwid+1+iscan;
0055     Double_t x=(ibm-0.5)*bwid+xmin;
0056     Double_t xj=x;
0057     Double_t xk=x;
0058     Int_t jbm=ibm;
0059     Int_t kbm=ibm;
0060     Double_t bin=hist->GetBinContent(ibm);
0061     total=bin;
0062     for(Int_t j=1;j<nb;j++){
0063       if(jbm < nb) {
0064         jbm++;
0065         xj+=bwid;
0066         bin=hist->GetBinContent(jbm);
0067         total+=bin;
0068         if(total > rlim) break;
0069       }
0070       else ierr=1;
0071       if(kbm > 0) {
0072         kbm--;
0073         xk-=bwid;
0074         bin=hist->GetBinContent(kbm);
0075         total+=bin;
0076         if(total > rlim) break;
0077       }
0078       else ierr=1;
0079     }
0080     Double_t dxf=(total-rlim)*bwid/bin;
0081     Double_t wid=(xj-xk+bwid-dxf)*0.5;
0082     if(wid < widmin) {
0083       widmin=wid;
0084       ismin=iscan;
0085     }
0086   }
0087   if(ismin == nrms || ismin == -nrms) ierr=3;
0088   if(ierr != 0) cout << "effsigma: Error of type " << ierr << endl;
0089
0090   return widmin;
0091
0092 }