Macros/RooFit/FitSlices.cc

0001 #ifndef FitSlices_cc
0002 #define FitSlices_cc
0003
0004 #include "FitXslices.cc"
0005 #include "TFile.h"
0006 #include "TH1F.h"
0007 #include "TH2F.h"
0008 #include "TH3F.h"
0009 #include "TROOT.h"
0010
0011 /**
0012  * This class can be used to fit the X slices of a TH1 histogram using RooFit.
0013  * It uses the FitXslices class to do the fitting.
0014  */
0015 class FitSlices {
0016 public:
0017   double xMean;
0018   double xMin;
0019   double xMax;
0020   double sigma;
0021   double sigmaMin;
0022   double sigmaMax;
0023   TString signalType;
0024   TString backgroundType;
0025
0026   unsigned int rebinX;
0027   unsigned int rebinY;
0028   unsigned int rebinZ;
0029   double sigma2, sigma2Min, sigma2Max;
0030   bool useChi2;
0031
0032   FitXslices fitXslices;
0033
0034   FitSlices(double xMean_,
0035             double xMin_,
0036             double xMax_,
0037             double sigma_,
0038             double sigmaMin_,
0039             double sigmaMax_,
0040             TString signalType_,
0041             TString backgroundType_)
0042       : xMean(xMean_),
0043         xMin(xMin_),
0044         xMax(xMax_),
0045         sigma(sigma_),
0046         sigmaMin(sigmaMin_),
0047         sigmaMax(sigmaMax_),
0048         signalType(signalType_),
0049         backgroundType(backgroundType_),
0050         rebinX(2),
0051         rebinY(2),
0052         rebinZ(2),
0053         sigma2(0.1),
0054         sigma2Min(0.),
0055         sigma2Max(10.),
0056         useChi2(false) {
0057     // Initialize fitXslices to some defaults, nothing un-overridable
0058     fitXslices.fitter()->useChi2_ = useChi2;
0059     fitXslices.fitter()->initMean(xMean, xMin, xMax);
0060     fitXslices.fitter()->initSigma(sigma, sigmaMin, sigmaMax);
0061     fitXslices.fitter()->initSigma2(sigma2, sigma2Min, sigma2Max);
0062     fitXslices.fitter()->initAlpha(1.5, 0.05, 10.);
0063     fitXslices.fitter()->initN(1, 0.01, 100.);
0064     fitXslices.fitter()->initFGCB(0.4, 0., 1.);
0065   }
0066
0067   // virtual void fit(const TString & inputFileName = "0_MuScleFit.root", const TString & outputFileName = "BiasCheck_0.root",
0068   //           const TString & signalType = "gaussian", const TString & backgroundType = "exponential",
0069   //           const double & xMean = 3.1, const double & xMin = 3., const double & xMax = 3.2,
0070   //           const double & sigma = 0.03, const double & sigmaMin = 0., const double & sigmaMax = 0.1,
0071   //           const TString & histoBaseName = "hRecBestResVSMu", const TString & histoBaseTitle = "MassVs") = 0;
0072
0073   void fitSlice(const TString& histoName, const TString& dirName, TFile* inputFile, TDirectory* outputFile) {
0074     //r.c. patch --------------
0075     if (histoName == "hRecBestResVSMu_MassVSEtaPhiPlus" || histoName == "hRecBestResVSMu_MassVSEtaPhiMinus" ||
0076         histoName == "hRecBestResVSMu_MassVSPhiPlusPhiMinus" || histoName == "hRecBestResVSMu_MassVSEtaPlusEtaMinus") {
0077       TH3* histoPt3 = (TH3*)inputFile->FindObjectAny(histoName);
0078       outputFile->mkdir(dirName);
0079       outputFile->cd(dirName);
0080
0081       //    histoPt3 = rebin3D(histoPt3);
0082       fitXslices(histoPt3, xMin, xMax, signalType, backgroundType, rebinZ);
0083
0084       //    histoPt3->RebinX(rebinX);
0085       //    histoPt3->RebinY(rebinX);
0086       //    histoPt3->RebinY(rebinY);
0087       //    (histoPt3->DoProject2D())->RebinX(rebinX);
0088       //    (histoPt3->DoProject2D())->RebinY(rebinY);
0089     } else {
0090       TH2* histoPt2 = (TH2*)inputFile->FindObjectAny(histoName);
0091       histoPt2->RebinX(rebinX);
0092       histoPt2->RebinY(rebinY);
0093       // TH2 * histoPt = 0;
0094       // inputFile->GetObject(histoName, histoPt);
0095       outputFile->mkdir(dirName);
0096       outputFile->cd(dirName);
0097       fitXslices(histoPt2, xMin, xMax, signalType, backgroundType, rebinZ);
0098     }
0099   }
0100
0101   TH3* rebin3D(const TH3* histo3D) {
0102     unsigned int zbins = histo3D->GetNbinsZ();
0103     // std::cout<< "number of bins in z (and tempHisto) --> "<<zbins<<std::endl;
0104     std::map<unsigned int, TH2*> twoDprojection;
0105     for (unsigned int z = 1; z < zbins; ++z) {
0106       TAxis* ax_tmp = const_cast<TAxis*>(histo3D->GetZaxis());
0107       ax_tmp->SetRange(z, z);
0108       TH2* tempHisto = (TH2*)histo3D->Project3D("xy");
0109       std::stringstream ss;
0110       ss << z;
0111       tempHisto->SetName((std::string{tempHisto->GetName()} + ss.str()).c_str());
0112       tempHisto->RebinX(rebinX);
0113       tempHisto->RebinY(rebinY);
0114       twoDprojection.insert(std::make_pair(z, tempHisto));
0115     }
0116     unsigned int xbins, ybins;
0117     TH3* rebinned3D = (TH3*)new TH3F(TString(histo3D->GetName()) + "_rebinned",
0118                                      histo3D->GetTitle(),
0119                                      xbins,
0120                                      histo3D->GetXaxis()->GetXmin(),
0121                                      histo3D->GetXaxis()->GetXmax(),
0122                                      ybins,
0123                                      histo3D->GetYaxis()->GetXmin(),
0124                                      histo3D->GetYaxis()->GetXmax(),
0125                                      zbins,
0126                                      histo3D->GetZaxis()->GetXmin(),
0127                                      histo3D->GetZaxis()->GetXmax());
0128     if (twoDprojection.size() != 0) {
0129       xbins = twoDprojection[1]->GetNbinsX();
0130       ybins = twoDprojection[1]->GetNbinsY();
0131       //std::cout<< "number of bins in x --> "<<xbins<<std::endl;
0132       //std::cout<< "number of bins in y --> "<<ybins<<std::endl;
0133       for (unsigned int z = 1; z < zbins; ++z) {
0134         for (unsigned int y = 1; y < ybins; ++y) {
0135           for (unsigned int x = 1; x < xbins; ++x) {
0136             std::cout << "x/y/z= " << x << "/" << y << "/" << z << std::endl;
0137             std::cout << "number of bins in x --> " << xbins << std::endl;
0138             std::cout << "number of bins in y --> " << ybins << std::endl;
0139             rebinned3D->Fill(x, y, twoDprojection[z]->GetBinContent(x, y));
0140           }
0141         }
0142       }
0143     }
0144     return rebinned3D;
0145   }
0146 };
0147
0148 #endif