Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:07:12 

0001 ////////////////////////////////////////////////////////////
0002 //
0003 // A group of helper functions used to compare histograms
0004 // and profiles between two different releases
0005 //
0006 ////////////////////////////////////////////////////////////
0007 
0008 /////
0009 // Some includes
0010 #include "TH1F.h"
0011 #include "TProfile.h"
0012 #include "TCanvas.h"
0013 #include "TMath.h"
0014 #include "TPad.h"
0015 #include "TPaveLabel.h"
0016 #include <vector>
0017 
0018 void NormalizeHistogramsToFirst(TH1* h1, TH1* h2);
0019 void NormalizeHistogramsTo1(TH1* h1, TH1* h2);
0020 TH1* PlotRatiosHistograms(TH1* h1, TH1* h2);
0021 
0022 int ratioCounter = 0;
0023 
0024 
0025 /////
0026 // Uncomment the following line to get more debuggin output
0027 // #define DEBUG
0028 
0029 ////////////////////////////////////////////////////////////
0030 //
0031 // Sets the global style for the whole system.
0032 //
0033 // Note: There might be redundancies in other places so
0034 //       things can be further simplified
0035 //
0036 void SetGlobalStyle() {
0037   gROOT->SetStyle("Plain");
0038   gStyle->SetPadGridX(kTRUE);
0039   gStyle->SetPadGridY(kTRUE);
0040   gStyle->SetPadRightMargin(0.07);
0041   gStyle->SetPadLeftMargin(0.13);
0042   //gStyle->SetTitleXSize(0.07); 
0043   //gStyle->SetTitleXOffset(0.6); 
0044   //tyle->SetTitleYSize(0.3);
0045   //gStyle->SetLabelSize(0.6) 
0046   //gStyle->SetTextSize(0.5);
0047   gStyle->SetOptStat(0);
0048 }
0049 //
0050 ////////////////////////////////////////////////////////////
0051 
0052 ////////////////////////////////////////////////////////////
0053 //
0054 // This function sets the style for a given histogram
0055 //
0056 // Note: There might be redundancies in other places so
0057 //       things can be further simplified
0058 //
0059 void SetHistogramStyle(TH1* h, Style_t mstyle, Color_t color, Size_t msize = 0.7, Width_t lwidth = 2,
0060                Float_t tsize = 0.05, Float_t toffset = 1.2) {
0061   if (!h)
0062     return;
0063   h->SetMarkerStyle(mstyle);
0064   h->SetMarkerColor(color);
0065   h->SetMarkerSize(msize);
0066   h->SetLineColor(color);
0067   h->SetLineWidth(lwidth);
0068   h->GetYaxis()->SetTitleSize(tsize);
0069   h->GetYaxis()->SetTitleOffset(toffset);
0070   h->GetXaxis()->SetLabelFont(63);
0071   h->GetXaxis()->SetLabelSize(14); // labels will be 14 pixels
0072   h->GetYaxis()->SetLabelFont(63);
0073   h->GetYaxis()->SetLabelSize(14);
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 (dataType == "HLT") {
0115     if(file->cd("DQMData/Run 1/HLT")) {
0116       file->cd("DQMData/Run 1/HLT/Run summary/Muon/MultiTrack");
0117     }
0118     else {
0119       file->cd("DQMData/HLT/Muon/MultiTrack");
0120     }
0121   }
0122   else if (dataType == "RECO") {
0123     if(file->cd("DQMData/Run 1/RecoMuonV")) {
0124       file->cd("DQMData/Run 1/RecoMuonV/Run summary/MultiTrack");
0125     }
0126     else if(file->cd("DQMData/Run 1/Muons/Run summary/RecoMuonV")) {
0127       file->cd("DQMData/Run 1/Muons/Run summary/RecoMuonV/MultiTrack");
0128     }
0129     else {
0130       file->cd("DQMData/RecoMuonV/MultiTrack");
0131     }
0132   }
0133   else {
0134     cout << "ERROR: Data type " << dataType << " not allowed: only RECO and HLT are considered" << endl;
0135     cerr << "ERROR: Data type " << dataType << " not allowed: only RECO and HLT are considered" << endl;
0136     return;
0137   }
0138   */
0139   
0140   if (TString(dataType) == "RECO") {
0141     if(! file->cd("DQMData/Run 1/Muons/Run summary")) {
0142       cout << "ERROR: Muon Histos for " << dataType << " not found" << endl;
0143       return 0;
0144     }
0145   }
0146   else {
0147     cout << "ERROR: Data type " << dataType << " not allowed: only RECO is considered" << endl;
0148     return 0;
0149   }
0150 
0151   TDirectory * dir=gDirectory;
0152   TList* sl = GetListOfDirectories(dir, branchContent);
0153 
0154   if (sl->GetSize() == 0) {
0155     cout << "ERROR: No DQM muon reco histos found in NEW file " << endl;
0156     delete sl;
0157     return 0;
0158   }
0159 
0160   return sl;
0161 }
0162 //
0163 ////////////////////////////////////////////////////////////
0164 
0165 ////////////////////////////////////////////////////////////
0166 //
0167 // This function performs a compatibility test between two
0168 // histogram based on the Kolmogorov-Smirnof algorithm. It
0169 // also prints the value in a TPaveLabel at the upper-right
0170 // corner.
0171 // The return value contains the result of the test
0172 //
0173 double KolmogorovTest(TH1 *h1, TH1 *h2){
0174 
0175   double mya_array[1300], myb_array[1300];
0176   vector<double> mya;
0177   vector<double> myb;
0178   
0179   
0180   for (int i=0; i<h1->GetNbinsX(); i++){
0181     mya.push_back(h1->GetBinContent(i+1));
0182     myb.push_back(h2->GetBinContent(i+1));
0183   }
0184 
0185   sort(mya.begin(),mya.end());
0186   sort(myb.begin(),myb.end()); 
0187   copy(mya.begin(),mya.end(),mya_array);
0188   copy(myb.begin(),myb.end(),myb_array);
0189   
0190   const int nbinsa = h1->GetNbinsX();
0191   const int nbinsb = h2->GetNbinsX();
0192    
0193   double kstest = TMath::KolmogorovTest(nbinsa, mya_array,
0194                     nbinsb, myb_array,
0195                     "UOX");
0196 #ifdef DEBUG
0197   cout << "DEBUG:   + KS value = " << kstest << endl;
0198 #endif
0199 
0200 
0201   // Create text with the value
0202   TString legend = Form("KS=%4.2f", kstest);
0203 
0204   // Create a pave text to put the value inside
0205 
0206   TPaveLabel* pl = new TPaveLabel(0.79,0.91,0.93,0.96, legend.Data(), "NDC");
0207 
0208   // Tune style
0209   //pl->SetTextSize(0.04);
0210   pl->SetLineColor(41);
0211   pl->SetLineWidth(1);
0212   pl->SetLineStyle(1);
0213   pl->SetFillColor(41);
0214   pl->SetBorderSize(3);
0215 
0216   if (kstest < 0.7)
0217     pl->SetTextColor(kRed);
0218 
0219   pl->Draw();
0220   
0221   return kstest;
0222 }
0223 //
0224 ////////////////////////////////////////////////////////////
0225 
0226 ////////////////////////////////////////////////////////////
0227 //
0228 // This function draws the stat box for the two plots
0229 //
0230 void setStats(TH1* s, TH1* r, double startingY, double startingX = .1, bool fit = false){
0231   if (startingY<0){
0232     s->SetStats(0);
0233     r->SetStats(0);
0234   } 
0235   else {
0236     //gStyle->SetOptStat(1001);
0237     
0238     if (fit){
0239       s->Fit("gaus");
0240       TF1* f1 = (TF1*) s->GetListOfFunctions()->FindObject("gaus");
0241       if (f1) {
0242     f1->SetLineColor(2);
0243     f1->SetLineWidth(1);
0244       }
0245     }
0246     s->Draw();
0247     gPad->Update(); 
0248     TPaveStats* st1 = (TPaveStats*) s->GetListOfFunctions()->FindObject("stats");
0249     if (st1) {
0250       if (fit) {st1->SetOptFit(0010);    st1->SetOptStat(1001);}
0251       st1->SetX1NDC(startingX);
0252       st1->SetX2NDC(startingX+0.30);
0253       st1->SetY1NDC(startingY+0.20);
0254       st1->SetY2NDC(startingY+0.35);
0255       st1->SetTextColor(2);
0256     }
0257     else s->SetStats(0);
0258     if (fit) {
0259       r->Fit("gaus");
0260       TF1* f2 = (TF1*) r->GetListOfFunctions()->FindObject("gaus");
0261       if (f2) {
0262     f2->SetLineColor(4);
0263     f2->SetLineWidth(1);
0264       }
0265     }
0266     r->Draw();
0267     gPad->Update(); 
0268     TPaveStats* st2 = (TPaveStats*) r->GetListOfFunctions()->FindObject("stats");
0269     if (st2) {
0270       if (fit) {st2->SetOptFit(0010);    st2->SetOptStat(1001);}
0271       st2->SetX1NDC(startingX);
0272       st2->SetX2NDC(startingX+0.30);
0273       st2->SetY1NDC(startingY);
0274       st2->SetY2NDC(startingY+0.15);
0275       st2->SetTextColor(4);
0276     }
0277     else r->SetStats(0);
0278   }
0279 }
0280 //
0281 ////////////////////////////////////////////////////////////
0282 
0283 ////////////////////////////////////////////////////////////
0284 //
0285 // Plot a page with several histograms
0286 //
0287 void PlotNHistograms(const TString& pdfFile,
0288              TDirectory* rdir, TDirectory* sdir,
0289              const TString& rcollname, const TString& scollname,
0290              const char* canvasName, const char* canvasTitle,
0291              const TString& refLabel, const TString& newLabel,
0292              unsigned int nhistos, const char** hnames, const char** htitles = 0,
0293              bool* logy = 0, bool* doKolmo = 0, Double_t* norm = 0, bool *resol = 0,
0294              Double_t* minx = 0, Double_t* maxx = 0,
0295              Double_t* miny = 0, Double_t* maxy = 0) {
0296 #ifdef DEBUG
0297   cout << "   + Plotting histograms for " << canvasTitle << endl;
0298   cerr << "   + Plotting histograms for " << canvasTitle << endl;
0299 #endif
0300   TH1* rh   = 0;
0301   TH1* sh   = 0;
0302   // TH2* rh2d = 0;
0303   // TH2* sh2d = 0;
0304   // TH1F* rproj = 0;
0305   // TH1F* sproj = 0;
0306   //  TString *hnames_tmp = hnames;
0307 #ifdef DEBUG
0308   cout << "DEBUG: Creating canvas..." << endl;
0309 #endif
0310   TCanvas* canvas =  0;
0311   if (nhistos >4)
0312     canvas = new TCanvas(canvasName, canvasTitle, 1000, 1400);
0313   else
0314     canvas = new TCanvas(canvasName, canvasTitle, 1000, 1050);
0315   canvas->Divide(2,(nhistos+1)/2); //This way we print in 2 columns
0316 
0317 #ifdef DEBUG
0318   cout << "DEBUG: Looping over histograms" << endl;
0319 #endif
0320   for (unsigned int i = 0; i < nhistos; i++) {
0321 #ifdef DEBUG
0322     cout << "DEBUG: [" << i << "] histogram name: " << flush << hnames[i] << endl;
0323     cout << "DEBUG: rdir:  " << rdir << endl;
0324     cout << "DEBUG: sdir:  " << sdir << endl;
0325     cout << "DEBUG: rcollname: " << rcollname << endl;
0326     cout << "DEBUG: scollname: " << scollname << endl;
0327     cout << "DEBUG: rname: " << (rcollname + "/" + hnames[i]) << endl;
0328     cout << "DEBUG: sname: " << (scollname + "/" + hnames[i]) << endl;
0329 #endif
0330 
0331     TString hnamesi = hnames[i];
0332     
0333     // Skip histogram if no name is provided
0334     if (hnamesi == "") continue;
0335     
0336     // Get Histograms
0337     // + Reference
0338     TString hnames_tmp = hnamesi;
0339     rdir->cd(rcollname);
0340     TIter next(gDirectory->GetListOfKeys());
0341     TKey *key;
0342     while ((key = (TKey*)next())) { 
0343     
0344       TObject *obj = key->ReadObj();
0345       TString temp = obj->GetName();
0346       if ( (temp.Contains("nhits_vs_eta_pfx")) &&
0347        hnamesi.Contains("hits_eta")
0348      ) {
0349     hnames_tmp = temp;      
0350       }
0351       if ( (temp.Contains("chi2_vs_eta_pfx")) &&
0352        hnamesi.Contains("chi2mean") 
0353      ) {
0354     hnames_tmp = temp;      
0355       }
0356     } 
0357     
0358 
0359     
0360 #ifdef DEBUG
0361     cout << "DEBUG: Getting object reference samples " << (rcollname + "/" + hnames_tmp) << " from " << rdir << endl;
0362 #endif
0363     rdir->GetObject(rcollname + "/" + hnames_tmp, rh);
0364     if (! rh) {
0365       cout << "WARNING: Could not find a reference histogram or profile named " << hnames_tmp
0366        << " in " << rdir->GetName() << endl;
0367       cout << "         Skipping" << endl;
0368       continue;
0369     }
0370 
0371     //If it is a 2D project it in Y... is this what we always want?
0372     if (TString(rh->IsA()->GetName()) == "TH2F") {
0373 #ifdef DEBUG
0374       cout << "DEBUG: It is a TH2F object... project in Y!" << endl;
0375 #endif
0376       TH1* proj = ((TH2F*) rh)->ProjectionY();
0377       rh = proj;
0378     }
0379 
0380 
0381     // + Signal
0382     
0383     hnames_tmp=hnamesi;
0384     
0385     sdir->cd(scollname);
0386     TIter next2(gDirectory->GetListOfKeys());
0387     TKey *key2;
0388     while ((key2 = (TKey*)next2())) { 
0389     
0390       TObject *obj = key2->ReadObj();
0391       TString temp = obj->GetName();
0392       if ( (temp.Contains("nhits_vs_eta_pfx")) &&
0393        hnamesi.Contains("hits_eta")
0394      ) {
0395     hnames_tmp = temp;      
0396       }
0397       if ( (temp.Contains("chi2_vs_eta_pfx")) &&
0398        hnamesi.Contains("chi2mean")
0399      ) {
0400     hnames_tmp = temp;      
0401       }
0402     } 
0403 
0404 #ifdef DEBUG
0405     cout << "DEBUG: Getting object for selected sample " << (scollname + "/" + hnames_tmp) 
0406      << " from " << sdir << endl;
0407 #endif
0408     sdir->GetObject(scollname + "/" + hnames_tmp, sh);
0409     if (! sh) {
0410       cout << "WARNING: Could not find a signal histogram or profile named " << hnames_tmp
0411        << " in " << rdir->GetName() << endl;
0412       cout << "         Skipping" << endl;
0413       continue;
0414     }
0415 
0416     //If it is a 2D project it in Y... is this what we always want?
0417     if (TString(sh->IsA()->GetName()) == "TH2F") {
0418 #ifdef DEBUG
0419       cout << "DEBUG: " << hnames_tmp << " is a TH2F object... project in Y!" << endl;
0420 #endif
0421       TH1* proj = ((TH2F*) sh)->ProjectionY();
0422       sh = proj;
0423     }
0424 
0425     if(TString(sh->GetName()).Contains(" vs "))norm[i]= -999.;
0426 
0427 
0428     // Normalize
0429     if (norm[i] == -1.)
0430       NormalizeHistogramsTo1(rh, sh);
0431     else if (norm[i] == 0.)
0432       NormalizeHistogramsToFirst(rh,sh);
0433     else if (norm[i] == -999.){
0434       cout << "DEBUG: Normalizing histograms to nothing" << "..." << endl;
0435     }
0436     /*    else {
0437 #ifdef DEBUG
0438       cout << "DEBUG: Normalizing histograms to " << norm[i] << "..." << endl;
0439 #endif
0440       sh->Scale(norm[i]);
0441       rh->Scale(norm[i]);
0442       }*/
0443 
0444 
0445     // Set styles
0446 #ifdef DEBUG
0447     cout << "DEBUG: Setting styles..." << endl;
0448 #endif
0449     SetHistogramStyle(rh, 21, 4);
0450     SetHistogramStyle(sh, 20, 2);
0451     //Change titles
0452     if (htitles) {
0453       rh->SetTitle(htitles[i]);
0454       sh->SetTitle(htitles[i]);
0455       rh->GetYaxis()->SetTitle(htitles[i]);
0456       sh->GetYaxis()->SetTitle(htitles[i]);
0457     }
0458     //Change x axis range
0459     if (minx) {
0460 
0461       //      if (minx < -1E99) {
0462 
0463         rh->GetXaxis()->SetRangeUser(minx[i],maxx[i]);
0464     sh->GetXaxis()->SetRangeUser(minx[i],maxx[i]);
0465 
0466     //}
0467     }
0468     //Change y axis range
0469     if (miny) {
0470       if (miny[i] < -1E99) {
0471     rh->GetYaxis()->SetRangeUser(miny[i],maxy[i]);
0472     sh->GetYaxis()->SetRangeUser(miny[i],maxy[i]);
0473       }
0474     }
0475 
0476 
0477 
0478     // Move to subpad
0479     canvas->cd(i+1);
0480     
0481     TPad* pad1 = NULL;
0482     TPad* pad2 = NULL;
0483 
0484     pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
0485     pad2 = new TPad("pad2", "pad2", 0, 0.0, 1, 0.3);
0486 
0487     pad1->SetTopMargin   (0.08);
0488     pad1->SetBottomMargin(0.01);
0489     pad1->Draw();
0490 
0491     pad2->SetTopMargin   (0.05);
0492     pad2->SetBottomMargin(0.45);
0493     pad2->Draw();
0494 
0495     pad1->cd();
0496 
0497     // Check Logy
0498     if (logy) {
0499       if (logy[i])
0500     gPad->SetLogy();
0501     
0502     }
0503 
0504     // Set stat boxes
0505 #ifdef DEBUG
0506     cout << "DEBUG: Setting statistics..." << endl;
0507 #endif
0508     setStats(sh, rh, -1, 0, false);
0509     
0510     // Draw histogram
0511 #ifdef DEBUG
0512     cout << "DEBUG: Drawing histograms..." << endl;
0513 #endif
0514     if (sh->GetMaximum() > rh->GetMaximum()) {
0515       sh->Draw();
0516       rh->Draw("sames");
0517     }
0518     else {
0519       rh->Draw();
0520       sh->Draw("sames");
0521     }
0522     
0523     // Perform Kolmogorov test if needed
0524     if (doKolmo) {
0525 #ifdef DEBUG
0526       cout << "DEBUG: Performing Kolmogorov test..." << endl;
0527 #endif
0528       if (doKolmo[i]) {
0529     //  TPad* c1_1 = canvas->GetPad(i+1);
0530     double kstest = KolmogorovTest(sh,rh);
0531     if(kstest<0.7)
0532       gPad->SetFillColor(kBlue-10);
0533       }
0534     }
0535 
0536     pad2->cd();
0537 
0538     TH1* ratioplot = PlotRatiosHistograms(rh, sh);
0539     SetHistogramStyle(ratioplot, 21, 4);                                                   
0540     ratioplot->Draw("ep");
0541   } // End loop
0542 
0543 
0544    // Draw Legend
0545 #ifdef DEBUG
0546   cout << "DEBUG: Drawing legend..." << endl;
0547 #endif
0548   canvas->cd();
0549   
0550   TLegend* l = 0;
0551   if (nhistos > 4)
0552     l = new TLegend(0.20,0.665,0.80,0.685);
0553   else
0554     l = new TLegend(0.20,0.52,0.80,0.55);
0555 
0556   l->SetTextSize(0.011);
0557   l->SetLineColor(1);
0558   l->SetLineWidth(1);
0559   l->SetLineStyle(1);
0560   l->SetFillColor(0);
0561   l->SetBorderSize(2);
0562   l->AddEntry(rh,refLabel,"LPF");
0563   l->AddEntry(sh,newLabel,"LPF");
0564   l->Draw();
0565   
0566   // Print Canvas
0567   canvas->SaveAs(pdfFile+".pdf");
0568   canvas->SaveAs(pdfFile+".png");
0569 
0570   // Clean memory
0571   // delete l;
0572   delete canvas;
0573   cout << "     ... plotted histograms for " << canvasTitle << endl;
0574 }
0575 //
0576 ////////////////////////////////////////////////////////////
0577 
0578 
0579 ////////////////////////////////////////////////////////////
0580 //
0581 // Plot a page with 4 histograms
0582 //
0583 void Plot4Histograms(const TString& pdfFile,
0584              TDirectory* rdir, TDirectory* sdir,
0585              const TString& rcollname, const TString& scollname,
0586              const char* canvasName, const char* canvasTitle,
0587              const TString& refLabel, const TString& newLabel,
0588              const char** hnames, const char** htitles = 0,
0589              bool* logy = 0, bool* doKolmo = 0, Double_t* norm = 0,bool *resol = 0,
0590              Double_t* minx = 0, Double_t* maxx = 0,
0591              Double_t* miny = 0, Double_t* maxy = 0) {
0592   PlotNHistograms(pdfFile,
0593           rdir, sdir,
0594           rcollname, scollname,
0595           canvasName, canvasTitle,
0596           refLabel, newLabel,
0597           4, hnames, htitles,
0598           logy, doKolmo, norm,resol,minx,maxx,miny,maxy);
0599 }
0600 //
0601 ////////////////////////////////////////////////////////////
0602 
0603 
0604 ////////////////////////////////////////////////////////////
0605 //
0606 // Plot a page with 6 histograms
0607 //
0608 void Plot6Histograms(const TString& pdfFile,
0609              TDirectory* rdir, TDirectory* sdir,
0610              const TString& rcollname, const TString& scollname,
0611              const char* canvasName, const char* canvasTitle,
0612              const TString& refLabel, const TString& newLabel,
0613              const char** hnames, const char** htitles = 0,
0614              bool* logy = 0, bool* doKolmo = 0, Double_t* norm = 0,bool *resol=0,
0615              Double_t* minx = 0, Double_t* maxx = 0,
0616              Double_t* miny = 0, Double_t* maxy = 0) {
0617   
0618 
0619   PlotNHistograms(pdfFile,
0620           rdir, sdir,
0621           rcollname, scollname,
0622           canvasName, canvasTitle,
0623           refLabel, newLabel,
0624           6, hnames, htitles,
0625           logy, doKolmo, norm,resol,
0626           minx, maxx,
0627           miny, maxy);
0628 }
0629 //
0630 ////////////////////////////////////////////////////////////
0631 
0632 void Plot5Histograms(const TString& pdfFile,
0633              TDirectory* rdir, TDirectory* sdir,
0634              const TString& rcollname, const TString& scollname,
0635              const char* canvasName, const char* canvasTitle,
0636              const TString& refLabel, const TString& newLabel,
0637              const char** hnames, const char** htitles = 0,
0638              bool* logy = 0, bool* doKolmo = 0, Double_t* norm = 0,bool *resol=0,
0639              Double_t* minx = 0, Double_t* maxx = 0,
0640              Double_t* miny = 0, Double_t* maxy = 0) {
0641   
0642 
0643   PlotNHistograms(pdfFile,
0644           rdir, sdir,
0645           rcollname, scollname,
0646           canvasName, canvasTitle,
0647           refLabel, newLabel,
0648           5, hnames, htitles,
0649           logy, doKolmo, norm,resol,
0650           minx, maxx,
0651           miny, maxy);
0652 }
0653 //
0654 ////////////////////////////////////////////////////////////
0655 
0656 
0657 ////////////////////////////////////////////////////////////
0658 //
0659 // Normalize the two histograms to the entries on the first histogram
0660 //
0661 void NormalizeHistogramsToFirst(TH1* h1, TH1* h2) {
0662   if (h1==0 || h2==0) return;
0663   
0664   if ( h1->Integral() > 0 && h2->Integral() > 0 ){
0665     Double_t scale2 = h1->Integral()/h2->Integral();
0666     h2->Scale(scale2);
0667   }
0668 }
0669 //
0670 ////////////////////////////////////////////////////////////
0671 
0672 
0673 ////////////////////////////////////////////////////////////
0674 //
0675 // Normalize the two histograms to unity
0676 //
0677 void NormalizeHistogramsTo1(TH1* h1, TH1* h2) {
0678   if (!h1 || !h2) return;
0679   
0680   if ( h1->Integral() != 0 && h2->Integral() != 0 ) {
0681     Double_t scale1 = 1.0/h1->Integral();
0682     Double_t scale2 = 1.0/h2->Integral();
0683     h1->Scale(scale1);
0684     h2->Scale(scale2);
0685   }
0686 }
0687 //
0688 ////////////////////////////////////////////////////////////
0689 
0690 /*
0691 void NormalizeHistograms(TH1F* h1, Double_t nrm)
0692 {
0693   if (h1==0) return;
0694   h1->Scale(nrm);
0695 }
0696 
0697 
0698 void plotHistoInCanvas(TCanvas* canvas, unsigned int ncanvas,
0699             TH1F* s, TH1F* r,
0700             Size_t markersize,
0701             bool logy = false,
0702             bool doKolmo = false) {
0703   if ((s==0) || (r==0)) {
0704     cerr << "ERROR: Histograms not found!" << endl;
0705     return;
0706   }
0707 
0708   s->SetMarkerStyle(20);
0709   r->SetMarkerStyle(21);
0710   s->SetMarkerColor(2);
0711   r->SetMarkerColor(4);
0712   s->SetMarkerSize(markersize);
0713   r->SetMarkerSize(markersize);
0714   s->SetLineColor(2);
0715   r->SetLineColor(4);
0716   s->SetLineWidth(2);
0717   r->SetLineWidth(2);
0718 
0719   //setStats(r,s, startingY, startingX, fit);
0720   canvas->cd(ncanvas);
0721   if (logy)
0722     gPad->SetLogy();
0723   setStats(s,r, -1, 0, false);
0724   s->Draw();
0725   r->Draw("sames");
0726   if (doKolmo) {
0727     TPad *c1_1 = canvas->GetPad(ncanvas);
0728     double kstest = KolmogorovTest(s,r);
0729     if(kstest<0.7){
0730       c1_1->SetFillColor(kYellow);}
0731   }
0732 }
0733 
0734 
0735 void plotProfileInCanvas(TCanvas* canvas, unsigned int ncanvas,
0736              TProfile* s, TProfile* r,
0737              Size_t markersize) {
0738   if ((s==0) || (r==0)) {
0739     cerr << "ERROR: TProfile not found!" << endl;
0740     return;
0741   }
0742 
0743   s->SetMarkerStyle(20);
0744   r->SetMarkerStyle(21);
0745   s->SetMarkerColor(2);
0746   r->SetMarkerColor(4);
0747   s->SetMarkerSize(markersize);
0748   r->SetMarkerSize(markersize);
0749   s->SetLineColor(2);
0750   r->SetLineColor(4);
0751   s->SetLineWidth(2);
0752   r->SetLineWidth(2);
0753 
0754   //setStats(r,s, startingY, startingX, fit);
0755   canvas->cd(ncanvas);
0756   setStats(s,r, -1, 0, false);
0757   s->Draw();
0758   r->Draw("sames");
0759   //  double kstest = KolmogorovTest(s1,r1);
0760 }
0761 
0762 void plot4histos(TCanvas *canvas, 
0763          TH1F *s1,TH1F *r1, TH1F *s2,TH1F *r2, 
0764          TH1F *s3,TH1F *r3, TH1F *s4,TH1F *r4,
0765          TText* te,
0766          char * option, 
0767          double startingY, double startingX = .1,
0768          bool fit = false){
0769 #ifdef DEBUG
0770   cout << "DEBUG: plot4histos for histos" << endl;
0771 #endif
0772 
0773   canvas->Divide(2,2);
0774   plotHistoInCanvas(canvas, 1, s1, r1, 0.7, false, true);
0775   plotHistoInCanvas(canvas, 2, s2, r2, 0.1, false, true);
0776   plotHistoInCanvas(canvas, 3, s3, r3, 0.7, false, true);
0777   plotHistoInCanvas(canvas, 4, s4, r4, 0.7, false, true);
0778 }
0779 */
0780 
0781 /*
0782 void plot4histos(TCanvas *canvas,
0783          TProfile *s1,TProfile *r1, TProfile *s2,TProfile *r2,
0784          TProfile *s3,TProfile *r3, TProfile *s4,TProfile *r4,
0785          TText* te,
0786          char * option, 
0787          double startingY, double startingX = .1,
0788          bool fit = false){
0789 
0790   cout << "DEBUG: plot4histos for profiles" << endl;
0791 
0792   canvas->Divide(2,2);
0793   
0794   plotProfileInCanvas(canvas, 1, s1, r1, 0.7);
0795   plotProfileInCanvas(canvas, 2, s2, r2, 0.1);
0796   plotProfileInCanvas(canvas, 3, s3, r3, 0.7);
0797   plotProfileInCanvas(canvas, 4, s4, r4, 0.7);
0798   
0799 }
0800 
0801 void plot6histos(TCanvas *canvas, 
0802          TH1F *s1,TH1F *r1, TH1F *s2,TH1F *r2, 
0803          TH1F *s3,TH1F *r3, TH1F *s4,TH1F *r4,
0804          TH1F *s5,TH1F *r5, TH1F *s6,TH1F *r6,
0805          TText* te,
0806          char * option, 
0807          double startingY, double startingX = .1,
0808          bool fit = false){
0809   canvas->Divide(2,3);
0810 
0811   plotHistoInCanvas(canvas, 1, s1, r1, 0.7, false, true);
0812   plotHistoInCanvas(canvas, 2, s2, r2, 0.1, true,  true);
0813   plotHistoInCanvas(canvas, 3, s3, r3, 0.7, false, true);
0814   plotHistoInCanvas(canvas, 4, s4, r4, 0.7, true,  true);
0815   plotHistoInCanvas(canvas, 5, s5, r5, 0.7, false, true);
0816   plotHistoInCanvas(canvas, 6, s6, r6, 0.7, true,  true);
0817 
0818 }
0819 */
0820 ////////////////////////////////////////////////////////
0821 //
0822 // ratio plot from the two histograms 
0823 //
0824 /////////////////////////////////////////////////////////
0825 
0826 TH1* PlotRatiosHistograms(TH1* h1, TH1* h2){
0827 
0828   ++ratioCounter;
0829 
0830   Int_t nbinsx = h1->GetNbinsX();
0831   Double_t xbins[nbinsx+1];
0832 
0833   //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)                                             
0834   for (Int_t ibin=1; ibin<=nbinsx+1; ibin++) {
0835     Float_t xmin = h1->GetBinLowEdge(ibin);
0836     xbins[ibin-1] = xmin;
0837   }
0838 
0839   TH1F* h_ratio = new TH1F(Form("h_ratio_%d", ratioCounter), "", nbinsx, xbins);
0840 
0841   for (Int_t ibin=1; ibin<=nbinsx; ibin++) {
0842 
0843     Float_t h1Value = h1->GetBinContent(ibin);
0844     Float_t h1Error = h1->GetBinError(ibin);
0845 
0846     Float_t h2Value = h2->GetBinContent(ibin);
0847     Float_t h2Error = h2->GetBinError(ibin);
0848 
0849     Float_t ratioVal = 999;
0850     Float_t ratioErr = 999;
0851 
0852     if (h2Value > 0 || h2Value < 0) {
0853       ratioVal = h1Value / h2Value;
0854       ratioErr = h1Error / h2Value;
0855     }
0856 
0857     h_ratio->SetBinContent(ibin, ratioVal);
0858     h_ratio->SetBinError  (ibin, ratioErr);
0859 
0860   }
0861 
0862   h_ratio->SetTitle("");
0863   h_ratio->GetYaxis()->SetTitle("");
0864   h_ratio->GetYaxis()->SetRangeUser(0.4, 1.6);
0865 
0866   return h_ratio;
0867 
0868 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!