Macros/RooFit/CompareBiasUpsilonValidation.cc

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 };