test/Macros/CompareRecoGenPt.C

0001 /**
0002  * This compiled macro compares the rms on the difference between reco and gen
0003  * pt for muons in two different files.
0004  * It prints the rms distributions superimposed.
0005  * It can be used to verify the effect of corrections.
0006  */
0007
0008 #include "TFile.h"
0009 #include "TDirectory.h"
0010 #include "TString.h"
0011 #include "TProfile.h"
0012 #include "TCanvas.h"
0013 #include "TLegend.h"
0014 #include "TF1.h"
0015 #include "TROOT.h"
0016 #include "TStyle.h"
0017
0018 #include <sstream>
0019 #include <iostream>
0020 #include <iomanip>
0021
0022 /// Helper function getting the histogram from file.
0023 TProfile *getHistogram(const TString &fileName) {
0024   TFile *file = new TFile(fileName, "READ");
0025   if (file == 0) {
0026     std::cout << "Wrong file: " << fileName << std::endl;
0027     exit(1);
0028   }
0029   TDirectory *dir = (TDirectory *)file->Get("hPtRecoVsPtGen");
0030   if (dir == 0) {
0031     std::cout << "Wrong directory for file: " << fileName << std::endl;
0032     exit(1);
0033   }
0034   TProfile *profile = (TProfile *)dir->Get("hPtRecoVsPtGenProf");
0035   if (profile == 0) {
0036     std::cout << "Wrong histogram for file: " << fileName << std::endl;
0037     exit(1);
0038   }
0039   return profile;
0040 }
0041
0042 /// Helper function building the histogram from the TProfile settings.
0043 TH1F *makeHistogram(const TProfile *profile, const TString &name) {
0044   return new TH1F(TString(profile->GetName()) + name,
0045                   TString(profile->GetTitle()) + " " + name,
0046                   profile->GetNbinsX(),
0047                   profile->GetXaxis()->GetXmin(),
0048                   profile->GetXaxis()->GetXmax());
0049 }
0050
0051 /// Helper function to write the histograms to file.
0052 void saveHistograms(TH1 *histo1, TH1 *histo2) {
0053   histo1->Draw();
0054   histo2->SetLineColor(kRed);
0055   histo2->Draw("Same");
0056   TLegend *leg = new TLegend(0.65, 0.85, 1, 1);
0057   leg->SetFillColor(0);
0058   leg->AddEntry(histo1, "before calibration", "L");
0059   leg->AddEntry(histo2, "after calibration", "L");
0060   leg->Draw("same");
0061 }
0062
0063 #include "TPaveText.h"
0064 /// Helper class holding a TPaveText for better formatting and predefined options
0065 class PaveText {
0066 public:
0067   PaveText(const double &textX = 0.7, const double &textY = 0.4) {
0068     paveText_ = new TPaveText(textX, textY, textX + 0.2, textY + 0.17, "NDC");
0069   }
0070   void AddText(const TString &text) { paveText_->AddText(text); }
0071   void Draw(const TString &option) {
0072     paveText_->SetFillColor(0);  // text is black on white
0073     paveText_->SetTextSize(0.03);
0074     paveText_->SetBorderSize(0);
0075     paveText_->SetTextAlign(12);
0076     paveText_->Draw(option);
0077   }
0078   void SetTextColor(const int color) { paveText_->SetTextColor(color); }
0079
0080 protected:
0081   TPaveText *paveText_;
0082 };
0083
0084 /**
0085  * Compute the precision to give to the stream operator so that the passed number
0086  * will be printed with two significant figures.
0087  */
0088 int precision(const double &value) {
0089   // Counter gives the precision
0090   int precision = 1;
0091   int k = 1;
0092   while (int(value * k) == 0) {
0093     k *= 10;
0094     ++precision;
0095   }
0096   return precision;
0097 }
0098
0099 /// Helper function to extract and format the text for the fitted parameters
0100 void getParameters(const TF1 *func, TString &fit1, TString &fit2, TString &fit3) {
0101   std::stringstream a;
0102
0103   double error = func->GetParError(0);
0104   a << std::setprecision(precision(error)) << std::fixed << func->GetParameter(0);
0105   fit1 += a.str() + "+-";
0106   a.str("");
0107   a << error;
0108   fit1 += a.str();
0109   a.str("");
0110
0111   error = func->GetParError(1);
0112
0113   a << std::setprecision(precision(error)) << std::fixed << func->GetParameter(1);
0114   fit2 += a.str() + "+-";
0115   a.str("");
0116   a << func->GetParError(1);
0117   fit2 += a.str();
0118   a.str("");
0119   a << std::setprecision(1) << std::fixed << func->GetChisquare();
0120   fit3 += a.str() + "/";
0121   a.str("");
0122   a << std::setprecision(0) << std::fixed << func->GetNDF();
0123   fit3 += a.str();
0124 }
0125
0126 void CompareRecoGenPt(const TString &fileNum1 = "0", const TString &fileNum2 = "1") {
0127   TFile *outputFile = new TFile("CompareRecoGenPt.root", "RECREATE");
0128
0129   TProfile *profile1 = getHistogram(fileNum1 + "_MuScleFit.root");
0130   profile1->SetXTitle("gen muon Pt (GeV)");
0131   profile1->SetYTitle("reco muon Pt (GeV)");
0132   TProfile *profile2 = getHistogram(fileNum2 + "_MuScleFit.root");
0133
0134   int xBins = profile1->GetNbinsX();
0135   if (xBins != profile2->GetNbinsX()) {
0136     std::cout << "Wrong number of bins" << std::endl;
0137     exit(1);
0138   }
0139
0140   // Loop on all bins and fill a histogram with the mean values.
0141   // Fill also a histogram with the rms values.
0142
0143   outputFile->cd();
0144
0145   TH1F *meanHisto1 = makeHistogram(profile1, "mean");
0146   TH1F *meanHisto2 = makeHistogram(profile2, "mean");
0147   TH1F *rmsHisto1 = makeHistogram(profile1, "rms");
0148   TH1F *rmsHisto2 = makeHistogram(profile2, "rms");
0149   for (int iBin = 1; iBin <= xBins; ++iBin) {
0150     //     if( profile1->GetBinError(iBin) != 0 ) {
0151     meanHisto1->SetBinContent(iBin, profile1->GetBinContent(iBin));
0152     meanHisto1->SetBinError(iBin, profile1->GetBinError(iBin));
0153     //     }
0154     //     if( profile2->GetBinError(iBin) ) {
0155     meanHisto2->SetBinContent(iBin, profile2->GetBinContent(iBin));
0156     meanHisto2->SetBinError(iBin, profile2->GetBinError(iBin));
0157     //     }
0158     rmsHisto1->SetBinContent(iBin, profile1->GetBinError(iBin));
0159     rmsHisto2->SetBinContent(iBin, profile2->GetBinError(iBin));
0160   }
0161
0162   // Setting all weigths to 1 ("W" option) because of Profile errors for low statistics bins biasing the fit
0163
0164   // meanHisto1->Fit("pol1", "W", "", 2, 1000);
0165   profile1->Fit("pol1", "W", "", 0, 1000);
0166   TF1 *func1 = profile1->GetFunction("pol1");
0167   // TF1 * func1 = meanHisto1->GetFunction("pol1");
0168   func1->SetLineWidth(1);
0169   func1->SetLineColor(kBlack);
0170
0171   profile2->Fit("pol1", "W", "", 0, 1000);
0172   // meanHisto2->Fit("pol1", "W", "", 2, 1000);
0173   TF1 *func2 = profile2->GetFunction("pol1");
0174   // TF1 * func2 = meanHisto2->GetFunction("pol1");
0175   func2->SetLineWidth(1);
0176   func2->SetLineColor(kRed);
0177
0178   TCanvas *canvas = new TCanvas("before", "before corrections", 1000, 800);
0179   // canvas->Divide(2);
0180   canvas->cd();
0181   // canvas->cd(1);
0182   // canvas->cd(2);
0183   // saveHistograms(rmsHisto1, rmsHisto2);
0184   // saveHistograms(meanHisto1, meanHisto2);
0185   saveHistograms(profile1, profile2);
0186   func1->Draw("same");
0187   func2->Draw("same");
0188
0189   TString fit11("a = ");
0190   TString fit12("b = ");
0191   TString fit13("#chi^2/ndf = ");
0192   getParameters(func1, fit11, fit12, fit13);
0193   PaveText pt1(0.45, 0.15);
0194   pt1.AddText("before:");
0195   pt1.AddText(fit11);
0196   pt1.AddText(fit12);
0197   pt1.AddText(fit13);
0198   pt1.Draw("same");
0199
0200   TString fit21("a = ");
0201   TString fit22("b = ");
0202   TString fit23("#chi^2/ndf = ");
0203   getParameters(func2, fit21, fit22, fit23);
0204   PaveText pt2(0.65, 0.15);
0205   pt2.SetTextColor(2);
0206   pt2.AddText("after:");
0207   pt2.AddText(fit21);
0208   pt2.AddText(fit22);
0209   pt2.AddText(fit23);
0210   pt2.Draw("same");
0211   gStyle->SetOptStat(0);
0212
0213   canvas->Write();
0214
0215   outputFile->Write();
0216   outputFile->Close();
0217
0218   //   TLegend *leg = new TLegend(0.2,0.4,0.4,0.6);
0219   //   leg->SetFillColor(0);
0220   //   leg->AddEntry(func1,"fit of before","L");
0221   //   leg->AddEntry(func2,"fit of after","L");
0222   //   leg->Draw("same");
0223 }