Back to home page

Project CMSSW displayed by LXR

 
 

    


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; // 2 for Z
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     // fitter1.fit(fileNum1+"_MuScleFit.root", "BiasCheck_"+fileNum1+".root", "voigtian", "", 91, 80, 100, 2, 0.1, 10);
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       // fitter2.fit(fileNum2+"_MuScleFit.root", "BiasCheck_"+fileNum2+".root", "voigtian", "", 91, 80, 100, 2, 0.1, 10); 
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     //r.c. ------------------
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     //    compare("MassVsEtaPhiPlus", "uniform", -2.8, 2.8, "muon - #eta #phi", "Mass (GeV)");   
0068     //r.c. -----------------
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     // compare("MassVsPhiPlus", "sinusoidal", -3.14, 3.14, "muon(+) #phi", "Mass (GeV)");
0075     // compare("MassVsPhiMinus", "sinusoidal", -3.14, 3.14, "muon(-) #phi", "Mass (GeV)");
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     // Fit using RooFit
0100     // The polynomial in RooFit is a pdf, so it is normalized to unity. This seems to give problems.
0101     // fitWithRooFit(histo1, histo2, histoName, fitType, xMin, xMax);
0102 
0103     // Fit with standard root, but then we also need to build the legends.
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     // fitter.initLinearTerm(0., -1., 1.);
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     // fitter->fit(histo, "", fitType, xMin, xMax, true);
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 };