EventGenerator/test/DetailedCompare.C

0001 // F. Cossutti
0002 //
0003 // ROOT macro for graphical compariosn of Monitor Elements in a user
0004 // supplied directory between two files with the same histogram content
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     // Set errors to sqrt(# entries)
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    // Draw a global picture title
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    // Draw reference
0251    pad1->cd();
0252    href_->SetYTitle("Events (A)");
0253    href_->DrawCopy("e1");
0254
0255
0256    // Draw new
0257    pad2->cd();
0258    hnew_->SetYTitle("Events (B)");
0259    hnew_->DrawCopy("e1");
0260
0261    gStyle->SetOptStat("nemruoi");
0262
0263    // Draw the two overlayed Normalized display region to 1
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      //     cout <<  v << "+/-" << e << " " << v/fabs(e) << endl;
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 }