Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:31:35

0001 #include "TCanvas.h"
0002 #include "TFile.h"
0003 #include "TH1.h"
0004 #include "TH2.h"
0005 #include "TString.h"
0006 #include "TROOT.h"
0007 #include "TLegend.h"
0008 
0009 #include "FitMassSlices.cc"
0010 #include "FitMass1D.cc"
0011 #include "Legend.h"
0012 
0013 
0014 class CompareBiasUpsilonValidation
0015 {
0016 public:
0017   CompareBiasUpsilonValidation(const int rebinXphi = 4, const int rebinXetadiff = 2, const int rebinXeta = 2, const int rebinXpt = 8)
0018   {
0019     gROOT->SetStyle("Plain");
0020 
0021     doFit_ = false;
0022 
0023 
0024     TString inputFileName("0_zmumuHisto.root");
0025     TString outputFileName("BiasCheck.root");
0026 
0027 
0028     FitMassSlices fitter;
0029 
0030     fitter.rebinX = 2; // for further rebinning for phi use rebinXphi in FitMassSlices.cc (L20)
0031     fitter.rebinY = 2; // default 2
0032     fitter.rebinZ = 1; // default 2
0033 
0034     fitter.useChi2 = false;
0035     fitter.sigma2 = 1.;
0036 
0037     double Mmin(9.1), Mmax(9.7);
0038     fitter.fit(
0039       inputFileName, outputFileName, "breitWignerTimesCB", "exponential", 9.46, Mmin, Mmax, 0.3, 0.001, 2.,
0040       rebinXphi, rebinXetadiff, rebinXeta, rebinXpt
0041       );
0042 
0043     FitMass1D fitMass1D;
0044     fitMass1D.fitter()->initMean(9.46, Mmin, Mmax);
0045     fitMass1D.fitter()->initGamma(5.4e-5, 0., 10.);
0046     fitMass1D.fitter()->gamma()->setConstant(kTRUE);
0047     fitMass1D.fitter()->initMean2(0., -20., 20.);
0048     fitMass1D.fitter()->mean2()->setConstant(kTRUE);
0049     fitMass1D.fitter()->initSigma(0.07, 0., 5.);
0050     fitMass1D.fitter()->initAlpha(1.5, 0.05, 10.);
0051     fitMass1D.fitter()->initN(1, 0.01, 100.);
0052     fitMass1D.fitter()->initExpCoeffA0(-1., -10., 10.);
0053     fitMass1D.fitter()->initExpCoeffA1(0., -10., 10.);
0054     fitMass1D.fitter()->initExpCoeffA2(0., -2., 2.);
0055     fitMass1D.fitter()->initFsig(0.9, 0., 1.);
0056     fitMass1D.fitter()->initA0(0., -10., 10.);
0057     fitMass1D.fitter()->initA1(0., -10., 10.);
0058     fitMass1D.fitter()->initA2(0., -10., 10.);
0059     fitMass1D.fitter()->initA3(0., -10., 10.);
0060     fitMass1D.fitter()->initA4(0., -10., 10.);
0061     fitMass1D.fitter()->initA5(0., -10., 10.);
0062     fitMass1D.fitter()->initA6(0., -10., 10.);
0063 
0064     /// Let's fit
0065     fitMass1D.fit(inputFileName, outputFileName, "UPDATE", Mmin, Mmax, "breitWignerTimesCB", "exponentialpol");
0066 
0067   }
0068 protected:
0069 
0070   void compare(const TString & histoName, const TString & fitType, const double & xMin, const double & xMax,
0071     const TString & xAxisTitle, const TString & yAxisTitle, const TString& leg)
0072   {
0073     gDirectory->mkdir(histoName);
0074     gDirectory->cd(histoName);
0075 
0076     TH1 * histo = (TH1*)getHisto(file_, histoName);
0077     histo->GetXaxis()->SetTitle(xAxisTitle);
0078     histo->GetYaxis()->SetTitle(yAxisTitle);
0079     histo->GetYaxis()->SetTitleOffset(1.25);
0080 
0081     // Fit using RooFit
0082     // The polynomial in RooFit is a pdf, so it is normalized to unity. This seems to give problems.
0083     // fitWithRooFit(histo, histoName, fitType, xMin, xMax);
0084     // Fit with standard root, but then we also need to build the legends.
0085     if (doFit_) {
0086       fitWithRoot(histo, xMin, xMax, fitType);
0087     }
0088     else {
0089       TCanvas * canvas = drawCanvas(histo, leg, true);
0090       canvas->Write();
0091     }
0092     gDirectory->GetMotherDir()->cd();
0093   }
0094 
0095   void compare(const TString & histoName,
0096     const double & xMin, const double & xMax, const TString & xAxisTitle,
0097     const double & yMin, const double & yMax, const TString & yAxisTitle,
0098     const TString & zAxisTitle, const TString& leg)
0099   {
0100     gDirectory->mkdir(histoName);
0101     gDirectory->cd(histoName);
0102 
0103     TH2 * histo = (TH2*)getHisto(file_, histoName);
0104     histo->GetXaxis()->SetTitle(xAxisTitle);
0105     histo->GetYaxis()->SetTitle(yAxisTitle);
0106     //histo->GetYaxis()->SetTitleOffset(1.25);
0107 
0108     // Fit using RooFit
0109     // The polynomial in RooFit is a pdf, so it is normalized to unity. This seems to give problems.
0110     // fitWithRooFit(histo, histo2, histoName, fitType, xMin, xMax);
0111     // Fit with standard root, but then we also need to build the legends.
0112 
0113     TCanvas * canvas = drawCanvas(histo);
0114     canvas->Write();
0115 
0116     gDirectory->GetMotherDir()->cd();
0117   }
0118 
0119 
0120   TH1* getHisto(TFile * file, const TString & histoName)
0121   {
0122     TDirectory* dir = (TDirectory*)file->Get(histoName);
0123     TCanvas * canvas = (TCanvas*)dir->Get("meanCanvas");
0124     return (TH1*)canvas->GetPrimitive("meanHisto");
0125   }
0126 
0127   void fitWithRoot(TH1 * histo, const double & xMin, const double & xMax, const TString & fitType)
0128   {
0129     TF1 * f1 = 0;
0130     if (fitType == "uniform") {
0131       f1 = new TF1("uniform1", "pol0", xMin, xMax);
0132     }
0133     else if (fitType == "sinusoidal") {
0134       f1 = new TF1("sinusoidal1", "[0] + [1]*sin([2]*x + [3])", xMin, xMax);
0135       f1->SetParameter(1, 2.);
0136       f1->SetParameter(2, 1.);
0137       f1->SetParameter(3, 1.);
0138     }
0139     else {
0140       std::cout << "Wrong fit type: " << fitType << std::endl;
0141       exit(1);
0142     }
0143 
0144     histo->Fit(f1, "", "", xMin, xMax);
0145 
0146     TCanvas * canvas = drawCanvas(histo);
0147 
0148     f1->Draw("same");
0149     //TLegend legend;
0150     //legend.setText(f1);
0151     //legend.Draw("same");
0152 
0153     canvas->Write();
0154   }
0155 
0156   TCanvas * drawCanvas(TH1 * histo, TString legText="geo 1", const bool addLegend = false)
0157   {
0158     TCanvas * canvas = new TCanvas(TString(histo->GetName())+"_canvas", TString(histo->GetName())+" canvas", 1000, 800);
0159     canvas->Draw();
0160     canvas->cd();
0161     histo->Draw();
0162     histo->SetMarkerStyle(24);
0163     histo->SetMarkerSize(0.5);
0164 
0165 
0166     if (addLegend) {
0167       TLegend * leg = new TLegend(0.1, 0.7, 0.48, 0.9);
0168       leg->AddEntry(histo, legText, "pl");
0169       leg->Draw("same");
0170     }
0171 
0172     return canvas;
0173   }
0174 
0175   void fitWithRooFit(TH1 * histo, const TString & histoName,
0176     const TString & fitType, const double & xMin, const double & xMax)
0177   {
0178     FitWithRooFit fitter;
0179     fitter.initA0(3.097, 3.05, 3.15);
0180     // fitter.initLinearTerm(0., -1., 1.);
0181 
0182     RooPlot * rooPlot1 = fit(histo, file_->GetName(), &fitter, fitType, xMin, xMax);
0183 
0184     RooRealVar * constant = fitter.a0();
0185     std::cout << "fitted value for constant 1 = " << constant->getVal() << std::endl;
0186 
0187     TCanvas * canvas = new TCanvas(histoName+"_canvas", histoName+" canvas", 1000, 800);
0188     canvas->Draw();
0189     canvas->cd();
0190     rooPlot1->Draw();
0191     canvas->Write();
0192   }
0193 
0194   RooPlot * fit(TH1 * histo, const TString & fileName, FitWithRooFit * fitter,
0195     const TString & fitType, const double & xMin, const double & xMax)
0196   {
0197     gDirectory->mkdir(fileName);
0198     gDirectory->cd(fileName);
0199 
0200     // fitter->fit(histo, "", fitType, xMin, xMax, true);
0201     fitter->fit(histo, "", fitType, xMin, xMax);
0202     RooPlot * rooPlot = (RooPlot*)gDirectory->Get(TString(histo->GetName())+"_frame");
0203 
0204     gDirectory->GetMotherDir()->cd();
0205 
0206     return rooPlot;
0207   }
0208 
0209   TFile * file_;
0210   bool doFit_;
0211 };