Macros/RooFit/CompareResol.cc

0001 #include "TCanvas.h"
0002 #include "TFile.h"
0003 #include "TH1.h"
0004 #include "TProfile.h"
0005 #include "TString.h"
0006 #include "TROOT.h"
0007
0008 #include "FitResolSlices.cc"
0009 #include "FitMassSlices.cc"
0010 #include "TLegend.h"
0011
0012 TH1 * getHisto(TDirectory * dir)
0013 {
0014   TCanvas * canvas = (TCanvas*)dir->Get("sigmaCanvas");
0015   return (TH1*)canvas->GetPrimitive("sigmaHisto");
0016 }
0017
0018 void drawCanvas(const TString & canvasName, TH1 * histo1, TH1 * histo2, const TString & histo1LegendTitle, const TString & histo2LegendTitle)
0019 {
0020   TCanvas * canvas = new TCanvas(canvasName, canvasName, 1000, 800);
0021   canvas->Draw();
0022   canvas->cd();
0023   histo1->Draw();
0024   histo1->SetMarkerStyle(24);
0025   histo1->SetMarkerSize(0.5);
0026
0027   if( histo2 != 0 ) {
0028     histo2->SetLineColor(kRed);
0029     histo2->SetMarkerColor(kRed);
0030     histo2->SetMarkerSize(0.5);
0031     histo2->SetMarkerStyle(24);
0032     histo2->Draw("same");
0033
0034     TLegend * legend = new TLegend(0.1, 0.7, 0.48, 0.9);
0035     legend->SetFillColor(0);
0036     legend->SetTextColor(1);
0037     legend->SetTextSize(0.02);
0038     legend->SetBorderSize(1);
0039     legend->AddEntry(histo1, histo1LegendTitle);
0040     legend->AddEntry(histo2, histo2LegendTitle);
0041     legend->Draw("same");
0042   }
0043
0044   canvas->Write();
0045 }
0046
0047 struct ConfCompare
0048 {
0049   ConfCompare() :
0050     file(0),
0051     takeProfile(false),
0052     MC(false)
0053   {}
0054
0055   TString inputFile;
0056   TString dirName;
0057   TString funcHistoName;
0058   TFile * file;
0059   TString xAxisTitle;
0060   TString yAxisTitle;
0061   TString mainDirName;
0062   TString subDirName;
0063   TString funcLegendTitle;
0064   TString histoLegendTitle;
0065   bool takeProfile;
0066   bool MC;
0067
0068   void compare()
0069   {
0070     gDirectory->mkdir(dirName);
0071     gDirectory->cd(dirName);
0072     TDirectory * outputDir = gDirectory;
0073
0074     // Resolutions from fitted function
0075     TFile * mainFile = new TFile(inputFile, "READ");
0076     TProfile * funcHisto = 0;
0077     if( takeProfile ) {
0078       TH2 * func2Dhisto = (TH2*)mainFile->FindObjectAny(funcHistoName);
0079       funcHisto = func2Dhisto->ProfileX();
0080     }
0081     else {
0082       funcHisto = (TProfile*)mainFile->FindObjectAny(funcHistoName);
0083     }
0084     outputDir->cd();
0085
0086     TH1 * histo = 0;
0087     // In this case we take the profile histogram
0088     if( MC || takeProfile ) {
0089       TDirectory* dir = (TDirectory*)file->Get(mainDirName);
0090       TDirectory * subDir = (TDirectory*)dir->Get(subDirName);
0091       histo = getHisto(subDir);
0092     }
0093
0094     funcHisto->GetXaxis()->SetTitle(xAxisTitle);
0095     funcHisto->GetYaxis()->SetTitle(yAxisTitle);
0096     funcHisto->GetYaxis()->SetTitleOffset(1.25);
0097
0098     drawCanvas(inputFile+dirName+funcHistoName+mainDirName+"_canvas", funcHisto, histo, funcLegendTitle, histoLegendTitle);
0099
0100     gDirectory->GetMotherDir()->cd();
0101   }
0102 };
0103
0104 class CompareResol
0105 {
0106 public:
0107   CompareResol() : file1_(0), file2_(0)
0108   {
0109     gROOT->SetStyle("Plain");
0110
0111     MC_ = true;
0112
0113     TString fileNum1("0");
0114     TString fileNum2("3");
0115
0116     TString inputFileName("_MuScleFit.root");
0117     TString outputFileName("ResolCheck_");
0118
0119     TString inputFile1(fileNum1+inputFileName);
0120     TString inputFile2(fileNum2+inputFileName);
0121     TString outputFileName1(outputFileName+fileNum1+".root");
0122     TString outputFileName2(outputFileName+fileNum2+".root");
0123
0124     TFile * outputFile1 = new TFile(outputFileName1, "RECREATE");
0125     TFile * outputFile2 = new TFile(outputFileName2, "RECREATE");
0126
0127     // Resolutions from MC comparison
0128     if( MC_ ) {
0129       FitResolSlices resolFitter1;
0130       FitResolSlices resolFitter2;
0131
0132       // muon Pt resolution vs muon Pt and Eta
0133       resolFitter1.fit(fileNum1+"_MuScleFit.root", "ResolCheck_"+fileNum1+".root", "gaussian",
0134                0., -0.1, 0.1, 0.03, 0., 0.1, "hResolPtGenVSMu_ResoVS", "ResolPtVs", outputFile1);
0135       if( fileNum2 != "" ) {
0136     resolFitter2.fit(fileNum2+"_MuScleFit.root", "ResolCheck_"+fileNum2+".root", "gaussian",
0137              0., -0.1, 0.1, 0.03, 0., 0.1, "hResolPtGenVSMu_ResoVS", "ResolPtVs", outputFile2);
0138       }
0139
0140       // Resonance mass resolution vs muon Pt and Eta
0141       resolFitter1.fit(fileNum1+"_MuScleFit.root", "ResolCheck_"+fileNum1+".root", "gaussian",
0142                0., -0.1, 0.1, 0.03, 0., 0.1, "DeltaMassOverGenMassVs", "ResolMassVsMuon", outputFile1);
0143       if( fileNum2 != "" ) {
0144     resolFitter2.fit(fileNum2+"_MuScleFit.root", "ResolCheck_"+fileNum2+".root", "gaussian",
0145              0., -0.1, 0.1, 0.03, 0., 0.1, "DeltaMassOverGenMassVs", "ResolMassVsMuon", outputFile2);
0146       }
0147     }
0148
0149     // Resolutions for mass from mass fits
0150     FitMassSlices massFitter1;
0151     FitMassSlices massFitter2;
0152
0153     TDirectory * massDir1 = outputFile1->mkdir("MassResol");
0154     massFitter1.rebinX = 4;
0155     massFitter1.fit(fileNum1+"_MuScleFit.root", "ResolCheck_"+fileNum1+".root", "gaussian", "exponential",
0156                 3.1, 3., 3.2, 0.03, 0., 0.1, massDir1);
0157     if( fileNum2 != "" ) {
0158       TDirectory * massDir2 = outputFile2->mkdir("MassResol");
0159       massFitter2.rebinX = 4;
0160       massFitter2.fit(fileNum2+"_MuScleFit.root", "ResolCheck_"+fileNum2+".root", "gaussian", "exponential",
0161               3.1, 3., 3.2, 0.03, 0., 0.1, massDir2);
0162     }
0163
0164     outputFile1->Write();
0165     outputFile1->Close();
0166     outputFile2->Write();
0167     outputFile2->Close();
0168
0169     // Reading back the closed files to do comparisons
0170     file1_ = new TFile(outputFileName1, "READ");
0171     if( fileNum2 != "" ) {
0172       file2_ = new TFile(outputFileName2, "READ");
0173     }
0174
0175     TFile * outputFile = new TFile("CompareResol.root", "RECREATE");
0176     outputFile->cd();
0177
0178     // Mass resolution
0179     // ---------------
0180     if( MC_ ) {
0181       // Comparison of MC resolution before-after the fit
0182       compareBeforeAfter("DeltaMassOverGenMassVs", "ResolMassVsMuonPt", "muon pt (GeV)", "#sigmaM/M");
0183       compareBeforeAfter("DeltaMassOverGenMassVs", "ResolMassVsMuonEta", "muon #eta", "#sigmaM/M");
0184     }
0185     // Comparison of mass resolution from data before and after the fit
0186     compareBeforeAfter("MassResol", "MassVsPt", "muon pt (GeV)", "#sigmaM");
0187     compareBeforeAfter("MassResol", "MassVsEta", "muon #eta", "#sigmaM");
0188
0189     // Make mass comparisons for both files
0190     makeAllMassComparisons(inputFile1, fileNum1, file1_);
0191     if( file2_ != 0 ) {
0192       makeAllMassComparisons(inputFile2, fileNum2, file2_);
0193     }
0194
0195     // Muon resolution
0196     // ---------------
0197     if( MC_ ) {
0198       // Muon Pt resolution vs Pt and Eta compared to MC resolution
0199       compareBeforeAfter("hResolPtGenVSMu_ResoVS", "ResolPtVsPt", "muon pt (GeV)", "#sigmaM");
0200       compareBeforeAfter("hResolPtGenVSMu_ResoVS", "ResolPtVsEta", "muon pt (GeV)", "#sigmaM");
0201     }
0202
0203     // Make sigmaPt/Pt comparisons for both files
0204     makeAllPtComparisons(inputFile1, fileNum1, file1_);
0205     if( file2_ != 0 ) {
0206       makeAllPtComparisons(inputFile2, fileNum2, file2_);
0207     }
0208
0209     outputFile->Write();
0210     outputFile->Close();
0211   }
0212 protected:
0213   void compareBeforeAfter(const TString & mainDirName, const TString & subDirName,
0214               const TString & xAxisTitle, const TString & yAxisTitle)
0215   {
0216     gDirectory->mkdir(mainDirName);
0217     gDirectory->cd(mainDirName);
0218
0219     TDirectory* dir1 = (TDirectory*)file1_->Get(mainDirName);
0220     TDirectory * subDir1 = (TDirectory*)dir1->Get(subDirName);
0221     TH1 * histo1 = getHisto(subDir1);
0222     histo1->GetXaxis()->SetTitle(xAxisTitle);
0223     histo1->GetYaxis()->SetTitle(yAxisTitle);
0224     histo1->GetYaxis()->SetTitleOffset(1.25);
0225
0226     TH1 * histo2 = 0;
0227     if( file2_ != 0 ) {
0228       TDirectory* dir2 = (TDirectory*)file2_->Get(mainDirName);
0229       TDirectory * subDir2 = (TDirectory*)dir2->Get(subDirName);
0230       histo2 = getHisto(subDir2);
0231       histo2->GetXaxis()->SetTitle(xAxisTitle);
0232       histo2->GetYaxis()->SetTitle(yAxisTitle);
0233       histo2->GetYaxis()->SetTitleOffset(1.25);
0234     }
0235
0236     drawCanvas(mainDirName+subDirName+"_canvas", histo1, histo2, histo1->GetTitle(), histo2->GetTitle());
0237
0238     gDirectory->GetMotherDir()->cd();
0239   }
0240
0241   void makeAllMassComparisons(const TString & inputFile, const TString & fileNum, TFile * mainFile)
0242   {
0243     // Mass relative resolution Vs Pt compared to MC resolution
0244
0245     ConfCompare conf;
0246     conf.inputFile = inputFile;
0247     conf.dirName = "FunctionMassRelativeResol_"+fileNum;
0248     conf.funcHistoName = "hFunctionResolMassVSMu_ResoVSPt_prof";
0249     conf.file = mainFile;
0250     conf.xAxisTitle = "muon pt (GeV)";
0251     conf.yAxisTitle = "#sigmaM/M";
0252     conf.mainDirName = "DeltaMassOverGenMassVs";
0253     conf.subDirName = "ResolMassVsMuonPt";
0254     conf.funcLegendTitle = "Mass relative resolution vs muon Pt from fitted function";
0255     conf.histoLegendTitle = "Mass relative resolution vs muon Pt from MC";
0256     conf.MC = MC_;
0257     conf.compare();
0258
0259     // Mass relative resolution Vs Eta compared to MC resolution
0260     conf.funcHistoName = "hFunctionResolMassVSMu_ResoVSEta_prof";
0261     conf.xAxisTitle = "muon #eta";
0262     conf.subDirName = "ResolMassVsMuonEta";
0263     conf.funcLegendTitle = "Mass relative resolution vs muon #eta from fitted function";
0264     conf.histoLegendTitle = "Mass relative resolution vs muon #eta from MC";
0265     conf.compare();
0266
0267     // Mass resolution Vs Pt compared to resolution from sigma of gaussian fits on data
0268     conf.dirName = "FunctionMassResol_"+fileNum;
0269     conf.funcHistoName = "hResolMassVSMu_ResoVSPt";
0270     conf.xAxisTitle = "muon pt (GeV)";
0271     conf.yAxisTitle = "#sigmaM";
0272     conf.mainDirName = "MassResol";
0273     conf.subDirName = "MassVsPt";
0274     conf.funcLegendTitle = "Mass resolution vs muon Pt from fitted function";
0275     conf.histoLegendTitle = "Mass relative resolution vs muon Pt from mass fits";
0276     conf.takeProfile = true;
0277     conf.compare();
0278
0279     // Mass resolution Vs Eta compared to resolution from sigma of gaussian fits on data
0280     conf.funcHistoName = "hResolMassVSMu_ResoVSEta";
0281     conf.xAxisTitle = "muon #eta";
0282     conf.yAxisTitle = "#sigmaM";
0283     conf.subDirName = "MassVsEta";
0284     conf.funcLegendTitle = "Mass resolution vs muon #eta from fitted function";
0285     conf.histoLegendTitle = "Mass relative resolution vs muon #eta from mass fits";
0286     conf.compare();
0287   }
0288
0289
0290   void makeAllPtComparisons(const TString & inputFile, const TString & fileNum, TFile * mainFile)
0291   {
0292     // Pt relative resolution Vs Pt compared to MC resolution
0293     ConfCompare conf;
0294     conf.inputFile = inputFile;
0295     conf.dirName = "FunctionPtRelativeResol_"+fileNum;
0296     conf.funcHistoName = "hFunctionResolPt_ResoVSPt_prof";
0297     conf.file = mainFile;
0298     conf.xAxisTitle = "muon pt (GeV)";
0299     conf.yAxisTitle = "#sigmaPt/Pt";
0300     conf.mainDirName = "hResolPtGenVSMu_ResoVS";
0301     conf.subDirName = "ResolPtVsPt";
0302     conf.funcLegendTitle = "Muon Pt relative resolution vs muon Pt from fitted function";
0303     conf.histoLegendTitle = "Muon Pt relative resolution vs muon Pt from MC";
0304     conf.MC = MC_;
0305     conf.compare();
0306
0307     // Pt relative resolution Vs Eta compared to MC resolution
0308     conf.dirName = "FunctionEtaRelativeResol_"+fileNum;
0309     conf.funcHistoName = "hFunctionResolPt_ResoVSEta_prof";
0310     conf.xAxisTitle = "muon #eta";
0311     conf.subDirName = "ResolPtVsEta";
0312     conf.funcLegendTitle = "Muon Pt relative resolution vs muon #eta from fitted function";
0313     conf.histoLegendTitle = "Muon Pt relative resolution vs muon #eta from MC";
0314     conf.compare();
0315   }
0316
0317   RooPlot * fit(TH1 * histo, const TString & fileName, FitWithRooFit * fitter,
0318         const TString & fitType, const double & xMin, const double & xMax)
0319   {
0320     gDirectory->mkdir(fileName);
0321     gDirectory->cd(fileName);
0322
0323     fitter->fit(histo, "", fitType, xMin, xMax);
0324     RooPlot * rooPlot = (RooPlot*)gDirectory->Get(TString(histo->GetName())+"_frame");
0325
0326     gDirectory->GetMotherDir()->cd();
0327
0328     return rooPlot;
0329   }
0330
0331   TFile * file1_;
0332   TFile * file2_;
0333
0334   bool MC_;
0335 };