File indexing completed on 2024-04-06 12:22:42
0001 #include "TCanvas.h"
0002 #include "TFile.h"
0003 #include "TH1.h"
0004 #include "TString.h"
0005 #include "TROOT.h"
0006 #include "TLegend.h"
0007
0008
0009 #include "FitMassSlices.cc"
0010 #include "Legend.h"
0011
0012 class CompareBias
0013 {
0014 public:
0015 CompareBias() : file1_(0), file2_(0)
0016 {
0017 gROOT->SetStyle("Plain");
0018
0019 doFit_ = false;
0020
0021 TString fileNum1("0");
0022 TString fileNum2("2");
0023
0024 TString inputFileName("_MuScleFit.root");
0025 TString outputFileName("BiasCheck_");
0026
0027 TString inputFile1(fileNum1+inputFileName);
0028 TString inputFile2(fileNum2+inputFileName);
0029 TString outputFile1(outputFileName+fileNum1+".root");
0030 TString outputFile2(outputFileName+fileNum2+".root");
0031
0032 FitMassSlices fitter1;
0033 FitMassSlices fitter2;
0034
0035 fitter1.rebinZ = 1;
0036
0037 fitter1.useChi2 = false;
0038
0039 fitter1.rebinX = 2;
0040 fitter1.rebinY = 2;
0041 fitter1.sigma2 = 1.;
0042 fitter1.fit(fileNum1+"_MuScleFit.root", "BiasCheck_"+fileNum1+".root", "crystalBall", "exponential", 3.095, 2.8, 3.4, 0.04, 0., 0.);
0043
0044
0045 if( fileNum2 != "" ) {
0046 fitter2.rebinX = 2;
0047 fitter2.rebinY = 2;
0048 fitter2.sigma2 = 1.;
0049 fitter2.fit(fileNum2+"_MuScleFit.root", "BiasCheck_"+fileNum2+".root", "crystalBall", "exponential", 3.095, 2.8, 3.4, 0.04, 0., 0.);
0050
0051 }
0052
0053 file1_ = new TFile(outputFile1, "READ");
0054 if( fileNum2 != "" ) {
0055 file2_ = new TFile(outputFile2, "READ");
0056 }
0057
0058 TFile * outputFile = new TFile("CompareBias.root", "RECREATE");
0059 outputFile->cd();
0060
0061 compare("MassVsPt", "uniform", 1., 8., "muon pt (GeV)", "Mass (GeV)");
0062 compare("MassVsEta", "uniform", -2.8, 2.8, "muon #eta", "Mass (GeV)");
0063
0064
0065 compare("MassVsEtaPlus", "uniform", -2.8, 2.8, "muon + #eta", "Mass (GeV)");
0066 compare("MassVsEtaMinus", "uniform", -2.8, 2.8, "muon - #eta", "Mass (GeV)");
0067
0068
0069
0070
0071 compare("MassVsPhiPlus", "uniform", -3.14, 3.14, "muon(+) #phi", "Mass (GeV)");
0072 compare("MassVsPhiMinus", "uniform", -3.14, 3.14, "muon(-) #phi", "Mass (GeV)");
0073
0074
0075
0076
0077 outputFile->Write();
0078 outputFile->Close();
0079 }
0080 protected:
0081 void compare(const TString & histoName, const TString & fitType, const double & xMin, const double & xMax,
0082 const TString & xAxisTitle, const TString & yAxisTitle)
0083 {
0084 gDirectory->mkdir(histoName);
0085 gDirectory->cd(histoName);
0086
0087 TH1 * histo1 = getHisto(file1_, histoName);
0088 histo1->GetXaxis()->SetTitle(xAxisTitle);
0089 histo1->GetYaxis()->SetTitle(yAxisTitle);
0090 histo1->GetYaxis()->SetTitleOffset(1.25);
0091 TH1 * histo2 = 0;
0092 if( file2_ != 0 ) {
0093 histo2 = getHisto(file2_, histoName);
0094 histo2->GetXaxis()->SetTitle(xAxisTitle);
0095 histo2->GetYaxis()->SetTitle(yAxisTitle);
0096 histo2->GetYaxis()->SetTitleOffset(1.25);
0097 }
0098
0099
0100
0101
0102
0103
0104 if( doFit_ ) {
0105 fitWithRoot(histo1, histo2, xMin, xMax, fitType);
0106 }
0107 else {
0108 TCanvas * canvas = drawCanvas(histo1, histo2, true);
0109 canvas->Write();
0110 }
0111 gDirectory->GetMotherDir()->cd();
0112 }
0113
0114 TH1 * getHisto(TFile * file, const TString & histoName)
0115 {
0116 TDirectory* dir = (TDirectory*)file->Get(histoName);
0117 TCanvas * canvas = (TCanvas*)dir->Get("meanCanvas");
0118 return (TH1*)canvas->GetPrimitive("meanHisto");
0119 }
0120
0121 void fitWithRoot(TH1 * histo1, TH1 * histo2, const double & xMin, const double & xMax, const TString & fitType)
0122 {
0123 TF1 * f1 = 0;
0124 TF1 * f2 = 0;
0125 if( fitType == "uniform" ) {
0126 f1 = new TF1("uniform1", "pol0", xMin, xMax);
0127 if( file2_ != 0 ) {
0128 f2 = new TF1("uniform2", "pol0", xMin, xMax);
0129 }
0130 }
0131 else if( fitType == "sinusoidal" ) {
0132 f1 = new TF1("sinusoidal1", "[0] + [1]*sin([2]*x + [3])", xMin, xMax);
0133 f1->SetParameter(1, 2.);
0134 f1->SetParameter(2, 1.);
0135 f1->SetParameter(3, 1.);
0136 if( file2_ != 0 ) {
0137 f2 = new TF1("sinusoidal2", "[0] + [1]*sin([2]*x + [3])", xMin, xMax);
0138 f2->SetParameter(1, 2.);
0139 f2->SetParameter(2, 1.);
0140 f2->SetParameter(3, 1.);
0141 }
0142 }
0143 else {
0144 std::cout << "Wrong fit type: " << fitType << std::endl;
0145 exit(1);
0146 }
0147
0148 histo1->Fit(f1, "", "", xMin, xMax);
0149 if( histo2 != 0 ) {
0150 histo2->Fit(f2, "", "", xMin, xMax);
0151 }
0152
0153 TCanvas * canvas = drawCanvas(histo1, histo2);
0154
0155 f1->Draw("same");
0156 if( histo2 != 0 ) {
0157 f2->Draw("same");
0158 f2->SetLineColor(kRed);
0159 }
0160 TwinLegend legends;
0161 legends.setText(f1, f2);
0162 legends.Draw("same");
0163
0164 canvas->Write();
0165 }
0166
0167 TCanvas * drawCanvas(TH1 * histo1, TH1 * histo2, const bool addLegend = false)
0168 {
0169 TCanvas * canvas = new TCanvas(TString(histo1->GetName())+"_canvas", TString(histo1->GetName())+" canvas", 1000, 800);
0170 canvas->Draw();
0171 canvas->cd();
0172 histo1->Draw();
0173 histo1->SetMarkerStyle(24);
0174 histo1->SetMarkerSize(0.5);
0175
0176 if( histo2 != 0 ) {
0177 histo2->SetLineColor(kRed);
0178 histo2->SetMarkerColor(kRed);
0179 histo2->SetMarkerSize(0.5);
0180 histo2->SetMarkerStyle(24);
0181 histo2->Draw("same");
0182 if( addLegend ) {
0183 TLegend * leg = new TLegend(0.1,0.7,0.48,0.9);
0184 leg->AddEntry(histo1,"Before calibration","pl");
0185 leg->AddEntry(histo2,"After calibration","pl");
0186 leg->Draw("same");
0187 }
0188 }
0189
0190 return canvas;
0191 }
0192
0193 void fitWithRooFit(TH1 * histo1, TH1 * histo2, const TString & histoName,
0194 const TString & fitType, const double & xMin, const double & xMax)
0195 {
0196 FitWithRooFit fitter;
0197 fitter.initConstant(3.097, 3.05, 3.15);
0198
0199
0200 RooPlot * rooPlot1 = fit( histo1, file1_->GetName(), &fitter, fitType, xMin, xMax );
0201
0202 RooRealVar * constant = fitter.constant();
0203 std::cout << "fitted value for constant 1 = " << constant->getVal() << std::endl;
0204
0205 RooPlot * rooPlot2 = fit( histo2, file2_->GetName(), &fitter, fitType, xMin, xMax );
0206
0207 constant = fitter.constant();
0208 std::cout << "fitted value for constant 2 = " << constant->getVal() << std::endl;
0209
0210 TCanvas * canvas = new TCanvas(histoName+"_canvas", histoName+" canvas", 1000, 800);
0211 canvas->Draw();
0212 canvas->cd();
0213 rooPlot1->Draw();
0214 rooPlot2->SetLineColor(kRed);
0215 rooPlot2->SetMarkerColor(kRed);
0216 rooPlot2->Draw("same");
0217 canvas->Write();
0218 }
0219
0220 RooPlot * fit(TH1 * histo, const TString & fileName, FitWithRooFit * fitter,
0221 const TString & fitType, const double & xMin, const double & xMax)
0222 {
0223 gDirectory->mkdir(fileName);
0224 gDirectory->cd(fileName);
0225
0226
0227 fitter->fit(histo, "", fitType, xMin, xMax);
0228 RooPlot * rooPlot = (RooPlot*)gDirectory->Get(TString(histo->GetName())+"_frame");
0229
0230 gDirectory->GetMotherDir()->cd();
0231
0232 return rooPlot;
0233 }
0234
0235 TFile * file1_;
0236 TFile * file2_;
0237 bool doFit_;
0238 };