Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:03:27

0001 #include "TFile.h"
0002 #include "TH1F.h"
0003 #include "TH2F.h"
0004 #include "TKey.h"
0005 #include "TDirectory.h"
0006 #include "TCanvas.h"
0007 #include "TLegend.h"
0008 #include "TObjArray.h"
0009 #include "TObject.h"
0010 #include "TClass.h"
0011 #include <iostream>
0012 #include <string>
0013 #include "TPaveText.h"
0014 
0015 #include <fstream>  // std::ofstream
0016 
0017 //*****************************************************//
0018 Double_t getMaximum(TObjArray *array)
0019 //*****************************************************//
0020 {
0021   Double_t theMaximum = (static_cast<TH1 *>(array->At(0)))->GetMaximum();
0022   for (Int_t i = 0; i < array->GetSize(); i++) {
0023     if ((static_cast<TH1 *>(array->At(i)))->GetMaximum() > theMaximum) {
0024       theMaximum = (static_cast<TH1 *>(array->At(i)))->GetMaximum();
0025       //cout<<"i= "<<i<<" theMaximum="<<theMaximum<<endl;
0026     }
0027   }
0028   return theMaximum;
0029 }
0030 
0031 /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032    MAIN
0033    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
0034 */
0035 
0036 //*****************************************************//
0037 void compareAll(TString file1, TString file2, TString leg1, TString leg2)
0038 //*****************************************************//
0039 {
0040   TFile *f1 = TFile::Open(file1);
0041   TFile *f2 = TFile::Open(file2);
0042 
0043   f1->cd();
0044 
0045   int i = 0;
0046   int j = 0;
0047 
0048   TCanvas dummyC;
0049   dummyC.Print("diff.pdf[");
0050 
0051   std::ofstream ofs("check.txt", std::ofstream::out);
0052 
0053   TIter nextkey(gDirectory->GetListOfKeys());
0054   while (TKey *key = (TKey *)nextkey()) {
0055     //std::cout << "i: " << i << std::endl;
0056     ++i;
0057     TObject *obj = key->ReadObj();
0058     std::string name = obj->GetName();
0059     std::cout << "name: " << name << std::endl;
0060     if (name.find("Residuals") != string::npos) {
0061       std::cout << "skipping:" << name << "folder" << std::endl;
0062       continue;
0063     }
0064     if (obj->IsA()->InheritsFrom("TDirectory")) {
0065       f1->cd((name).c_str());
0066       TIter nextkey(gDirectory->GetListOfKeys());
0067       while (key = (TKey *)nextkey()) {
0068         obj = key->ReadObj();
0069         if (obj->IsA()->InheritsFrom("TH1")) {
0070           TH1 *h = (TH1 *)obj;
0071           TString fullpath = name + "/" + h->GetName();
0072           //std::cout << "j: " << j << " "<< h->GetName() <<" "<< fullpath << std::endl;
0073           ++j;
0074           if (TString(h->GetName()).Contains("Charge_Vs_Index"))
0075             continue;
0076           TCanvas *c1 = new TCanvas(h->GetName(), h->GetName(), 800, 600);
0077           c1->cd();
0078 
0079           TPad *pad1;
0080           if (obj->IsA()->InheritsFrom("TH2")) {
0081             pad1 = new TPad("pad1", "pad1", 0, 0., 0.5, 1.0);
0082           } else {
0083             pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
0084             pad1->SetBottomMargin(0.01);  // Upper and lower plot are joined
0085           }
0086           //      pad1->SetGridx();         // Vertical grid
0087           pad1->Draw();  // Draw the upper pad: pad1
0088           pad1->cd();    // pad1 becomes the current pad
0089 
0090           h->SetMarkerColor(kRed);
0091           h->SetStats(kFALSE);
0092           h->SetLineColor(kRed);
0093           h->SetMarkerStyle(kOpenSquare);
0094           h->GetXaxis()->SetLabelOffset(0.2);
0095 
0096           h->GetYaxis()->SetTitleSize(0.05);
0097           //h->GetYaxis()->SetTitleFont(43);
0098           h->GetYaxis()->SetTitleOffset(0.8);
0099 
0100           TH1 *h2 = (TH1 *)f2->Get(fullpath.Data());
0101           h2->SetStats(kFALSE);
0102 
0103           if (h2 == nullptr) {
0104             std::cout << "WARNING!!!!!! " << fullpath << " does NOT exist in second file!!!!!" << std::endl;
0105             continue;
0106           }
0107 
0108           h2->SetMarkerColor(kBlue);
0109           h2->SetLineColor(kBlue);
0110           h2->SetMarkerStyle(kOpenCircle);
0111           h2->GetXaxis()->SetLabelOffset(0.2);
0112           h2->GetYaxis()->SetTitleOffset(0.8);
0113 
0114           TObjArray *arrayHistos = new TObjArray();
0115           arrayHistos->Expand(2);
0116           arrayHistos->Add(h);
0117           arrayHistos->Add(h2);
0118 
0119           if (!obj->IsA()->InheritsFrom("TH2")) {
0120             Double_t theMaximum = getMaximum(arrayHistos);
0121             h->GetYaxis()->SetRangeUser(0., theMaximum * 1.30);
0122             h->Draw();
0123             h2->Draw("same");
0124           } else {
0125             c1->cd();
0126             pad1->cd();
0127             h->SetMarkerSize(0.1);
0128             h->Draw("colz");
0129 
0130             TLegend *lego = new TLegend(0.12, 0.80, 0.29, 0.88);
0131             lego->SetFillColor(10);
0132             lego->SetTextSize(0.042);
0133             lego->SetTextFont(42);
0134             lego->SetFillColor(10);
0135             lego->SetLineColor(10);
0136             lego->SetShadowColor(10);
0137             lego->AddEntry(h, leg1.Data());
0138             lego->Draw();
0139 
0140             c1->cd();
0141             TPad *pad2 = new TPad("pad2", "pad2", 0.5, 0., 1.0, 1.0);
0142             pad2->Draw();
0143             pad2->cd();
0144             h2->SetMarkerSize(0.1);
0145             h2->Draw("colz");
0146 
0147             TLegend *lego2 = new TLegend(0.12, 0.80, 0.29, 0.88);
0148             lego2->SetFillColor(10);
0149             lego2->SetTextSize(0.042);
0150             lego2->SetTextFont(42);
0151             lego2->SetFillColor(10);
0152             lego2->SetLineColor(10);
0153             lego2->SetShadowColor(10);
0154             lego2->AddEntry(h2, leg2.Data());
0155             lego2->Draw("same");
0156           }
0157 
0158           TString savename = fullpath.ReplaceAll("/", "_");
0159           double ksProb = 0;
0160 
0161           // lower plot will be in pad
0162           if (!obj->IsA()->InheritsFrom("TH2")) {
0163             TLegend *lego = new TLegend(0.12, 0.80, 0.29, 0.88);
0164             lego->SetFillColor(10);
0165             lego->SetTextSize(0.042);
0166             lego->SetTextFont(42);
0167             lego->SetFillColor(10);
0168             lego->SetLineColor(10);
0169             lego->SetShadowColor(10);
0170             lego->AddEntry(h, leg1.Data());
0171             lego->AddEntry(h2, leg2.Data());
0172 
0173             lego->Draw("same");
0174 
0175             c1->cd();  // Go back to the main canvas before defining pad2
0176 
0177             TPad *pad2 = new TPad("pad2", "pad2", 0, 0.05, 1, 0.3);
0178             pad2->SetTopMargin(0.01);
0179             pad2->SetBottomMargin(0.35);
0180             pad2->SetGridy();  // horizontal grid
0181             pad2->Draw();
0182             pad2->cd();  // pad2 becomes the current pad
0183 
0184             // Define the ratio plot
0185             TH1F *h3 = (TH1F *)h->Clone("h3");
0186             h3->SetLineColor(kBlack);
0187             h3->SetMarkerColor(kBlack);
0188             h3->SetTitle("");
0189             h3->SetMinimum(0.55);  // Define Y ..
0190             h3->SetMaximum(1.55);  // .. range
0191             h3->SetStats(0);       // No statistics on lower plot
0192             h3->Divide(h2);
0193             h3->SetMarkerStyle(20);
0194             h3->Draw("ep");  // Draw the ratio plot
0195 
0196             // Y axis ratio plot settings
0197             h3->GetYaxis()->SetTitle("ratio");
0198             h3->GetYaxis()->SetNdivisions(505);
0199             h3->GetYaxis()->SetTitleSize(20);
0200             h3->GetYaxis()->SetTitleFont(43);
0201             h3->GetYaxis()->SetTitleOffset(1.2);
0202             h3->GetYaxis()->SetLabelFont(43);  // Absolute font size in pixel (precision 3)
0203             h3->GetYaxis()->SetLabelSize(15);
0204 
0205             // X axis ratio plot settings
0206             h3->GetXaxis()->SetTitleSize(20);
0207             h3->GetXaxis()->SetTitleFont(43);
0208             h3->GetXaxis()->SetTitleOffset(4.);
0209             h3->GetXaxis()->SetLabelFont(43);  // Absolute font size in pixel (precision 3)
0210             h3->GetXaxis()->SetLabelSize(15);
0211             h3->GetXaxis()->SetLabelOffset(0.02);
0212           }
0213 
0214           // avoid spurious false positive due to empty histogram
0215           if (h->GetEntries() == 0 && h2->GetEntries() == 0) {
0216             ofs << "histogram # " << j << ": " << fullpath << " |has zero entries in both files: " << std::endl;
0217             delete c1;
0218             delete h;
0219             delete h2;
0220             continue;
0221           }
0222           // avoid doing k-test on histograms with different binning
0223           if (h->GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
0224               h->GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin()) {
0225             ofs << "histogram # " << j << ": " << fullpath << " |mismatched bins!!!!!!: " << std::endl;
0226             delete c1;
0227             delete h;
0228             delete h2;
0229             continue;
0230           }
0231 
0232           ksProb = h->KolmogorovTest(h2);
0233           //if(ksProb!=1.){
0234           //c1->SaveAs(savename+".pdf");
0235           TPaveText ksPt(0, 0.02, 0.35, 0.043, "NDC");
0236           ksPt.SetBorderSize(0);
0237           ksPt.SetFillColor(0);
0238           ksPt.AddText(Form("P(KS)=%g, ered %g eblue %g", ksProb, h->GetEntries(), h2->GetEntries()));
0239           if (!obj->IsA()->InheritsFrom("TH2")) {
0240             ksPt.SetTextSize(0.1);
0241             ksPt.Draw();
0242           } else {
0243             ksPt.SetTextSize(0.03);
0244             ksPt.Draw();
0245           }
0246           c1->Print("diff.pdf");
0247           std::cout << "histogram # " << j << ": " << fullpath << " |kolmogorov: " << ksProb << std::endl;
0248           ofs << "histogram # " << j << ": " << fullpath << " |kolmogorov: " << ksProb << std::endl;
0249           // }
0250 
0251           delete c1;
0252           delete h;
0253           delete h2;
0254         }
0255       }
0256     } else if (obj->IsA()->InheritsFrom("TH1")) {
0257       obj = key->ReadObj();
0258       TH1 *h = (TH1 *)obj;
0259       TString fullpath = (TString)h->GetName();
0260       //std::cout << "j: " << j << " "<< h->GetName() <<" "<< fullpath << std::endl;
0261       ++j;
0262       if (obj->IsA()->InheritsFrom("TH2"))
0263         continue;
0264       TCanvas *c1 = new TCanvas(h->GetName(), h->GetName(), 800, 600);
0265       c1->cd();
0266       TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
0267       pad1->SetBottomMargin(0.01);  // Upper and lower plot are joined
0268       //      pad1->SetGridx();         // Vertical grid
0269       pad1->Draw();  // Draw the upper pad: pad1
0270       pad1->cd();    // pad1 becomes the current pad
0271 
0272       h->SetMarkerColor(kRed);
0273       h->SetStats(kFALSE);
0274       h->SetLineColor(kRed);
0275       h->SetMarkerStyle(kOpenSquare);
0276       h->GetXaxis()->SetLabelOffset(0.2);
0277 
0278       h->GetYaxis()->SetTitleSize(0.05);
0279       //h->GetYaxis()->SetTitleFont(43);
0280       h->GetYaxis()->SetTitleOffset(0.8);
0281 
0282       TH1 *h2 = (TH1 *)f2->Get(fullpath.Data());
0283       h2->SetStats(kFALSE);
0284 
0285       if (h2 == nullptr) {
0286         std::cout << "WARNING!!!!!! " << fullpath << " does NOT exist in second file!!!!!" << std::endl;
0287         continue;
0288       }
0289 
0290       h2->SetMarkerColor(kBlue);
0291       h2->SetLineColor(kBlue);
0292       h2->SetMarkerStyle(kOpenCircle);
0293       h2->GetXaxis()->SetLabelOffset(0.2);
0294       h2->GetYaxis()->SetTitleOffset(0.8);
0295 
0296       TObjArray *arrayHistos = new TObjArray();
0297       arrayHistos->Expand(2);
0298       arrayHistos->Add(h);
0299       arrayHistos->Add(h2);
0300 
0301       Double_t theMaximum = getMaximum(arrayHistos);
0302       h->GetYaxis()->SetRangeUser(0., theMaximum * 1.30);
0303 
0304       h->Draw();
0305       h2->Draw("same");
0306       TString savename = fullpath.ReplaceAll("/", "_");
0307       double ksProb = 0;
0308 
0309       TLegend *lego = new TLegend(0.12, 0.80, 0.29, 0.88);
0310       lego->SetFillColor(10);
0311       lego->SetTextSize(0.042);
0312       lego->SetTextFont(42);
0313       lego->SetFillColor(10);
0314       lego->SetLineColor(10);
0315       lego->SetShadowColor(10);
0316       lego->AddEntry(h, leg1.Data());
0317       lego->AddEntry(h2, leg2.Data());
0318 
0319       lego->Draw("same");
0320 
0321       // lower plot will be in pad
0322       c1->cd();  // Go back to the main canvas before defining pad2
0323       TPad *pad2 = new TPad("pad2", "pad2", 0, 0.05, 1, 0.3);
0324       pad2->SetTopMargin(0.01);
0325       pad2->SetBottomMargin(0.35);
0326       pad2->SetGridy();  // horizontal grid
0327       pad2->Draw();
0328       pad2->cd();  // pad2 becomes the current pad
0329 
0330       // Define the ratio plot
0331       TH1F *h3 = (TH1F *)h->Clone("h3");
0332       h3->SetLineColor(kBlack);
0333       h3->SetMarkerColor(kBlack);
0334       h3->SetTitle("");
0335       h3->SetMinimum(0.55);  // Define Y ..
0336       h3->SetMaximum(1.55);  // .. range
0337       h3->SetStats(0);       // No statistics on lower plot
0338       h3->Divide(h2);
0339       h3->SetMarkerStyle(20);
0340       h3->Draw("ep");  // Draw the ratio plot
0341 
0342       // Y axis ratio plot settings
0343       h3->GetYaxis()->SetTitle("ratio");
0344       h3->GetYaxis()->SetNdivisions(505);
0345       h3->GetYaxis()->SetTitleSize(20);
0346       h3->GetYaxis()->SetTitleFont(43);
0347       h3->GetYaxis()->SetTitleOffset(1.2);
0348       h3->GetYaxis()->SetLabelFont(43);  // Absolute font size in pixel (precision 3)
0349       h3->GetYaxis()->SetLabelSize(15);
0350 
0351       // X axis ratio plot settings
0352       h3->GetXaxis()->SetTitleSize(20);
0353       h3->GetXaxis()->SetTitleFont(43);
0354       h3->GetXaxis()->SetTitleOffset(4.);
0355       h3->GetXaxis()->SetLabelFont(43);  // Absolute font size in pixel (precision 3)
0356       h3->GetXaxis()->SetLabelSize(15);
0357       h3->GetXaxis()->SetLabelOffset(0.02);
0358 
0359       // avoid spurious false positive due to empty histogram
0360       if (h->GetEntries() == 0 && h2->GetEntries() == 0) {
0361         ofs << "histogram # " << j << ": " << fullpath << " |has zero entries in both files: " << std::endl;
0362         delete c1;
0363         delete h;
0364         delete h2;
0365         continue;
0366       }
0367       // avoid doing k-test on histograms with different binning
0368       if (h->GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
0369           h->GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin()) {
0370         ofs << "histogram # " << j << ": " << fullpath << " |mismatched bins!!!!!!: " << std::endl;
0371         delete c1;
0372         delete h;
0373         delete h2;
0374         continue;
0375       }
0376 
0377       ksProb = h->KolmogorovTest(h2);
0378       //if(ksProb!=1.){
0379       //c1->SaveAs(savename+".pdf");
0380       TPaveText ksPt(0, 0.02, 0.35, 0.043, "NDC");
0381       ksPt.SetBorderSize(0);
0382       ksPt.SetFillColor(0);
0383       ksPt.AddText(Form("P(KS)=%g, ered %g eblue %g", ksProb, h->GetEntries(), h2->GetEntries()));
0384       ksPt.SetTextSize(0.1);
0385       ksPt.Draw();
0386       c1->Print("diff.pdf");
0387       std::cout << "histogram # " << j << ": " << fullpath << " |kolmogorov: " << ksProb << std::endl;
0388       ofs << "histogram # " << j << ": " << fullpath << " |kolmogorov: " << ksProb << std::endl;
0389       // }
0390 
0391       delete c1;
0392       delete h;
0393       delete h2;
0394     }
0395   }  //while
0396   f1->Close();
0397   f2->Close();
0398   dummyC.Print("diff.pdf]");
0399 
0400   ofs.close();
0401 }