Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:31:36

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