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