test/BSMTriggerCheck/CompHisto1D.cc

0001 #include <TH1D.h>
0002 #include <TCanvas.h>
0003 #include <TLegend.h>
0004 #include "CompHisto1D.hh"
0005
0006 CompHisto1D::CompHisto1D(TH1D* histo1_v, TH1D* histo2_v){
0007
0008   histo1 = histo1_v;
0009   histo2 = histo2_v;
0010   label1 = "first set";
0011   label1 = "second set";
0012 }
0013
0014 double CompHisto1D::Compare() {
0015   title = string(histo1->GetName());
0016   if(title != string(histo2->GetName()) ){
0017     cout << "CompHisto1D::Compare WARNING: histograms have different" << endl;
0018     cout << "                              labels. Using the label of first plot" << endl;
0019   }
0020
0021   compCanvas = new TCanvas(title.c_str(),title.c_str(),300,300);
0022
0023   double norm = 1./histo1->Integral();
0024   //  histo1->Scale(norm);
0025   for(int i=1; i<=histo1->GetNbinsX(); i++) {
0026     double content = histo1->GetBinContent(i);
0027     histo1->SetBinContent(i,content*norm);
0028     histo1->SetBinError(i,sqrt(content)*norm);
0029   }
0030
0031   // space for the TLegend
0032   histo1->SetMaximum(1.5*histo1->GetMaximum());
0033
0034   histo1->SetLineColor(2);
0035   //  histo1->SetMarkerStyle(20);
0036   histo1->SetStats(kFALSE);
0037   histo1->Draw("pe");
0038
0039   norm = 1./histo2->Integral();
0040   //  histo2->Scale(norm);
0041   for(int i=1; i<=histo2->GetNbinsX(); i++) {
0042     double content = histo2->GetBinContent(i);
0043     histo2->SetBinContent(i,content*norm);
0044     histo2->SetBinError(i,sqrt(content)*norm);
0045   }
0046
0047   histo2->SetLineColor(4);
0048   //  histo2->SetMarkerStyle(30);
0049   histo2->SetStats(kFALSE);
0050   histo2->Draw("pesame");
0051
0052   TLegend *legend =new TLegend(0.4,0.72,0.95,0.87);
0053   legend->SetTextSize(0.03);
0054   legend->AddEntry(histo1,label1.c_str());
0055   legend->AddEntry(histo2,label2.c_str());
0056   legend->Draw();
0057
0058   // calculate JetMET-like compatibility
0059   double chisq = myChisq();
0060   double kolmo = histo1->KolmogorovTest(histo2);
0061
0062   return min(chisq,kolmo);
0063
0064 }
0065
0066 void CompHisto1D::SaveAsEps(){
0067   compCanvas->SaveAs((title+".eps").c_str());
0068 }
0069
0070 double CompHisto1D::myChisq() {
0071   double chisq = 0.;
0072   int  nDOF = 0;
0073   for(int i=1; i<=histo1->GetXaxis()->GetNbins(); i++) {
0074     if(histo1->GetBinContent(i) != 0. ||
0075        histo2->GetBinContent(i) != 0.) {
0076       chisq +=
0077     std::pow(histo1->GetBinContent(i)-histo2->GetBinContent(i),2.)/
0078     (std::pow(histo1->GetBinError(i),2.) + std::pow(histo2->GetBinError(i),2.));
0079       nDOF += 1;
0080     }
0081   }
0082
0083   return
0084     1./(std::pow(2.,nDOF/2)*TMath::Gamma(nDOF/2))*
0085     std::pow(chisq, nDOF/2-1)*exp(-chisq/2.);
0086 }