/* * * Program : ztest.cpp * Task : Compute the mean, rms, mean error and rms error of a log-equidistant histogram and display the corresponding (corrected) statbox * Compil/Exec : .x ztest.cpp, .x ztest.cpp+ * Author : Mathieu Trocmé (mathieu.trocme@gmail.com) * Date : 4 nov 2012 * Framework : None * */ #include #include #include #include #include #include #include ///// void MRF_histoGetStatsForEquidistantLogBinning ( TH1 * histo, Double_t * stats ) ; /// DESCRIPTION: Computes stats /// USAGE : See functions below Double_t MRF_histoGetMeanForEquidistantLogBinning ( TH1 * histo, Bool_t hasBeenSetBinContentedInsteadOfFilled ) ; Double_t MRF_histoGetRMSForEquidistantLogBinning ( TH1 * histo, Bool_t hasBeenSetBinContentedInsteadOfFilled ) ; Double_t MRF_histoGetMeanErrorForEquidistantLogBinning ( TH1 * histo, Bool_t hasBeenSetBinContentedInsteadOfFilled ) ; Double_t MRF_histoGetRMSErrorForEquidistantLogBinning ( TH1 * histo, Bool_t hasBeenSetBinContentedInsteadOfFilled ) ; /// DESCRIPTION: Returns mean, stddev, mean error and stddev error taking into account log binning if hasBeenSetBinContentedInsteadOfFilled is true i.e. if method TH1::Fill*() has not been used. /// In the latter case stats are computed from exact values (not bin centers) i.e. a specific log treatment does not make any sense. /// No default value is given for the boolean in order people keep this in mind. /// USAGE : Double_t mean = MRF_histoGetMeanForEquidistantLogBinning (histo,kTRUE) ; /// Double_t stddev = MRF_histoGetRMSForEquidistantLogBinning (histo,kTRUE) ; /// Double_t mean_err = MRF_histoGetMeanErrorForEquidistantLogBinning (histo,kTRUE), /// Double_t stddev_err = MRF_histoGetRMSErrorForEquidistantLogBinning (histo,kTRUE), TPaveStats * MRF_createStatBox ( TH1 * histo, Float_t x1, Float_t y1, Float_t x2, Float_t y2, Int_t color=-1, Int_t style=-1, Int_t width=-1, Float_t textSize=0.04, const TString & title_="", Bool_t displayErrors=kTRUE, Bool_t displayUnderFlow=kTRUE, Bool_t displayOverFlow=kTRUE, Bool_t displayIntegral=kTRUE, Bool_t displayEntries=kTRUE, Bool_t fill=kTRUE, Bool_t histoLogEquidistantAndNotFilled=kFALSE ) ; /// DESCRIPTION: Creates the statbox corresponding to histo /// USAGE : MRF_createStatBox (histo, x1,y1,x2,y2, color,-1,-1, 0.04,"", 1,1,1,1,1, 1) ; void MRF_drawStatBox ( TH1 * histo, const TString & name, Float_t x1, Float_t y1, Float_t x2, Float_t y2, Int_t color=-1, Int_t style=-1, Int_t width=-1, Float_t textSize=0.04, const TString & title="", Bool_t displayErrors=kTRUE, Bool_t displayUnderFlow=kTRUE, Bool_t displayOverFlow=kTRUE, Bool_t displayIntegral=kTRUE, Bool_t displayEntries=kTRUE, Bool_t fill=kTRUE, Bool_t histoLogEquidistantAndNotFilled=kFALSE ) ; /// DESCRIPTION: Draws the statbox corresponding to histo /// USAGE : MRF_drawStatBox (histo,"", x1,y1,x2,y2, color,-1,-1, 0.04,"", 1,1,1,1,1, 1) ; ///// void MRF_histoGetStatsForEquidistantLogBinning ( TH1 * histo, Double_t * stats ) { /// See TH1::GetStats() /// For "setBinContented" histos, one cannot take into account overflows since the x information is lost! Double_t x=-1., w=-1., err=-1. ; for (Int_t k=0; k<4; k++) { stats[k] = 0. ; } for (Int_t ibin=1; ibin<=histo->GetNbinsX(); ibin++) { //x = histo->GetXaxis()->GetBinCenter (ibin) ; /// debug x = histo->GetXaxis()->GetBinCenterLog(ibin) ; w = histo->GetBinContent(ibin) ; err = TMath::Abs(histo->GetBinError(ibin)) ; stats[0] += w ; stats[1] += err*err ; stats[2] += w*x ; stats[3] += w*x*x ; } return ; } Double_t MRF_histoGetMeanForEquidistantLogBinning ( TH1 * histo, Bool_t hasBeenSetBinContentedInsteadOfFilled ) { /// See TH1::GetMean() Double_t stats[4] ; histo->GetStats(stats) ; if ( stats[0] == 0. ) { return (0.) ; } if ( hasBeenSetBinContentedInsteadOfFilled ) { MRF_histoGetStatsForEquidistantLogBinning (histo,stats) ; } return ( stats[2]/stats[0] ) ; } Double_t MRF_histoGetRMSForEquidistantLogBinning ( TH1 * histo, Bool_t hasBeenSetBinContentedInsteadOfFilled ) { /// See TH1::GetRMS() Double_t stats[4] ; histo->GetStats(stats) ; if ( stats[0] == 0. ) { return (0.) ; } if ( hasBeenSetBinContentedInsteadOfFilled ) { MRF_histoGetStatsForEquidistantLogBinning (histo,stats) ; } Double_t mean = stats[2]/stats[0] ; Double_t rms2 = TMath::Abs( stats[3]/stats[0] - mean*mean ) ; return ( TMath::Sqrt(rms2) ) ; } Double_t MRF_histoGetMeanErrorForEquidistantLogBinning ( TH1 * histo, Bool_t hasBeenSetBinContentedInsteadOfFilled ) { /// See TH1::GetMeanError(): mean error = RMS / sqrt( Neff ) Double_t rms=-1., neff=-1., mean_err=-1. ; rms = MRF_histoGetRMSForEquidistantLogBinning (histo,hasBeenSetBinContentedInsteadOfFilled) ; neff = histo->GetEffectiveEntries() ; if ( neff > 0. ) { mean_err = rms / TMath::Sqrt(neff) ; } else { mean_err = 0. ; } return ( mean_err ) ; } Double_t MRF_histoGetRMSErrorForEquidistantLogBinning ( TH1 * histo, Bool_t hasBeenSetBinContentedInsteadOfFilled ) { /// See TH1::GetRMSError(): The right formula for RMS error depends on 4th momentum (see Kendall-Stuart Vol 1 pag 243). The formula below is valid for only gaussian distribution ( 4-th momentum = 3 * sigma^4 ) Double_t rms=-1., neff=-1., rms_err=-1. ; rms = MRF_histoGetRMSForEquidistantLogBinning (histo,hasBeenSetBinContentedInsteadOfFilled) ; neff = histo->GetEffectiveEntries() ; if ( neff > 0. ) { rms_err = TMath::Sqrt((rms*rms)/(2.*neff)) ; } else { rms_err = 0. ; } return ( rms_err ) ; } /// TPaveStats * MRF_createStatBox ( TH1 * histo, Float_t x1, Float_t y1, Float_t x2, Float_t y2, Int_t color, Int_t style, Int_t width, Float_t textSize, const TString & title_, Bool_t displayErrors, Bool_t displayUnderFlow, Bool_t displayOverFlow, Bool_t displayIntegral, Bool_t displayEntries, Bool_t fill, Bool_t histoLogEquidistantAndNotFilled ) { if ( x2 <= x1 ) { gSystem->Error("MRF_createStatBox",TString::Format("x2=%e <= x1=%e!",x2,x1)) ; gSystem->Abort() ; } if ( y2 <= y1 ) { gSystem->Error("MRF_createStatBox",TString::Format("y2=%e <= y1=%e!",y2,y1)) ; gSystem->Abort() ; } TPaveStats * statbox = new TPaveStats (x1,y1,x2,y2,"NDC") ; statbox->SetBorderSize(2) ; statbox->SetTextAlign(12) ; statbox->SetTextFont(42) ; statbox->SetTextSize(textSize) ; statbox->SetShadowColor(kWhite) ; if ( color == -1 ) { color = histo->GetLineColor() ; } statbox->SetLineColor(color) ; statbox->SetTextColor(color) ; if ( style == -1 ) { style = histo->GetLineStyle() ; } statbox->SetLineStyle(style) ; if ( width == -1 ) { width = histo->GetLineWidth() ; } statbox->SetLineWidth(width) ; if ( fill ) { statbox->SetFillColor(kWhite) ; } else { statbox->SetFillStyle(0) ; } TString title = title_ ; if ( title == "-1" ) { title = histo->GetName() ; } if ( title != "" ) { statbox->AddText(title) ; statbox->InsertLine() ; } if ( displayEntries ) { statbox->AddText(TString::Format("Entries = %d",Int_t(histo->GetEntries()))) ; } statbox->AddText(TString(Form("Mean = %s%s", TString(Form(TString(Form("%%%s",gStyle->GetStatFormat())).Data(), histoLogEquidistantAndNotFilled?MRF_histoGetMeanForEquidistantLogBinning(histo,kTRUE):histo->GetMean())).Data(), displayErrors?TString(Form(" #pm %s",TString(Form(TString(Form("%%%s",gStyle->GetStatFormat())).Data(),histoLogEquidistantAndNotFilled?MRF_histoGetMeanErrorForEquidistantLogBinning(histo,kTRUE):histo->GetMeanError())).Data())).Data():"" ))) ; statbox->AddText(TString(Form("RMS = %s%s", TString(Form(TString(Form("%%%s",gStyle->GetStatFormat())).Data(), histoLogEquidistantAndNotFilled?MRF_histoGetRMSForEquidistantLogBinning (histo,kTRUE):histo->GetRMS ())).Data(), displayErrors?TString(Form(" #pm %s",TString(Form(TString(Form("%%%s",gStyle->GetStatFormat())).Data(),histoLogEquidistantAndNotFilled?MRF_histoGetRMSErrorForEquidistantLogBinning (histo,kTRUE):histo->GetRMSError ())).Data())).Data():"" ))) ; if ( displayUnderFlow ) { statbox->AddText(TString::Format("Underflow = %d",Int_t(histo->GetBinContent(0)))) ; } if ( displayOverFlow ) { statbox->AddText(TString::Format("Overflow = %d",Int_t(histo->GetBinContent(histo->GetNbinsX()+1)))) ;} if ( displayIntegral ) { statbox->AddText(TString::Format("Integral = %s",TString::Format(TString::Format("%%%s",gStyle->GetStatFormat()).Data(),histo->Integral()).Data())) ; } return (statbox) ; } void MRF_drawStatBox ( TH1 * histo, const TString & name, Float_t x1, Float_t y1, Float_t x2, Float_t y2, Int_t color, Int_t style, Int_t width, Float_t textSize, const TString & title, Bool_t displayErrors, Bool_t displayUnderFlow, Bool_t displayOverFlow, Bool_t displayIntegral, Bool_t displayEntries, Bool_t fill, Bool_t histoLogEquidistantAndNotFilled ) { TPaveStats * statbox = MRF_createStatBox (histo,x1,y1,x2,y2,color,style,width,textSize,title,displayErrors,displayUnderFlow,displayOverFlow,displayIntegral,displayEntries,fill,histoLogEquidistantAndNotFilled) ; statbox->SetName(name) ; statbox->Draw() ; return ; } ///// void ztest () { gStyle->SetStatFormat(".5f") ; gStyle->SetOptStat("neMRoui") ; //TH1::StatOverflows(kTRUE) ; Double_t xmin=0.1, xmax=10. ; Int_t npoints=50 ; Double_t logxmin=TMath::Log10(xmin), logxmax=TMath::Log10(xmax) ; Int_t nbins = npoints-1 ; Double_t binwidth=(logxmax-logxmin)/(nbins*1.) ; Double_t * xbins = new Double_t [npoints] ; xbins[0] = xmin ; for (Int_t i=1; iSetBins(nbins,xbins) ; TH1D * h2 = new TH1D ("h2_setBinContent_changed","h2",nbins,xmin,xmax) ; h2->SetBins(nbins,xbins) ; h1->FillRandom("gaus",10000) ; //h1->Fill(0.01) ; h1->Fill(8756.) ; h1->Fill(2312.) ; for (Int_t ibin=1; ibin<=h2->GetNbinsX(); ibin++) { h2->SetBinContent(ibin,h1->GetBinContent(ibin)/1000.) ; } /// Any histogram read from a file //printf ("\nh1->GetEffectiveEntries() = %g",h1->GetEffectiveEntries()) ; fflush(stdout) ; //printf ("\nh2->GetEffectiveEntries() = %g",h2->GetEffectiveEntries()) ; fflush(stdout) ; TCanvas * canvas = new TCanvas ("canvas") ; canvas->Divide(1,2) ; canvas->cd(1) ; gPad->SetLogx() ; h1->Draw() ; MRF_drawStatBox (h1, "", 0.70,0.18,0.92,0.50, kRed,-1,-1, 0.04,"-1", 1,1,1,1,1, 1, 0) ; canvas->cd(2) ; gPad->SetLogx() ; h2->Draw() ; MRF_drawStatBox (h2, "", 0.70,0.18,0.92,0.50, kRed,-1,-1, 0.04,"-1", 1,1,1,1,1, 1, 1) ; Double_t mean_h1_mt = MRF_histoGetMeanForEquidistantLogBinning (h1,kFALSE), mean_h1_rt = h1->GetMean(), mean_err_h1_mt = MRF_histoGetMeanErrorForEquidistantLogBinning (h1,kFALSE), mean_err_h1_rt = h1->GetMeanError() ; Double_t mean_h2_mt = MRF_histoGetMeanForEquidistantLogBinning (h2,kTRUE ), mean_h2_rt = h2->GetMean(), mean_err_h2_mt = MRF_histoGetMeanErrorForEquidistantLogBinning (h2,kTRUE ), mean_err_h2_rt = h2->GetMeanError() ; Double_t rms_h1_mt = MRF_histoGetRMSForEquidistantLogBinning (h1,kFALSE), rms_h1_rt = h1->GetRMS (), rms_err_h1_mt = MRF_histoGetRMSErrorForEquidistantLogBinning (h1,kFALSE), rms_err_h1_rt = h1->GetRMSError () ; Double_t rms_h2_mt = MRF_histoGetRMSForEquidistantLogBinning (h2,kTRUE ), rms_h2_rt = h2->GetRMS (), rms_err_h2_mt = MRF_histoGetRMSErrorForEquidistantLogBinning (h2,kTRUE ), rms_err_h2_rt = h2->GetRMSError () ; printf ("\n") ; printf ("\nh1(fill): mean MT/RT = %.8f/%.8f = %.8f, mean_err MT/RT = %.8f/%.8f = %.8f", mean_h1_mt,mean_h1_rt,mean_h1_mt/mean_h1_rt, mean_err_h1_mt,mean_err_h1_rt,mean_err_h1_mt/mean_err_h1_rt) ; fflush(stdout) ; printf ("\nh2(setc): mean MT/RT = %.8f/%.8f = %.8f, mean_err MT/RT = %.8f/%.8f = %.8f", mean_h2_mt,mean_h2_rt,mean_h2_mt/mean_h2_rt, mean_err_h2_mt,mean_err_h2_rt,mean_err_h2_mt/mean_err_h2_rt) ; fflush(stdout) ; printf ("\n") ; printf ("\nh1(fill): rms MT/RT = %.8f/%.8f = %.8f, rms_err MT/RT = %.8f/%.8f = %.8f", rms_h1_mt ,rms_h1_rt ,rms_h1_mt /rms_h1_rt , rms_err_h1_mt ,rms_err_h1_rt ,rms_err_h1_mt /rms_err_h1_rt ) ; fflush(stdout) ; printf ("\nh2(setc): rms MT/RT = %.8f/%.8f = %.8f, rms_err MT/RT = %.8f/%.8f = %.8f", rms_h2_mt ,rms_h2_rt ,rms_h2_mt /rms_h2_rt , rms_err_h2_mt ,rms_err_h2_rt ,rms_err_h2_mt /rms_err_h2_rt ) ; fflush(stdout) ; printf ("\n\n") ; return ; }