Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:33:46

0001 ////////////////////////////////////////////////////////////
0002 //
0003 // A group of helper functions used to compare histograms
0004 // and profiles between two different releases
0005 //
0006 ////////////////////////////////////////////////////////////
0007 
0008 #include "TH1F.h"
0009 #include "TProfile.h"
0010 #include "TCanvas.h"
0011 #include "TMath.h"
0012 #include "TPad.h"
0013 #include "TPaveLabel.h"
0014 #include <vector>
0015 
0016 void NormalizeHistogramsToFirst(TH1* h1, TH1* h2);
0017 void NormalizeHistogramsToOne(TH1* h1, TH1* h2);
0018 void NormalizeHistogramsAsDensity(TH1* h1, TH1* h2);
0019 TH1* PlotRatiosHistograms(TH1* h1, TH1* h2);
0020 TH1* AddOverflow(TH1* h);
0021 
0022 int ratioCounter = 0;
0023 int overflowCounter = 0;
0024 
0025 // debugging printouts
0026 bool DEBUGP = false;
0027 
0028 ////////////////////////////////////////////////////////////
0029 //
0030 // Sets the global style for the whole system.
0031 //
0032 // Note: There might be redundancies in other places so
0033 //       things can be further simplified
0034 //
0035 void SetGlobalStyle() {
0036   gROOT->SetStyle("Plain");
0037   gStyle->SetPadGridX(kTRUE);
0038   gStyle->SetPadGridY(kTRUE);
0039   gStyle->SetPadRightMargin(0.07);
0040   gStyle->SetPadLeftMargin(0.13);
0041   //gStyle->SetTitleXSize(0.07); 
0042   //gStyle->SetTitleXOffset(0.6); 
0043   //tyle->SetTitleYSize(0.3);
0044   //gStyle->SetLabelSize(0.6) 
0045   //gStyle->SetTextSize(0.5);
0046   gStyle->SetOptStat(0);
0047 }
0048 //
0049 ////////////////////////////////////////////////////////////
0050 
0051 ////////////////////////////////////////////////////////////
0052 //
0053 // This function sets the style for a given histogram
0054 //
0055 // Note: There might be redundancies in other places so
0056 //       things can be further simplified
0057 //
0058 void SetHistogramStyle(TH1* h, Style_t mstyle, Color_t color, Size_t msize = 0.7, Width_t lwidth = 2,
0059                Float_t tsize = 0.05, Float_t toffset = 1.2) {
0060   if (!h)
0061     return;
0062   h->SetMarkerStyle(mstyle);
0063   h->SetMarkerColor(color);
0064   h->SetMarkerSize(msize);
0065   h->SetLineColor(color);
0066   h->SetLineWidth(lwidth);
0067   h->GetYaxis()->SetTitleSize(tsize);
0068   h->GetYaxis()->SetTitleOffset(toffset);
0069   h->GetXaxis()->SetLabelFont(63);
0070   h->GetXaxis()->SetLabelSize(14); // labels will be 14 pixels
0071   h->GetYaxis()->SetLabelFont(63);
0072   h->GetYaxis()->SetLabelSize(14);
0073 }
0074 //
0075 ////////////////////////////////////////////////////////////
0076 
0077 ////////////////////////////////////////////////////////////
0078 //
0079 // This function finds the list of TDirectories in a branch
0080 // that match a given string
0081 
0082 TList* GetListOfDirectories(TDirectory* dir, const char* match = 0) {
0083   TIter  nextkey(dir->GetListOfKeys());
0084   TList* sl     = new TList();
0085   TKey*  key    = 0;
0086   TKey*  oldkey = 0;
0087 
0088   while (( key = (TKey*)nextkey() )) {
0089     TObject *obj = key->ReadObj();
0090     if ( obj->IsA()->InheritsFrom( "TDirectory" ) ) {
0091       TString theName = obj->GetName();
0092       if (!match) {
0093     cout << " -> " << theName << endl;
0094     sl->Add(obj);
0095       }
0096       else if (theName.Contains(match)) {
0097     cout << " -> " << theName << endl;
0098     sl->Add(obj);
0099       }
0100     }
0101   } // End of while
0102 
0103   return sl;
0104 }
0105 
0106 ////////////////////////////////////////////////////////////
0107 //
0108 // This function goes to the right branch inside the file
0109 // looking for branches having branchContent.
0110 // It returns a list with all those branches,
0111 //
0112 TList* getListOfBranches(const char* dataType, TFile* file, const char* branchContent) {
0113   
0114   if (TString(dataType) == "RECO") {
0115     if(! file->cd("DQMData/Run 1/Muons/Run summary")) {
0116       cout << "ERROR: Muon Histos for " << dataType << " not found" << endl;
0117       return 0;
0118     }
0119   }
0120   else {
0121     cout << "ERROR: Data type " << dataType << " not allowed: only RECO is considered" << endl;
0122     return 0;
0123   }
0124 
0125   TDirectory * dir=gDirectory;
0126   TList* sl = GetListOfDirectories(dir, branchContent);
0127 
0128   if (sl->GetSize() == 0) {
0129     cout << "ERROR: No DQM muon reco histos found in NEW file " << endl;
0130     delete sl;
0131     return 0;
0132   }
0133 
0134   return sl;
0135 }
0136 //
0137 ////////////////////////////////////////////////////////////
0138 
0139 ////////////////////////////////////////////////////////////
0140 //
0141 // This function performs a compatibility test between two
0142 // histogram based on the Kolmogorov-Smirnof algorithm. It
0143 // also prints the value in a TPaveLabel at the upper-right
0144 // corner.
0145 // The return value contains the result of the test
0146 //
0147 double KolmogorovTest(TH1 *h1, TH1 *h2){
0148 
0149   double mya_array[1300], myb_array[1300];
0150   vector<double> mya;
0151   vector<double> myb;
0152   
0153   
0154   for (int i=0; i<h1->GetNbinsX(); i++){
0155     mya.push_back(h1->GetBinContent(i+1));
0156     myb.push_back(h2->GetBinContent(i+1));
0157   }
0158 
0159   sort(mya.begin(),mya.end());
0160   sort(myb.begin(),myb.end()); 
0161   copy(mya.begin(),mya.end(),mya_array);
0162   copy(myb.begin(),myb.end(),myb_array);
0163   
0164   const int nbinsa = h1->GetNbinsX();
0165   const int nbinsb = h2->GetNbinsX();
0166    
0167   double kstest = TMath::KolmogorovTest(nbinsa, mya_array,
0168                     nbinsb, myb_array,
0169                     "UOX");
0170   if (DEBUGP) cout << "   + KS value = " << kstest << endl;
0171 
0172   // Create text with the value
0173   TString legend = Form("KS=%4.2f", kstest);
0174 
0175   // Create a pave text to put the value inside
0176 
0177   TPaveLabel* pl = new TPaveLabel(0.79,0.91,0.93,0.96, legend.Data(), "NDC");
0178 
0179   // Tune style
0180   //pl->SetTextSize(0.04);
0181   pl->SetLineColor(41);
0182   pl->SetLineWidth(1);
0183   pl->SetLineStyle(1);
0184   pl->SetFillColor(41);
0185   pl->SetBorderSize(3);
0186 
0187   if (kstest < 0.7)
0188     pl->SetTextColor(kRed);
0189 
0190   pl->Draw();
0191   
0192   return kstest;
0193 }
0194 //
0195 ////////////////////////////////////////////////////////////
0196 
0197 ////////////////////////////////////////////////////////////
0198 //
0199 // This function draws the stat box for the two plots
0200 //
0201 void setStats(TH1* s, TH1* r, double startingY, double startingX = .1, bool fit = false){
0202   if (startingY<0){
0203     s->SetStats(0);
0204     r->SetStats(0);
0205   } 
0206   else {
0207     //gStyle->SetOptStat(1001);
0208     
0209     if (fit){
0210       s->Fit("gaus");
0211       TF1* f1 = (TF1*) s->GetListOfFunctions()->FindObject("gaus");
0212       if (f1) {
0213     f1->SetLineColor(2);
0214     f1->SetLineWidth(1);
0215       }
0216     }
0217     s->Draw();
0218     gPad->Update(); 
0219     TPaveStats* st1 = (TPaveStats*) s->GetListOfFunctions()->FindObject("stats");
0220     if (st1) {
0221       if (fit) {st1->SetOptFit(0010);    st1->SetOptStat(1001);}
0222       st1->SetX1NDC(startingX);
0223       st1->SetX2NDC(startingX+0.30);
0224       st1->SetY1NDC(startingY+0.20);
0225       st1->SetY2NDC(startingY+0.35);
0226       st1->SetTextColor(2);
0227     }
0228     else s->SetStats(0);
0229     if (fit) {
0230       r->Fit("gaus");
0231       TF1* f2 = (TF1*) r->GetListOfFunctions()->FindObject("gaus");
0232       if (f2) {
0233     f2->SetLineColor(4);
0234     f2->SetLineWidth(1);
0235       }
0236     }
0237     r->Draw();
0238     gPad->Update(); 
0239     TPaveStats* st2 = (TPaveStats*) r->GetListOfFunctions()->FindObject("stats");
0240     if (st2) {
0241       if (fit) {st2->SetOptFit(0010);    st2->SetOptStat(1001);}
0242       st2->SetX1NDC(startingX);
0243       st2->SetX2NDC(startingX+0.30);
0244       st2->SetY1NDC(startingY);
0245       st2->SetY2NDC(startingY+0.15);
0246       st2->SetTextColor(4);
0247     }
0248     else r->SetStats(0);
0249   }
0250 }
0251 //
0252 ////////////////////////////////////////////////////////////
0253 
0254 ////////////////////////////////////////////////////////////
0255 //
0256 // Plot a page with several histograms
0257 //
0258 void PlotNHistograms(const TString& pdfFile,
0259              TDirectory* rdir, TDirectory* sdir,
0260              const TString& rcollname, const TString& scollname,
0261              const char* canvasName, const char* canvasTitle,
0262              const TString& refLabel, const TString& newLabel,
0263              unsigned int nhistos, const TString* hnames,
0264              const TString* htitles, const char** drawopt,
0265              bool* logy = 0,  bool* logx = 0, bool* doKolmo = 0, Double_t* norm = 0,
0266              Double_t* minx = 0, Double_t* maxx = 0,
0267              Double_t* miny = 0, Double_t* maxy = 0) {
0268 
0269   if (DEBUGP) {
0270     cout << "   + Plotting histograms for " << canvasTitle << endl;
0271     cerr << "   + Plotting histograms for " << canvasTitle << endl;
0272   }
0273 
0274   TH1* rh_raw   = 0;
0275   TH1* rh       = 0;
0276   TH1* sh_raw   = 0;
0277   TH1* sh       = 0;
0278 
0279   TCanvas* canvas =  0;
0280   if (nhistos >4)
0281     canvas = new TCanvas(canvasName, canvasTitle, 1000, 1400);
0282   else
0283     canvas = new TCanvas(canvasName, canvasTitle, 1000, 1050);
0284 
0285   canvas->Draw();
0286   canvas->Divide(2,(nhistos+1)/2); //This way we print in 2 columns
0287 
0288   for (unsigned int i = 0; i < nhistos; i++) {
0289 
0290     if (DEBUGP) cout << " [" << i << "] histogram name: " << flush << hnames[i] << endl;
0291 
0292     //draw option for the new histogram
0293     TString drawoption = drawopt[i]; 
0294 
0295     // Skip histogram if no name is provided
0296     if (hnames[i] == "") continue;
0297     
0298     // Get Histograms
0299     // + Reference release
0300     rdir->cd(rcollname);
0301     
0302     if (DEBUGP) cout << " Getting object for reference sample " << (rcollname + "/" + hnames[i]) << endl;
0303 
0304     rdir->GetObject(rcollname + "/" + hnames[i], rh_raw);
0305     if (! rh_raw) {
0306       cout << "WARNING: Could not find a reference histogram or profile named " << hnames[i]
0307        << " in " << rdir->GetName() << endl;
0308       cout << "         Skipping" << endl;
0309       continue;
0310     }
0311 
0312     //If it is a 2D project it in Y... is this what we always want?
0313     if (TString(rh_raw->IsA()->GetName()) == "TH2F") {
0314 
0315       if (DEBUGP) cout << " It is a TH2F object... project in Y!" << endl;
0316 
0317       TH1* proj = ((TH2F*) rh_raw)->ProjectionY();
0318       rh_raw = proj;
0319     }
0320 
0321     // + New release    
0322     sdir->cd(scollname);
0323 
0324     if (DEBUGP) cout << " Getting object for target sample " << (scollname + "/" + hnames[i]) << endl;
0325 
0326     sdir->GetObject(scollname + "/" + hnames[i], sh_raw);
0327     if (! sh_raw) {
0328       cout << "WARNING: Could not find a signal histogram or profile named " << hnames[i] 
0329        << " in " << sdir->GetName() << endl;
0330       cout << "         Skipping" << endl;
0331       continue;
0332     }
0333 
0334     //If it is a 2D project it in Y... is this what we always want?
0335     if (TString(sh_raw->IsA()->GetName()) == "TH2F") {
0336 
0337       if (DEBUGP) cout << hnames[i] << " is a TH2F object... project in Y!" << endl;
0338 
0339       TH1* proj = ((TH2F*) sh_raw)->ProjectionY();
0340       sh_raw = proj;
0341     }
0342 
0343     rh = AddOverflow(rh_raw);
0344     sh = AddOverflow(sh_raw);
0345 
0346     // Set styles
0347 
0348     if (DEBUGP) cout << " Setting style..." << endl;
0349 
0350     SetHistogramStyle(rh, 21, 4);
0351     SetHistogramStyle(sh, 20, 2);
0352 
0353     //Change titles
0354     if (htitles) {
0355       rh->SetTitle(htitles[i]);
0356       sh->SetTitle(htitles[i]);
0357       rh->GetYaxis()->SetTitle(htitles[i]);
0358       sh->GetYaxis()->SetTitle(htitles[i]);
0359     }
0360 
0361     // SET X AXIS RANGE in plots
0362     Bool_t ChangeXRange = false ;
0363     Double_t Xleft = rh->GetXaxis()->GetXmin();
0364     Double_t Xright = rh->GetXaxis()->GetXmax(); 
0365 
0366     if (DEBUGP) cout << "ref histo Xleft, Xright = "<< Xleft << ", "<< Xright  << endl;
0367    
0368     if (sh->GetXaxis()->GetXmin() < rh->GetXaxis()->GetXmin()) {
0369       Xleft = sh->GetXaxis()->GetXmin();
0370       ChangeXRange = true;
0371       if (DEBUGP) cout << "automatic reset MIN (new < ref) Xleft = "<< Xleft << endl;
0372     }
0373     if (sh->GetXaxis()->GetXmax() > rh->GetXaxis()->GetXmax()) {
0374       Xright = sh->GetXaxis()->GetXmax();
0375       ChangeXRange = true;
0376       if (DEBUGP) cout << "automatic reset MAX (new > ref) Xright = "<< Xright << endl;
0377     }
0378 
0379     if (minx[i]!=0) {
0380       ChangeXRange = true ;
0381       Xleft = minx[i];
0382       if (DEBUGP) cout << "user reset Xleft = "<< Xleft << endl;
0383     }
0384     
0385     if (maxx[i]!=0) {
0386       ChangeXRange = true ;
0387       Xright = maxx[i];
0388       if (DEBUGP) cout << "user reset Xright = "<< Xleft << endl;
0389     }
0390         
0391     if (ChangeXRange) {
0392       if (DEBUGP) {
0393     cout << "Ref histo Xmin, Xmax = "<< rh->GetXaxis()->GetXmin()  << ", " << rh->GetXaxis()->GetXmax() <<endl;
0394     cout << "New histo Xmin, Xmax = "<< sh->GetXaxis()->GetXmin()  << ", " << sh->GetXaxis()->GetXmax() <<endl;
0395       }
0396       
0397       rh->GetXaxis()->SetRangeUser(Xleft,Xright);
0398       sh->GetXaxis()->SetRangeUser(Xleft,Xright);      
0399 
0400       if (DEBUGP) {
0401     cout << "reset Ref histo Xmin, Xmax = "<< rh->GetXaxis()->GetXmin()  << ", " << rh->GetXaxis()->GetXmax() <<endl;
0402     cout << "reset New histo Xmin, Xmax = "<< sh->GetXaxis()->GetXmin()  << ", " << sh->GetXaxis()->GetXmax() <<endl;
0403     cout << "resetting Ref and New histo Xleft, Xright = "<< Xleft << ", " << Xright <<endl;
0404       }
0405     }
0406 
0407     // ===============================================================================================
0408     // Normalize
0409     if (norm[i] < 0.) ; 
0410       // Default: do not normalize at all !
0411     else if (norm[i] == 0.)
0412       NormalizeHistogramsToFirst(rh,sh);
0413     else if (norm[i] == 1.)
0414       NormalizeHistogramsToOne(rh,sh);
0415     else if (norm[i] == 2.)
0416       NormalizeHistogramsAsDensity(rh,sh);
0417     // ===============================================================================================
0418 
0419     // ===============================================================================================
0420     // SET Y AXIS RANGE in plots
0421     //
0422     // MINIMUM
0423     //
0424     Double_t Ybottom;
0425 
0426     // if user-defined range force it !
0427     if (miny[i]!=0) {
0428       Ybottom = miny[i];
0429       if (DEBUGP) cout << "setting Minimum Y to user defined value: "<< miny[i] << endl;
0430     }
0431     else if (logy[i]) {
0432       // automatic setting for log scale
0433       Double_t yminr = rh->GetMinimum(0.); // min value larger than zero
0434       Double_t ymins = sh->GetMinimum(0.);
0435       Ybottom = yminr < ymins ? yminr*0.5 : ymins*0.5;
0436       if (DEBUGP) cout << "LOG scale, yminr, ymins: "<<yminr<<", "<<ymins <<"  ==>> Ybottom = "<<Ybottom<< endl;
0437     }
0438     else {
0439       // automatic setting for linear scale
0440       Double_t yminr = rh->GetMinimum(); // min value larger than zero
0441       Double_t ymins = sh->GetMinimum();
0442       Ybottom = yminr < ymins ? yminr-0.1*abs(yminr) : ymins-0.1*abs(ymins) ;
0443       // limit the scale to -1,+1 for relative pt bias to avoid failing fits
0444       if ((hnames[i] == "ptres_vs_eta_Mean") && (Ybottom <-1.)) Ybottom = -1.;
0445       if ((hnames[i] == "ptres_vs_pt_Mean") && (Ybottom <-1.)) Ybottom = -1.;
0446       if (DEBUGP) cout << "LIN scale, yminr, ymins: "<<yminr<<", "<<ymins <<"  ==>> Ybottom = "<<Ybottom<< endl;    
0447     }
0448 
0449     ///////////////////
0450     // MAXIMUM
0451     //
0452     Double_t Ytop;
0453 
0454     // if user-defined range force it !
0455     if (maxy[i]!=0) {
0456       Ytop = maxy[i];
0457       if (DEBUGP) cout << "setting Maximum Y to user defined value: "<< maxy[i] << endl;
0458     }
0459     else {
0460       Double_t ymaxr = rh->GetMaximum(); // max value
0461       Double_t ymaxs = sh->GetMaximum();
0462       Ytop = ymaxr > ymaxs ? ymaxr : ymaxs ;
0463       // automatic setting for log scale
0464       if (logy[i]) {
0465     Ytop = Ytop*2;
0466     if (DEBUGP) cout << "LOG scale, ymaxr, ymaxs: "<<ymaxr<<", "<<ymaxs <<"  ==>> Ytop = "<<Ytop<< endl;
0467       }
0468       else {
0469     Ytop = Ytop+0.1*abs(Ytop);
0470     // limit the scale to -1,+1 for relative pt bias to avoid failing fits
0471     if ((hnames[i] == "ptres_vs_eta_Mean") && (Ytop >1.)) Ytop = 1.;
0472     if ((hnames[i] == "ptres_vs_pt_Mean") && (Ytop >1.)) Ytop = 1.;
0473     if (DEBUGP) cout << "LIN scale, ymaxr, ymaxs: "<<ymaxr<<", "<<ymaxs <<"  ==>> Ytop = "<<Ytop<< endl;     
0474       }
0475     }
0476 
0477     // +++++++++++++++++++++++++++++++++++++++++
0478     rh->GetYaxis()->SetRangeUser(Ybottom,Ytop);
0479     sh->GetYaxis()->SetRangeUser(Ybottom,Ytop); 
0480     // +++++++++++++++++++++++++++++++++++++++++     
0481 
0482     // Move to subpad
0483     canvas->cd(i+1);
0484     
0485     TPad* pad1 = NULL;
0486     TPad* pad2 = NULL;
0487 
0488     pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
0489     pad2 = new TPad("pad2", "pad2", 0, 0.0, 1, 0.3);
0490 
0491     pad1->SetTopMargin   (0.08);
0492     pad1->SetBottomMargin(0.01);
0493     pad1->Draw();
0494 
0495     pad2->SetTopMargin   (0.05);
0496     pad2->SetBottomMargin(0.45);
0497     pad2->Draw();// Set stat boxes                                                                      
0498     pad1->cd();
0499 
0500     // Check Logy                                                                                       
0501     if (logy[i]) gPad->SetLogy();
0502     if (logx[i]) {gPad->SetLogx(); pad2->SetLogx();}
0503 
0504     // Set stat boxes
0505 
0506     if (DEBUGP) cout << " Setting statistics..." << endl;
0507     setStats(sh, rh, -1, 0, false);
0508    
0509     // /////////// DRAW  histograms //////////////////////////////////////
0510     //
0511     // FIRST plot: reference (blue)   SECOND plot: new (red)
0512     if (DEBUGP) cout << " Drawing histograms..." << endl;
0513 
0514     if (ChangeXRange) {
0515       sh->Draw(drawoption);
0516       rh->Draw("same"+drawoption);
0517       sh->Draw("same"+drawoption);
0518     }
0519     else {
0520       rh->Draw(drawoption);
0521       sh->Draw("same"+drawoption);
0522     }
0523 
0524     // Perform Kolmogorov test if needed
0525     if (doKolmo) {
0526       if (doKolmo[i]) {
0527     if (DEBUGP) cout << " Performing Kolmogorov test..." << endl;
0528     //  TPad* c1_1 = canvas->GetPad(i+1);
0529     double kstest = KolmogorovTest(sh,rh);
0530     if(kstest<0.7)
0531       gPad->SetFillColor(kBlue-10);
0532       }
0533     }
0534  
0535     pad2->cd();
0536 
0537     TH1* ratioplot = PlotRatiosHistograms(rh, sh);
0538     SetHistogramStyle(ratioplot, 21, 4);
0539     ratioplot->Draw("ep");
0540  } // End loop
0541   
0542    // Draw Legend
0543 
0544   if (DEBUGP) cout << " Drawing legend..." << endl;
0545 
0546   canvas->cd();
0547   
0548   TLegend* l = 0;
0549   if (nhistos > 4)
0550     l = new TLegend(0.20,0.665,0.80,0.685);
0551   else
0552     l = new TLegend(0.20,0.50,0.80,0.53);
0553 
0554   l->SetTextSize(0.011);
0555   l->SetLineColor(1);
0556   l->SetLineWidth(1);
0557   l->SetLineStyle(1);
0558   l->SetFillColor(0);
0559   l->SetBorderSize(2);
0560   l->AddEntry(rh,refLabel,"LPF");
0561   l->AddEntry(sh,newLabel,"LPF");
0562   l->Draw();
0563   
0564   // Print Canvas
0565   canvas->SaveAs(pdfFile+".pdf");
0566   canvas->SaveAs(pdfFile+".png");
0567 
0568   // Clean memory
0569   // delete l;
0570   delete canvas;
0571   if (DEBUGP) cout << "     ... plotted histograms for " << canvasTitle << endl;
0572 }
0573 //
0574 ////////////////////////////////////////////////////////////
0575 
0576 
0577 ////////////////////////////////////////////////////////////
0578 //
0579 // Plot a page with 4 histograms
0580 //
0581 void Plot4Histograms(const TString& pdfFile,
0582              TDirectory* rdir, TDirectory* sdir,
0583              const TString& rcollname, const TString& scollname,
0584              const char* canvasName, const char* canvasTitle,
0585              const TString& refLabel, const TString& newLabel,
0586              const TString* hnames, const TString* htitles, const char** drawopt,
0587              bool* logy = 0, bool* logx = 0, bool* doKolmo = 0, Double_t* norm = 0,
0588              Double_t* minx = 0, Double_t* maxx = 0,
0589              Double_t* miny = 0, Double_t* maxy = 0) {
0590   PlotNHistograms(pdfFile,
0591           rdir, sdir,
0592           rcollname, scollname,
0593           canvasName, canvasTitle,
0594           refLabel, newLabel,
0595           4, hnames, htitles, drawopt,
0596           logy, logx, doKolmo, norm, minx, maxx, miny, maxy);
0597 
0598 }
0599 //
0600 ////////////////////////////////////////////////////////////
0601 
0602 
0603 ////////////////////////////////////////////////////////////
0604 //
0605 // Plot a page with 6 histograms
0606 //
0607 void Plot6Histograms(const TString& pdfFile,
0608              TDirectory* rdir, TDirectory* sdir,
0609              const TString& rcollname, const TString& scollname,
0610              const char* canvasName, const char* canvasTitle,
0611              const TString& refLabel, const TString& newLabel,
0612              const TString* hnames, const TString* htitles, const char** drawopt,
0613              bool* logy = 0, bool* logx = 0, bool* doKolmo = 0, Double_t* norm = 0,
0614              Double_t* minx = 0, Double_t* maxx = 0,
0615              Double_t* miny = 0, Double_t* maxy = 0) {
0616   
0617 
0618   PlotNHistograms(pdfFile,
0619           rdir, sdir,
0620           rcollname, scollname,
0621           canvasName, canvasTitle,
0622           refLabel, newLabel,
0623           6, hnames, htitles, drawopt,
0624           logy, logx, doKolmo, norm,
0625           minx, maxx,
0626           miny, maxy);
0627 }
0628 //
0629 ////////////////////////////////////////////////////////////
0630 
0631 void Plot5Histograms(const TString& pdfFile,
0632              TDirectory* rdir, TDirectory* sdir,
0633              const TString& rcollname, const TString& scollname,
0634              const char* canvasName, const char* canvasTitle,
0635              const TString& refLabel, const TString& newLabel,
0636              const TString* hnames, const TString* htitles, const char** drawopt,
0637              bool* logy = 0,  bool* logx = 0, bool* doKolmo = 0, Double_t* norm = 0,
0638              Double_t* minx = 0, Double_t* maxx = 0,
0639              Double_t* miny = 0, Double_t* maxy = 0) {
0640   
0641 
0642   PlotNHistograms(pdfFile,
0643           rdir, sdir,
0644           rcollname, scollname,
0645           canvasName, canvasTitle,
0646           refLabel, newLabel,
0647           5, hnames, htitles, drawopt,
0648           logy, logx, doKolmo, norm,
0649           minx, maxx,
0650           miny, maxy);
0651 }
0652 //
0653 ////////////////////////////////////////////////////////////
0654 
0655 
0656 ////////////////////////////////////////////////////////////
0657 //
0658 // Normalize the two histograms to the entries on the first histogram
0659 //
0660 void NormalizeHistogramsToFirst(TH1* h1, TH1* h2) {
0661   if (h1==0 || h2==0) return;
0662   
0663   if ( h1->Integral() > 0 && h2->Integral() > 0 ){
0664     Double_t scale2 = h1->Integral()/h2->Integral();
0665     h2->Scale(scale2);
0666   }
0667 }
0668 //
0669 ////////////////////////////////////////////////////////////
0670 
0671 
0672 ////////////////////////////////////////////////////////////
0673 //
0674 // Normalize the two histograms to unity
0675 //
0676 void NormalizeHistogramsToOne(TH1* h1, TH1* h2) {
0677   if (!h1 || !h2) return;
0678   
0679   if ( h1->Integral() != 0 && h2->Integral() != 0 ) {
0680     Double_t scale1 = 1.0/h1->Integral();
0681     Double_t scale2 = 1.0/h2->Integral();
0682     h1->Scale(scale1);
0683     h2->Scale(scale2);
0684   }
0685 }
0686 //
0687 ////////////////////////////////////////////////////////////
0688 
0689 // Normalize the two histograms as probability density functions:
0690 // Normalize areas to unity and consider the bin width
0691 //
0692 void NormalizeHistogramsAsDensity(TH1* h1, TH1* h2) {
0693   if (h1==0 || h2==0) return;
0694   if (h1->Integral() != 0 && h2->Integral() != 0 ) {
0695     Double_t scale1 = 1.0/h1->Integral();
0696     Double_t scale2 = 1.0/h2->Integral();
0697     h1->Scale(scale1, "width");
0698     h2->Scale(scale2, "width");
0699   }
0700 }
0701 ///////////////////////////////////////////////////////////
0702 // 
0703 // ratio plot from the two histograms
0704 //  
0705 /////////////////////////////////////////////////////////
0706 
0707 TH1* PlotRatiosHistograms(TH1* h1, TH1* h2){
0708 
0709   ++ratioCounter;
0710 
0711   const Int_t nbinsx = h1->GetNbinsX();
0712 
0713   Double_t xbins[nbinsx+1];
0714 
0715   //Loop necessary since some histograms are TH1 while some are TProfile and do not allow us to use the clone function (SetBinContent do not exists for TProfiles)
0716   for (Int_t ibin=1; ibin<=nbinsx+1; ibin++) {
0717     Float_t xmin = h1->GetBinLowEdge(ibin);
0718     xbins[ibin-1] = xmin;
0719   }
0720 
0721 
0722   TH1F* h_ratio = new TH1F(Form("h_ratio_%d", ratioCounter), "", nbinsx, xbins);
0723 
0724   for (Int_t ibin=1; ibin<=nbinsx; ibin++) {
0725 
0726     Float_t h1Value = h1->GetBinContent(ibin);
0727     Float_t h2Value = h2->GetBinContent(ibin);
0728 
0729     Float_t h1Error = h1->GetBinError(ibin);
0730     Float_t h2Error = h2->GetBinError(ibin);
0731 
0732     Float_t ratioVal = 999;
0733     Float_t ratioErr = 999;
0734 
0735     if (h2Value > 0 || h2Value < 0) {
0736       ratioVal = h1Value / h2Value;
0737       ratioErr = h1Error / h2Value;
0738     }
0739 
0740     h_ratio->SetBinContent(ibin, ratioVal);
0741     h_ratio->SetBinError  (ibin, ratioErr);
0742 
0743   }
0744 
0745   h_ratio->SetTitle("");
0746   h_ratio->GetYaxis()->SetTitle("");
0747   h_ratio->GetYaxis()->SetRangeUser(0.4, 1.6);
0748 
0749   return h_ratio;
0750 }
0751 
0752 TH1* AddOverflow(TH1* h) {
0753 
0754   ++overflowCounter;
0755 
0756   TString  name = h->GetName();
0757   const Int_t nx = h->GetNbinsX();
0758   Double_t x1   = h->GetBinLowEdge(1);
0759 
0760   Double_t xbins[nx+1];
0761 
0762   //Loop necessary since some histograms are TH1 while some are TProfile and do not allow us to use the clone function (SetBinContent do not exists for TProfiles)
0763   for (Int_t ibin=1; ibin<=nx+1; ibin++) {
0764     Float_t xmin = h->GetBinLowEdge(ibin);
0765     xbins[ibin-1] = xmin;
0766   }
0767 
0768   TH1F* htmp = new TH1F(Form(name + "_overflow_%d", overflowCounter), "", nx, xbins);
0769 
0770   // Fill the new histogram including the extra bin for overflows
0771   for (Int_t i=1; i<=nx; i++) {
0772     htmp->Fill(htmp->GetBinCenter(i), h->GetBinContent(i));
0773     htmp->SetBinError(i, h->GetBinError(i));
0774   }
0775 
0776   // Fill the underflow
0777   htmp->Fill(x1-1, h->GetBinContent(0));
0778 
0779   // Restore the number of entries
0780   htmp->SetEntries(h->GetEntries());
0781 
0782   // Cosmetics
0783   htmp->SetLineColor(h->GetLineColor());
0784   htmp->SetLineWidth(h->GetLineWidth());
0785   htmp->GetXaxis()->SetTitleOffset(1.5);
0786 
0787   return htmp;
0788 }