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