File indexing completed on 2024-04-06 12:32:12
0001
0002
0003
0004
0005
0006 #include <iostream.h>
0007
0008 class HistoCompare {
0009
0010 public:
0011
0012 HistoCompare() { std::cout << "Initializing HistoCompare... " << std::endl;printresy=0.96; } ;
0013
0014 void PVCompute(TH1 * oldHisto , TH1 * newHisto , TText * te );
0015 void PVCompute(TH2 * oldHisto , TH2 * newHisto , TText * te );
0016 void PVCompute(TProfile * oldHisto , TProfile * newHisto , TText * te );
0017
0018 private:
0019
0020 Double_t mypv;
0021
0022 TH1 * myoldHisto1;
0023 TH1 * mynewHisto1;
0024
0025 TH2 * myoldHisto2;
0026 TH2 * mynewHisto2;
0027
0028 TProfile * myoldProfile;
0029 TProfile * mynewProfile;
0030
0031 TText * myte;
0032
0033 void printRes(TString theName, Double_t thePV, TText * te);
0034 double printresy;
0035 };
0036
0037 HistoCompare::printRes(TString theName, Double_t thePV, TText * te)
0038 {
0039 myte->DrawTextNDC(0.1,printresy,theName );
0040 std::cout << "[Compatibility test] " << theName << std::endl;
0041 printresy+=-0.04;
0042 }
0043
0044
0045 HistoCompare::PVCompute(TH1 * oldHisto , TH1 * newHisto , TText * te )
0046 {
0047
0048 myoldHisto1 = oldHisto;
0049 mynewHisto1 = newHisto;
0050 myte = te;
0051
0052 Double_t *res;
0053
0054 Double_t mypvchi = myoldHisto1->Chi2Test(mynewHisto1,"WW",res);
0055
0056 char str [128];
0057 sprintf(str,"chi^2 P Value: %f",mypvchi);
0058 TString title = str;
0059 printRes(title, mypvchi, myte);
0060 Double_t mypvKS = myoldHisto1->KolmogorovTest(mynewHisto1,"");
0061 sprintf(str,"KS (prob): %f",mypvKS);
0062 TString title = str;
0063 std::strstream buf;
0064 std::string value;
0065
0066 printRes(title, mypvKS, myte);
0067
0068 return;
0069
0070 }
0071
0072 HistoCompare::PVCompute(TH2 * oldHisto , TH2 * newHisto , TText * te )
0073 {
0074
0075 myoldHisto2 = oldHisto;
0076 mynewHisto2 = newHisto;
0077 myte = te;
0078
0079 Double_t *res ;
0080 Double_t mypvchi = myoldHisto2->Chi2Test(mynewHisto2,"WW",res);
0081 char str [128];
0082 sprintf(str,"chi^2 P Value: %f",mypvchi);
0083 TString title = str;
0084 printRes(title, mypvchi, myte);
0085
0086 return;
0087
0088 }
0089
0090
0091 HistoCompare::PVCompute(TProfile * oldHisto , TProfile * newHisto , TText * te )
0092 {
0093
0094 myoldProfile = oldHisto;
0095 mynewProfile = newHisto;
0096 myte = te;
0097
0098 Double_t *res ;
0099
0100 Double_t mypv = myoldProfile->Chi2Test(mynewProfile,"WW",res);
0101 TString title = "chi^2 P Value: ";
0102 printRes(title, mypv, myte);
0103 return;
0104
0105 }
0106
0107 #include "TObject.h"
0108 #include "TDirectory.h"
0109 #include "TKey.h"
0110 #include "TFile.h"
0111 #include "TTree.h"
0112 #include "TText.h"
0113
0114
0115 void DetailedCompare( TString currentfile = "new.root",
0116 TString referencefile = "ref.root",
0117 TString theDir = "DQMData/Run 1/Generator/Run summary/Tau")
0118 {
0119 std::cout << "Note: This code correct the Histograms errors to sqrt(N) - this is not correct for samples with weights" << std::endl;
0120 std::vector<TString> theList = histoList(currentfile, theDir);
0121
0122 gROOT ->Reset();
0123 char* rfilename = referencefile ;
0124 char* sfilename = currentfile ;
0125
0126 delete gROOT->GetListOfFiles()->FindObject(rfilename);
0127 delete gROOT->GetListOfFiles()->FindObject(sfilename);
0128
0129 TFile * rfile = new TFile(rfilename);
0130 TFile * sfile = new TFile(sfilename);
0131
0132 char* baseDir=theDir;
0133
0134 rfile->cd(baseDir);
0135 gDirectory->ls();
0136
0137 sfile->cd(baseDir);
0138 gDirectory->ls();
0139
0140 gStyle->SetLabelFont(63,"X");
0141 gStyle->SetLabelSize(30,"X");
0142 gStyle->SetTitleOffset(1.25,"X");
0143 gStyle->SetTitleFont(63,"X");
0144 gStyle->SetTitleSize(35,"X");
0145 gStyle->SetLabelFont(63,"Y");
0146 gStyle->SetLabelSize(30,"Y");
0147 gStyle->SetTitleOffset(3.0,"Y");
0148 gStyle->SetTitleFont(63,"Y");
0149 gStyle->SetTitleSize(35,"Y");
0150 gStyle->SetLabelFont(63,"Z");
0151 gStyle->SetLabelSize(30,"Z");
0152
0153
0154 for ( unsigned int index = 0; index < theList.size() ; index++ ) {
0155
0156 std::cout << index << std::endl;
0157
0158 TString theName = theDir+"/"+theList[index];
0159 std::cout << theName << std::endl;
0160
0161 TH1* href_;
0162 rfile->GetObject(theName,href_);
0163 href_;
0164 double nentries=href_->GetEntries();
0165 double integral=href_->Integral(0,href_->GetNbinsX()+1);
0166 if(integral!=0)href_->Scale(nentries/integral);
0167 href_->Sumw2();
0168
0169
0170 TH1* hnew_;
0171 sfile->GetObject(theName,hnew_);
0172 hnew_;
0173
0174 double nentries=hnew_->GetEntries();
0175 double integral=hnew_->Integral(0,hnew_->GetNbinsX()+1);
0176 if(integral!=0)hnew_->Scale(nentries/integral);
0177 hnew_->Sumw2();
0178 cout << referencefile << " " << nentries << " " << integral << endl;
0179
0180
0181 DetailedComparePlot(href_, hnew_, currentfile, referencefile, theDir, theList[index]);
0182
0183 }
0184
0185 }
0186
0187 void DetailedComparePlot(TH1 * href_, TH1 * hnew_, TString currentfile, TString referencefile, TString theDir, TString theHisto )
0188 {
0189
0190 gStyle->SetOptTitle(0);
0191
0192 TString theName = theDir+"/"+theHisto;
0193 std::cout << "Histogram name = " << theName << std::endl;
0194
0195 HistoCompare * myPV = new HistoCompare();
0196
0197 int rcolor = 2;
0198 int scolor = 4;
0199
0200 int rmarker = 21;
0201 int smarker = 20;
0202
0203 Double_t markerSize = 0.75;
0204
0205 href_->SetLineColor(rcolor);
0206 href_->SetMarkerStyle(rmarker);
0207 href_->SetMarkerSize(markerSize);
0208 href_->SetMarkerColor(rcolor);
0209
0210 hnew_->SetLineColor(scolor);
0211 hnew_->SetMarkerStyle(smarker);
0212 hnew_->SetMarkerSize(markerSize);
0213 hnew_->SetMarkerColor(scolor);
0214
0215 if ( href_ && hnew_ ) {
0216
0217
0218 TCanvas *myPlot = new TCanvas("myPlot","Histogram comparison",200,10,700,900);
0219 TPad *pad0 = new TPad("pad0","The pad with the function ",0.0,0.9,0.0,1.0);
0220 TPad *pad1 = new TPad("pad1","The pad with the function ",0.0,0.6,0.5,0.9);
0221 TPad *pad2 = new TPad("pad2","The pad with the histogram",0.5,0.6,1.0,0.9);
0222 TPad *pad3 = new TPad("pad3","The pad with the histogram",0.0,0.3,0.5,0.6);
0223 TPad *pad4 = new TPad("pad4","The pad with the histogram",0.5,0.3,1.0,0.6);
0224 TPad *pad5 = new TPad("pad5","The pad with the histogram",0.0,0.0,0.5,0.3);
0225 TPad *pad6 = new TPad("pad6","The pad with the histogram",0.5,0.0,1.0,0.3);
0226
0227 pad1->Draw();
0228 pad2->Draw();
0229 pad3->Draw();
0230 pad4->Draw();
0231 pad5->Draw();
0232 pad6->Draw();
0233
0234
0235
0236 TText titte;
0237 titte.SetTextSize(0.02);
0238 titte.DrawTextNDC(0.1,0.98,theName);
0239 titte.DrawTextNDC(0.1,0.96,"Reference File (A):");
0240 titte.DrawTextNDC(0.1,0.94,referencefile);
0241 titte.DrawTextNDC(0.1,0.92,"Current File (B):");
0242 titte.DrawTextNDC(0.1,0.90,currentfile);
0243
0244 hnew_->SetXTitle(href_->GetTitle());
0245 href_->SetXTitle(href_->GetTitle());
0246 hnew_->SetTitleOffset(1.5,"Y");
0247 href_->SetTitleOffset(1.5,"Y");
0248
0249
0250
0251 pad1->cd();
0252 href_->SetYTitle("Events (A)");
0253 href_->DrawCopy("e1");
0254
0255
0256
0257 pad2->cd();
0258 hnew_->SetYTitle("Events (B)");
0259 hnew_->DrawCopy("e1");
0260
0261 gStyle->SetOptStat("nemruoi");
0262
0263
0264 pad3->cd();
0265
0266 pad3->SetGridx();
0267 pad3->SetGridy();
0268 TH1 *hnew_c=hnew_->Clone();
0269 TH1 *href_c=href_->Clone();
0270 double integral=href_c->Integral(1,href_c->GetNbinsX());
0271 if(integral!=0)href_c->Scale(1/integral);
0272 double integral=hnew_c->Integral(1,href_c->GetNbinsX());
0273 if(integral!=0)hnew_c->Scale(1/integral);
0274 href_c->SetYTitle("Comparison of A and B");
0275 href_c->DrawCopy("e1");
0276 hnew_c->DrawCopy("e1same");
0277 TText* te = new TText();
0278 te->SetTextSize(0.04);
0279 myPV->PVCompute( href_ , hnew_ , te );
0280
0281
0282
0283 TH1 *ratio_=hnew_c->Clone();
0284 ratio_->Divide(href_c);
0285 ratio_->SetYTitle("Ratio: B/A");
0286 pad4->cd();
0287 pad4->SetGridx();
0288 pad4->SetGridy();
0289 ratio_->DrawCopy("e1");
0290
0291 TH1 *diff_=hnew_c->Clone();
0292 diff_->Add(href_c,-1);
0293 diff_->SetYTitle("Difference: B-A");
0294 pad5->cd();
0295 pad5->SetGridx();
0296 pad5->SetGridy();
0297 diff_->DrawCopy("e1");
0298
0299
0300 TH1 *sigma_=diff_->Clone();
0301 for(unsigned int i=1; i<=sigma_->GetNbinsX(); i++ ){
0302 double v=sigma_->GetBinContent(i);
0303 double e=sigma_->GetBinError(i);
0304
0305 if(e!=0)sigma_->SetBinContent(i,v/fabs(e));
0306 else sigma_->SetBinContent(i,0);
0307 }
0308 sigma_->SetYTitle("Sigma on Difference: (B-A)/sigma(B-A)");
0309 pad6->cd();
0310 pad6->SetGridx();
0311 pad6->SetGridy();
0312 sigma_->DrawCopy("phisto");
0313
0314
0315
0316 pad0->cd();
0317
0318 gStyle->SetOptStat(0000000);
0319
0320 }
0321 TString plotFile = theHisto+".eps";
0322 myPlot->Print(plotFile);
0323
0324 delete myPV;
0325 delete myPlot;
0326
0327 }
0328
0329 std::vector<TString> histoList( TString currentfile, TString theDir )
0330 {
0331
0332 gROOT ->Reset();
0333 char* sfilename = currentfile ;
0334
0335 delete gROOT->GetListOfFiles()->FindObject(sfilename);
0336
0337 TFile * sfile = new TFile(sfilename);
0338
0339 char* baseDir=theDir;
0340
0341 sfile->cd(baseDir);
0342
0343 TDirectory * d = gDirectory;
0344
0345 std::vector<TString> theHistList;
0346
0347 TIter i( d->GetListOfKeys() );
0348 TKey *k;
0349 while( (k = (TKey*)i())) {
0350 TClass * c1 = gROOT->GetClass(k->GetClassName());
0351 if ( !c1->InheritsFrom("TH1")) continue;
0352 theHistList.push_back(k->GetName());
0353 }
0354
0355 std::cout << "Histograms considered: " << std::endl;
0356 for (unsigned int index = 0; index < theHistList.size() ; index++ ) {
0357 std::cout << index << " " << theHistList[index] << std::endl;
0358 }
0359
0360 return theHistList;
0361
0362 }