Macros/RooFit/FitXslices.cc

0001 #ifndef FitXslices_cc
0002 #define FitXslices_cc
0003
0004 #include <iostream>
0005 #include <sstream>
0006 #include "TColor.h"
0007 #include "TH2F.h"
0008 #include "TH3F.h"
0009 #include "TDirectory.h"
0010 #include "TROOT.h"
0011 #include "TStyle.h"
0012 #include "FitWithRooFit.cc"
0013
0014 /**
0015  * This class performs the following actions: <br>
0016  * - take a TH2* as input and fit slices of it in the x coordinate (one slice per bin) <br>
0017  * - store the result of each slice fit and draw them on a canvas <br>
0018  * - draw plots of each parameter result with the corresponding error <br>
0019  * It uses RooFit for the fitting.
0020  */
0021
0022 class FitXslices {
0023 public:
0024   FitXslices() {
0025     /// Initialize suitable parameter values for Z fitting
0026     /// N.B.: mean and sigma of gaussian, as well as mean and range of the peak are defined in the CompareBiasZValidation class
0027
0028     /// Fraction of signal
0029     fitter_.initFsig(0.9, 0., 1.);
0030
0031     /// CB
0032     fitter_.initMean2(0., -20., 20.);
0033     fitter_.mean2()->setConstant(kTRUE);
0034     fitter_.initSigma(1.5, 0.5, 3.);
0035     fitter_.initSigma2(1.2, 0., 5.);
0036     fitter_.initAlpha(1.5, 0.05, 10.);
0037     fitter_.initN(1, 0.01, 100.);
0038
0039     // BW Fix the gamma for the Z
0040     fitter_.initGamma(2.4952, 0., 10.);
0041     fitter_.gamma()->setConstant(kTRUE);  // This can be overriden as well.
0042
0043     fitter_.initGaussFrac(0.5, 0., 1.);
0044
0045     /// Exponential with pol2 argument
0046     fitter_.initExpCoeffA0(-1., -10., 10.);
0047     fitter_.initExpCoeffA1(0., -10., 10.);
0048     fitter_.initExpCoeffA2(0., -2., 2.);
0049
0050     /// Polynomial
0051     fitter_.initA0(0., -10., 10.);
0052     fitter_.initA1(0., -10., 10.);
0053     fitter_.initA2(0., -10., 10.);
0054     fitter_.initA3(0., -10., 10.);
0055     fitter_.initA4(0., -10., 10.);
0056     fitter_.initA5(0., -10., 10.);
0057     fitter_.initA6(0., -10., 10.);
0058
0059     fitter_.useChi2_ = false;
0060   }
0061
0062   FitWithRooFit* fitter() { return (&fitter_); }
0063
0064   //  void fitSlices( std::map<unsigned int, TH1*> & slices, const double & xMin, const double & xMax, const TString & signalType, const TString & backgroundType, const bool twoD ){    }
0065
0066   void operator()(TH2* histo,
0067                   const double& xMin,
0068                   const double& xMax,
0069                   const TString& signalType,
0070                   const TString& backgroundType,
0071                   unsigned int rebinY) {
0072     // Create and move in a subdir
0073     gDirectory->mkdir("allHistos");
0074     gDirectory->cd("allHistos");
0075
0076     gStyle->SetPalette(1);
0077     // Loop on all X bins, project on Y and fit the resulting TH1
0078     TString name = histo->GetName();
0079     unsigned int binsX = histo->GetNbinsX();
0080
0081     // The canvas for the results of the fit (the mean values for the gaussians +- errors)
0082     TCanvas* meanCanvas = new TCanvas("meanCanvas", "meanCanvas", 1000, 800);
0083     TH1D* meanHisto =
0084         new TH1D("meanHisto", "meanHisto", binsX, histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0085
0086     TCanvas* sigmaCanvas = new TCanvas("sigmaCanvas", "sigmaCanvas", 1000, 800);
0087     TH1D* sigmaHisto =
0088         new TH1D("sigmaHisto", "sigmaHisto", binsX, histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0089
0090     TCanvas* backgroundCanvas = new TCanvas("backgroundCanvas", "backgroundCanvas", 1000, 800);
0091     TH1D* backgroundHisto = new TH1D(
0092         "backgroundHisto", "backgroundHisto", binsX, histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0093
0094     TCanvas* backgroundCanvas2 = new TCanvas("backgroundCanvas2", "backgroundCanvas2", 1000, 800);
0095     TH1D* backgroundHisto2 =
0096         new TH1D("backgroundHisto2", "exp a1", binsX, histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0097
0098     TCanvas* backgroundCanvas3 = new TCanvas("backgroundCanvas3", "backgroundCanvas3", 1000, 800);
0099     TH1D* backgroundHisto3 =
0100         new TH1D("backgroundHisto3", "exp a2", binsX, histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0101
0102     TCanvas* signalFractionCanvas = new TCanvas("signalFractionCanvas", "signalFractionCanvas", 1000, 800);
0103     TH1D* signalFractionHisto = new TH1D(
0104         "signalFractionHisto", "signalFractionHisto", binsX, histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0105
0106     TCanvas* probChi2Canvas = new TCanvas("probChi2Canvas", "probChi2Canvas", 1000, 800);
0107     TH1D* probChi2Histo =
0108         new TH1D("probChi2Histo", "probChi2Histo", binsX, histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
0109
0110     // Store all the non-empty slices
0111     std::map<unsigned int, TH1*> slices;
0112     for (unsigned int x = 1; x <= binsX; ++x) {
0113       std::stringstream ss;
0114       ss << x;
0115       TH1* sliceHisto = histo->ProjectionY((std::string{name.Data()} + ss.str()).c_str(), x, x);
0116       if (sliceHisto->GetEntries() != 0) {
0117         // std::cout << "filling for x = " << x << endl;
0118         slices.insert(std::make_pair(x, sliceHisto));
0119         sliceHisto->Rebin(rebinY);
0120       }
0121     }
0122
0123     // Create the canvas for all the fits
0124     TCanvas* fitsCanvas = new TCanvas("fitsCanvas", "fits canvas", 1000, 800);
0125     // cout << "slices.size = " << slices.size() << endl;
0126     unsigned int x = sqrt(slices.size());
0127     unsigned int y = x;
0128     if (x * y < slices.size()) {
0129       x += 1;
0130       y += 1;
0131     }
0132     fitsCanvas->Divide(x, y);
0133
0134     // Loop on the saved slices and fit
0135     std::map<unsigned int, TH1*>::iterator it = slices.begin();
0136     unsigned int i = 1;
0137     for (; it != slices.end(); ++it, ++i) {
0138       fitsCanvas->cd(i);
0139       fitter_.fit(it->second, signalType, backgroundType, xMin, xMax);
0140       //      fitsCanvas->GetPad(i)->SetLogy();
0141       // FIXME: prob(chi2) needs to be computed properly inside FitWithRooFit.cc
0142       // probChi2Histo->SetBinContent(it->first, mean->getVal());
0143
0144       RooRealVar* mean = fitter_.mean();
0145
0146       meanHisto->SetBinContent(it->first, mean->getVal());
0147       meanHisto->SetBinError(it->first, mean->getError());
0148
0149       RooRealVar* sigma = fitter_.sigma();
0150       sigmaHisto->SetBinContent(it->first, sigma->getVal());
0151       sigmaHisto->SetBinError(it->first, sigma->getError());
0152
0153       std::cout << "backgroundType = " << backgroundType << std::endl;
0154       if (backgroundType == "exponential") {
0155         RooRealVar* expCoeff = fitter_.expCoeffa1();
0156         backgroundHisto->SetBinContent(it->first, expCoeff->getVal());
0157         backgroundHisto->SetBinError(it->first, expCoeff->getError());
0158       } else if (backgroundType == "exponentialpol") {
0159         RooRealVar* expCoeffa0 = fitter_.expCoeffa0();
0160         backgroundHisto->SetBinContent(it->first, expCoeffa0->getVal());
0161         backgroundHisto->SetBinError(it->first, expCoeffa0->getError());
0162
0163         RooRealVar* expCoeffa1 = fitter_.expCoeffa1();
0164         backgroundHisto2->SetBinContent(it->first, expCoeffa1->getVal());
0165         backgroundHisto2->SetBinError(it->first, expCoeffa1->getError());
0166
0167         RooRealVar* expCoeffa2 = fitter_.expCoeffa2();
0168         backgroundHisto3->SetBinContent(it->first, expCoeffa2->getVal());
0169         backgroundHisto3->SetBinError(it->first, expCoeffa2->getError());
0170       }
0171
0172       else if (backgroundType == "linear") {
0173         RooRealVar* linearTerm = fitter_.a1();
0174         backgroundHisto->SetBinContent(it->first, linearTerm->getVal());
0175         backgroundHisto->SetBinError(it->first, linearTerm->getError());
0176
0177         RooRealVar* constant = fitter_.a0();
0178         backgroundHisto2->SetBinContent(it->first, constant->getVal());
0179         backgroundHisto2->SetBinError(it->first, constant->getError());
0180       }
0181       RooRealVar* fsig = fitter_.fsig();
0182       signalFractionHisto->SetBinContent(it->first, fsig->getVal());
0183       signalFractionHisto->SetBinError(it->first, fsig->getError());
0184     }
0185     // Go back to the main dir before saving the canvases
0186     gDirectory->GetMotherDir()->cd();
0187     meanCanvas->cd();
0188     meanHisto->Draw();
0189     sigmaCanvas->cd();
0190     sigmaHisto->Draw();
0191     backgroundCanvas->cd();
0192     backgroundHisto->Draw();
0193     if (backgroundType == "linear") {
0194       backgroundCanvas2->cd();
0195       backgroundHisto2->Draw();
0196     }
0197     if (backgroundType == "exponentialpol") {
0198       backgroundCanvas2->cd();
0199       backgroundHisto2->Draw();
0200       backgroundCanvas3->cd();
0201       backgroundHisto3->Draw();
0202     }
0203     signalFractionCanvas->cd();
0204     signalFractionHisto->Draw();
0205     probChi2Canvas->cd();
0206     probChi2Histo->Draw();
0207
0208     fitsCanvas->Write();
0209     meanCanvas->Write();
0210     sigmaCanvas->Write();
0211     backgroundCanvas->Write();
0212     signalFractionCanvas->Write();
0213     if (backgroundType == "linear") {
0214       backgroundCanvas2->Write();
0215     }
0216     probChi2Canvas->Write();
0217
0218     fitsCanvas->Close();
0219     probChi2Canvas->Close();
0220     signalFractionCanvas->Close();
0221     backgroundCanvas3->Close();
0222     backgroundCanvas2->Close();
0223     backgroundCanvas->Close();
0224     sigmaCanvas->Close();
0225     meanCanvas->Close();
0226   }
0227
0228   void operator()(TH3* histo,
0229                   const double& xMin,
0230                   const double& xMax,
0231                   const TString& signalType,
0232                   const TString& backgroundType,
0233                   unsigned int rebinZ) {
0234     // Create and move in a subdir
0235     gDirectory->mkdir("allHistos");
0236     gDirectory->cd("allHistos");
0237
0238     // Loop on all X bins, project on Y and fit the resulting TH2
0239     TString name = histo->GetName();
0240     unsigned int binsX = histo->GetNbinsX();
0241     unsigned int binsY = histo->GetNbinsY();
0242
0243     // std::cout<< "number of bins in x --> "<<binsX<<std::endl;
0244     // std::cout<< "number of bins in y --> "<<binsY<<std::endl;
0245
0246     // The canvas for the results of the fit (the mean values for the gaussians +- errors)
0247     TCanvas* meanCanvas = new TCanvas("meanCanvas", "meanCanvas", 1000, 800);
0248     TH2D* meanHisto = new TH2D("meanHisto",
0249                                "meanHisto",
0250                                binsX,
0251                                histo->GetXaxis()->GetXmin(),
0252                                histo->GetXaxis()->GetXmax(),
0253                                binsY,
0254                                histo->GetYaxis()->GetXmin(),
0255                                histo->GetYaxis()->GetXmax());
0256
0257     TCanvas* errorMeanCanvas = new TCanvas("errorMeanCanvas", "errorMeanCanvas", 1000, 800);
0258     TH2D* errorMeanHisto = new TH2D("errorMeanHisto",
0259                                     "errorMeanHisto",
0260                                     binsX,
0261                                     histo->GetXaxis()->GetXmin(),
0262                                     histo->GetXaxis()->GetXmax(),
0263                                     binsY,
0264                                     histo->GetYaxis()->GetXmin(),
0265                                     histo->GetYaxis()->GetXmax());
0266
0267     TCanvas* sigmaCanvas = new TCanvas("sigmaCanvas", "sigmaCanvas", 1000, 800);
0268     TH2D* sigmaHisto = new TH2D("sigmaHisto",
0269                                 "sigmaHisto",
0270                                 binsX,
0271                                 histo->GetXaxis()->GetXmin(),
0272                                 histo->GetXaxis()->GetXmax(),
0273                                 binsY,
0274                                 histo->GetYaxis()->GetXmin(),
0275                                 histo->GetYaxis()->GetXmax());
0276
0277     TCanvas* backgroundCanvas = new TCanvas("backgroundCanvas", "backgroundCanvas", 1000, 800);
0278     TH2D* backgroundHisto = new TH2D("backgroundHisto",
0279                                      "backgroundHisto",
0280                                      binsX,
0281                                      histo->GetXaxis()->GetXmin(),
0282                                      histo->GetXaxis()->GetXmax(),
0283                                      binsY,
0284                                      histo->GetYaxis()->GetXmin(),
0285                                      histo->GetYaxis()->GetXmax());
0286     TCanvas* backgroundCanvas2 = new TCanvas("backgroundCanvas2", "backgroundCanvas2", 1000, 800);
0287     TH2D* backgroundHisto2 = new TH2D("backgroundHisto2",
0288                                       "a1",
0289                                       binsX,
0290                                       histo->GetXaxis()->GetXmin(),
0291                                       histo->GetXaxis()->GetXmax(),
0292                                       binsY,
0293                                       histo->GetYaxis()->GetXmin(),
0294                                       histo->GetYaxis()->GetXmax());
0295     TCanvas* backgroundCanvas3 = new TCanvas("backgroundCanvas3", "backgroundCanvas3", 1000, 800);
0296     TH2D* backgroundHisto3 = new TH2D("backgroundHisto3",
0297                                       "a2",
0298                                       binsX,
0299                                       histo->GetXaxis()->GetXmin(),
0300                                       histo->GetXaxis()->GetXmax(),
0301                                       binsY,
0302                                       histo->GetYaxis()->GetXmin(),
0303                                       histo->GetYaxis()->GetXmax());
0304
0305     TCanvas* signalFractionCanvas = new TCanvas("signalFractionCanvas", "signalFractionCanvas", 1000, 800);
0306     TH2D* signalFractionHisto = new TH2D("signalFractionHisto",
0307                                          "signalFractionHisto",
0308                                          binsX,
0309                                          histo->GetXaxis()->GetXmin(),
0310                                          histo->GetXaxis()->GetXmax(),
0311                                          binsY,
0312                                          histo->GetYaxis()->GetXmin(),
0313                                          histo->GetYaxis()->GetXmax());
0314
0315     // Store all the non-empty slices
0316     std::map<unsigned int, TH1*> slices;
0317     for (unsigned int x = 1; x <= binsX; ++x) {
0318       for (unsigned int y = 1; y <= binsY; ++y) {
0319         std::stringstream ss;
0320         ss << x << "_" << y;
0321         TH1* sliceHisto = histo->ProjectionZ((std::string{name.Data()} + ss.str()).c_str(), x, x, y, y);
0322         if (sliceHisto->GetEntries() != 0) {
0323           sliceHisto->Rebin(rebinZ);
0324           // std::cout << "filling for x = " << x << endl;
0325           slices.insert(std::make_pair(x + (binsX + 1) * y, sliceHisto));
0326         }
0327       }
0328     }
0329
0330     // Create the canvas for all the fits
0331     TCanvas* fitsCanvas = new TCanvas("fitsCanvas", "canvas of all fits", 1000, 800);
0332     // cout << "slices.size = " << slices.size() << endl;
0333     unsigned int x = sqrt(slices.size());
0334     unsigned int y = x;
0335     if (x * y < slices.size()) {
0336       x += 1;
0337       y += 1;
0338     }
0339     fitsCanvas->Divide(x, y);
0340
0341     // Loop on the saved slices and fit
0342     std::map<unsigned int, TH1*>::iterator it = slices.begin();
0343     unsigned int i = 1;
0344     for (; it != slices.end(); ++it, ++i) {
0345       fitsCanvas->cd(i);
0346
0347       fitter_.fit(it->second, signalType, backgroundType, xMin, xMax);
0348
0349       RooRealVar* mean = fitter_.mean();
0350       meanHisto->SetBinContent(it->first % (binsX + 1), int(it->first / (binsX + 1)), mean->getVal());
0351       errorMeanHisto->SetBinContent(it->first % (binsX + 1), int(it->first / (binsX + 1)), mean->getError());
0352       //      meanHisto->SetBinError(it->first%binsX, int(it->first/binsX), mean->getError());
0353       //std::cout<<"int i -->"<<i<<std::endl;
0354       //std::cout<< " it->first%(binsX+1) --> "<<it->first%(binsX+1)<<std::endl;
0355       //std::cout<< " it->first/(binsX+1) --> "<<int(it->first/(binsX+1))<<std::endl;
0356
0357       RooRealVar* sigma = fitter_.sigma();
0358       sigmaHisto->SetBinContent(it->first % binsX, int(it->first / binsX), sigma->getVal());
0359       sigmaHisto->SetBinError(it->first % binsX, int(it->first / binsX), sigma->getError());
0360
0361       std::cout << "backgroundType = " << backgroundType << std::endl;
0362       if (backgroundType == "exponential") {
0363         RooRealVar* expCoeff = fitter_.expCoeffa1();
0364         backgroundHisto->SetBinContent(it->first % binsX, int(it->first / binsX), expCoeff->getVal());
0365         backgroundHisto->SetBinError(it->first % binsX, int(it->first / binsX), expCoeff->getError());
0366       } else if (backgroundType == "exponentialpol") {
0367         RooRealVar* expCoeffa0 = fitter_.expCoeffa0();
0368         backgroundHisto->SetBinContent(it->first % binsX, int(it->first / binsX), expCoeffa0->getVal());
0369         backgroundHisto->SetBinError(it->first % binsX, int(it->first / binsX), expCoeffa0->getError());
0370
0371         RooRealVar* expCoeffa1 = fitter_.expCoeffa1();
0372         backgroundHisto2->SetBinContent(it->first % binsX, int(it->first / binsX), expCoeffa1->getVal());
0373         backgroundHisto2->SetBinError(it->first % binsX, int(it->first / binsX), expCoeffa1->getError());
0374
0375         RooRealVar* expCoeffa2 = fitter_.expCoeffa2();
0376         backgroundHisto3->SetBinContent(it->first % binsX, int(it->first / binsX), expCoeffa2->getVal());
0377         backgroundHisto3->SetBinError(it->first % binsX, int(it->first / binsX), expCoeffa2->getError());
0378       } else if (backgroundType == "linear") {
0379         RooRealVar* linearTerm = fitter_.a1();
0380         backgroundHisto->SetBinContent(it->first % binsX, int(it->first / binsX), linearTerm->getVal());
0381         backgroundHisto->SetBinError(it->first % binsX, int(it->first / binsX), linearTerm->getError());
0382
0383         RooRealVar* constant = fitter_.a0();
0384         backgroundHisto2->SetBinContent(it->first % binsX, int(it->first / binsX), constant->getVal());
0385         backgroundHisto2->SetBinError(it->first % binsX, int(it->first / binsX), constant->getError());
0386       }
0387
0388       RooRealVar* fsig = fitter_.fsig();
0389       signalFractionHisto->SetBinContent(it->first % binsX, int(it->first / binsX), fsig->getVal());
0390       signalFractionHisto->SetBinError(it->first % binsX, int(it->first / binsX), fsig->getError());
0391     }
0392     // Go back to the main dir before saving the canvases
0393     gDirectory->GetMotherDir()->cd();
0394
0395     meanCanvas->cd();
0396     meanHisto->GetXaxis()->SetRangeUser(-3.14, 3.14);
0397     meanHisto->GetYaxis()->SetRangeUser(-2.5, 2.5);
0398     meanHisto->GetXaxis()->SetTitle("#phi (rad)");
0399     meanHisto->GetYaxis()->SetTitle("#eta");
0400     meanHisto->Draw("COLZ");
0401
0402     sigmaCanvas->cd();
0403     sigmaHisto->GetXaxis()->SetRangeUser(-3.14, 3.14);
0404     sigmaHisto->GetYaxis()->SetRangeUser(-2.5, 2.5);
0405     sigmaHisto->GetXaxis()->SetTitle("#phi (rad)");
0406     sigmaHisto->GetYaxis()->SetTitle("#eta");
0407     sigmaHisto->Draw("COLZ");
0408
0409     backgroundCanvas->cd();
0410     backgroundHisto->Draw("COLZ");
0411     if (backgroundType == "linear") {
0412       backgroundCanvas2->cd();
0413       backgroundHisto2->Draw("COLZ");
0414     }
0415     signalFractionCanvas->cd();
0416     signalFractionHisto->Draw("COLZ");
0417
0418     fitsCanvas->Write();
0419     meanCanvas->Write();
0420     sigmaCanvas->Write();
0421     backgroundCanvas->Write();
0422     signalFractionCanvas->Write();
0423     if (backgroundType == "linear") {
0424       backgroundCanvas2->Write();
0425     }
0426
0427     fitsCanvas->Close();
0428     signalFractionCanvas->Close();
0429     backgroundCanvas3->Close();
0430     backgroundCanvas2->Close();
0431     backgroundCanvas->Close();
0432     sigmaCanvas->Close();
0433     errorMeanCanvas->Close();
0434     meanCanvas->Close();
0435   }
0436
0437 protected:
0438   struct Parameter {
0439     Parameter(const double& inputValue, const double& inputError) : value(inputValue), error(inputError) {}
0440
0441     double value;
0442     double error;
0443   };
0444
0445   FitWithRooFit fitter_;
0446 };
0447
0448 #endif