Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-17 01:47:51

0001 //================================================================================================
0002 //
0003 //
0004 //________________________________________________________________________________________________
0005 
0006 #if !defined(__CINT__) || defined(__MAKECINT__)
0007 #include <TROOT.h>       // access to gROOT, entry point to ROOT system
0008 #include <TSystem.h>     // interface to OS
0009 #include <TStyle.h>      // class to handle ROOT plotting style
0010 #include <TCanvas.h>     // class for drawing
0011 #include <TBenchmark.h>  // class to track macro running statistics
0012 #include <iostream>      // standard I/O
0013 
0014 // RooFit headers
0015 #include "RooRealVar.h"
0016 #include "RooDataSet.h"
0017 #include "RooCategory.h"
0018 #include "RooArgList.h"
0019 #include "RooDataHist.h"
0020 #include "RooFormulaVar.h"
0021 #include "RooHistPdf.h"
0022 #include "RooGenericPdf.h"
0023 #include "RooAddPdf.h"
0024 #include "RooSimultaneous.h"
0025 #include "RooGaussian.h"
0026 #include "RooNLLVar.h"
0027 #include "RooConstVar.h"
0028 #include "RooFitResult.h"
0029 #include "RooExponential.h"
0030 
0031 #include "TFile.h"
0032 #include "TH1D.h"
0033 #include "TPaveText.h"
0034 #include "RooPlot.h"
0035 
0036 #endif
0037 
0038 #define LUMINOSITY 2.88  //(in pb^-1)
0039 #define NBINSPASS 24
0040 #define NBINSFAIL 24
0041 
0042 using namespace RooFit;
0043 
0044 double ErrorInProduct(double x, double errx, double y, double erry, double corr) {
0045   double xFrErr = errx / x;
0046   double yFrErr = erry / y;
0047   return sqrt(pow(xFrErr, 2) + pow(yFrErr, 2) + 2.0 * corr * xFrErr * yFrErr) * x * y;
0048 }
0049 
0050 //=== MAIN MACRO =================================================================================================
0051 
0052 void performFit(string inputDir,
0053                 string PassInputDataFilename,
0054                 string FailInputDataFilename,
0055                 string PassSignalTemplateHistName,
0056                 string FailSignalTemplateHistName) {
0057   gBenchmark->Start("fitZCat");
0058 
0059   //--------------------------------------------------------------------------------------------------------------
0060   // Settings
0061   //==============================================================================================================
0062 
0063   const Double_t mlow = 60;
0064   const Double_t mhigh = 120;
0065   const Int_t nbins = 20;
0066 
0067   // Which fit parameters to set constant
0068   // MuTrk
0069   Bool_t cBgA1MuTrk = kFALSE;
0070   Bool_t cBgA2MuTrk = kTRUE;
0071   Bool_t cBgAlphaMuTrk = kFALSE;
0072   // MuMuNoIso
0073   Bool_t cBgA1MuMuNoIso = kFALSE;
0074   Bool_t cBgAlphaMuMuNoIso = kFALSE;
0075   // MuSta
0076   Bool_t cBgA1MuSta = kFALSE;
0077   Bool_t cBgAlphaMuSta = kFALSE;
0078 
0079   //--------------------------------------------------------------------------------------------------------------
0080   // Main analysis code
0081   //==============================================================================================================
0082   RooRealVar mass("mass", "mass", mlow, mhigh);
0083   mass.setBins(nbins);
0084 
0085   //
0086   // Prepare data for the fits
0087   //
0088   //Define categories
0089   RooCategory sample("sample", "");
0090   sample.defineType("Pass", 1);
0091   sample.defineType("Fail", 2);
0092 
0093   char datname[100];
0094   RooDataSet* dataPass = RooDataSet::read((inputDir + PassInputDataFilename).c_str(), RooArgList(mass));
0095   RooDataSet* dataFail = RooDataSet::read((inputDir + FailInputDataFilename).c_str(), RooArgList(mass));
0096 
0097   RooDataSet* dataCombined = new RooDataSet("dataCombined",
0098                                             "dataCombined",
0099                                             RooArgList(mass),
0100                                             Index(sample),
0101                                             Import("Pass", *dataPass),
0102                                             Import("Fail", *dataFail));
0103 
0104   //*********************************************************************************************
0105   //Define Free Parameters
0106   //*********************************************************************************************
0107   RooRealVar* ParEfficiency = new RooRealVar("ParEfficiency", "ParEfficiency", 0.9, 0.0, 1.0);
0108   RooRealVar* ParNumSignal = new RooRealVar("ParNumSignal", "ParNumSignal", 4000.0, 0.0, 100000.0);
0109   RooRealVar* ParNumBkgPass = new RooRealVar("ParNumBkgPass", "ParNumBkgPass", 1000.0, 0.0, 10000000.0);
0110   RooRealVar* ParNumBkgFail = new RooRealVar("ParNumBkgFail", "ParNumBkgFail", 1000.0, 0.0, 10000000.0);
0111 
0112   RooRealVar* ParPassBackgroundExpCoefficient =
0113       new RooRealVar("ParPassBackgroundExpCoefficient",
0114                      "ParPassBackgroundExpCoefficient",
0115                      -0.2,
0116                      -1.0,
0117                      0.0);  //ParPassBackgroundExpCoefficient.setConstant(kTRUE);
0118   RooRealVar* ParFailBackgroundExpCoefficient =
0119       new RooRealVar("ParFailBackgroundExpCoefficient",
0120                      "ParFailBackgroundExpCoefficient",
0121                      -0.2,
0122                      -1.0,
0123                      0.0);  //ParFailBackgroundExpCoefficient.setConstant(kTRUE);
0124   RooRealVar* ParPassSignalMassShift = new RooRealVar(
0125       "ParPassSignalMassShift", "ParPassSignalMassShift", -10.0, 10.0);  //ParPassSignalMassShift.setConstant(kTRUE);
0126   RooRealVar* ParFailSignalMassShift = new RooRealVar(
0127       "ParFailSignalMassShift", "ParFailSignalMassShift", -10.0, 10.0);  //ParFailSignalMassShift.setConstant(kTRUE);
0128   RooRealVar* ParPassSignalResolution = new RooRealVar(
0129       "ParPassSignalResolution", "ParPassSignalResolution", 0.0, 10.0);  //ParPassSignalResolution.setConstant(kTRUE);
0130   RooRealVar* ParFailSignalResolution = new RooRealVar(
0131       "ParFailSignalResolution", "ParFailSignalResolution", 0.0, 10.0);  //ParFailSignalResolution.setConstant(kTRUE);
0132 
0133   //*********************************************************************************************
0134   //
0135   //Load Signal PDFs
0136   //
0137   //*********************************************************************************************
0138   //Load histogram templates
0139   TFile* Zeelineshape_file = new TFile("MitFitter/templates/60To120Range/ZMassLineshape.root", "READ");
0140   TH1* histTemplatePass = (TH1D*)Zeelineshape_file->Get(PassSignalTemplateHistName.c_str());
0141   TH1* histTemplateFail = (TH1D*)Zeelineshape_file->Get(FailSignalTemplateHistName.c_str());
0142 
0143   //Introduce mass shift coordinate transformation
0144   RooFormulaVar PassShiftedMass("PassShiftedMass", "@0-@1", RooArgList(mass, *ParPassSignalMassShift));
0145   RooFormulaVar FailShiftedMass("FailShiftedMass", "@0-@1", RooArgList(mass, *ParFailSignalMassShift));
0146 
0147   RooDataHist* dataHistPass = new RooDataHist("dataHistPass", "dataHistPass", RooArgSet(mass), histTemplatePass);
0148   RooDataHist* dataHistFail = new RooDataHist("dataHistFail", "dataHistFail", RooArgSet(mass), histTemplateFail);
0149   RooHistPdf* signalShapePassPdf = new RooHistPdf("signalShapePassPdf", "signalShapePassPdf", mass, *dataHistPass, 1);
0150   RooHistPdf* signalShapeFailPdf = new RooHistPdf("signalShapeFailPdf", "signalShapeFailPdf", mass, *dataHistFail, 1);
0151 
0152   RooFormulaVar* NumSignalPass =
0153       new RooFormulaVar("NumSignalPass", "ParEfficiency*ParNumSignal", RooArgList(*ParEfficiency, *ParNumSignal));
0154   RooFormulaVar* NumSignalFail =
0155       new RooFormulaVar("NumSignalFail", "(1.0-ParEfficiency)*ParNumSignal", RooArgList(*ParEfficiency, *ParNumSignal));
0156 
0157   //*********************************************************************************************
0158   //
0159   // Create Background PDFs
0160   //
0161   //*********************************************************************************************
0162   RooExponential* bkgPassPdf = new RooExponential("bkgPassPdf", "bkgPassPdf", mass, *ParPassBackgroundExpCoefficient);
0163   RooExponential* bkgFailPdf = new RooExponential("bkgFailPdf", "bkgFailPdf", mass, *ParFailBackgroundExpCoefficient);
0164 
0165   //*********************************************************************************************
0166   //
0167   // Create Total PDFs
0168   //
0169   //*********************************************************************************************
0170   RooAddPdf pdfPass(
0171       "pdfPass", "pdfPass", RooArgList(*signalShapePassPdf, *bkgPassPdf), RooArgList(*NumSignalPass, *ParNumBkgPass));
0172   RooAddPdf pdfFail(
0173       "pdfFail", "pdfFail", RooArgList(*signalShapeFailPdf, *bkgFailPdf), RooArgList(*NumSignalFail, *ParNumBkgFail));
0174 
0175   // PDF for simultaneous fit
0176   RooSimultaneous pdfTotal("pdfTotal", "pdfTotal", sample);
0177   pdfTotal.addPdf(pdfPass, "pdfPass");
0178   pdfTotal.addPdf(pdfFail, "pdfFail");
0179   pdfTotal.Print();
0180 
0181   //
0182   // Define likelihood, add constraints, and run the fit
0183   //
0184   RooFitResult* fitResult =
0185       pdfTotal.fitTo(*dataCombined, RooFit::Save(true), RooFit::Extended(true), RooFit::PrintLevel(-1));
0186   fitResult->Print("v");
0187 
0188   double nSignalPass = NumSignalPass->getVal();
0189   double nSignalFail = NumSignalFail->getVal();
0190 
0191   printf("\nFit results:\n");
0192   printf("    Efficiency = %.4f +- %.4f\n", ParEfficiency->getVal(), ParEfficiency->getPropagatedError(*fitResult));
0193   cout << "Signal Pass: " << nSignalPass << endl;
0194   cout << "Signal Fail: " << nSignalFail << endl;
0195 
0196   //--------------------------------------------------------------------------------------------------------------
0197   // Make plots
0198   //==============================================================================================================
0199   RooAbsData::ErrorType errorType = RooAbsData::Poisson;
0200 
0201   TString cname = TString("fit_Pass");
0202   TCanvas* c = new TCanvas(cname, cname, 800, 600);
0203   RooPlot* frame1 = mass.frame();
0204   frame1->SetMinimum(0);
0205   dataPass->plotOn(frame1, RooFit::DataError(errorType));
0206   pdfPass.plotOn(frame1, RooFit::ProjWData(*dataPass), RooFit::Components(*bkgPassPdf), RooFit::LineColor(kRed));
0207   pdfPass.plotOn(frame1, RooFit::ProjWData(*dataPass));
0208   frame1->Draw("e0");
0209 
0210   TPaveText* plotlabel = new TPaveText(0.23, 0.87, 0.43, 0.92, "NDC");
0211   plotlabel->SetTextColor(kBlack);
0212   plotlabel->SetFillColor(kWhite);
0213   plotlabel->SetBorderSize(0);
0214   plotlabel->SetTextAlign(12);
0215   plotlabel->SetTextSize(0.03);
0216   plotlabel->AddText("CMS Preliminary 2010");
0217   TPaveText* plotlabel2 = new TPaveText(0.23, 0.82, 0.43, 0.87, "NDC");
0218   plotlabel2->SetTextColor(kBlack);
0219   plotlabel2->SetFillColor(kWhite);
0220   plotlabel2->SetBorderSize(0);
0221   plotlabel2->SetTextAlign(12);
0222   plotlabel2->SetTextSize(0.03);
0223   plotlabel2->AddText("#sqrt{s} = 7 TeV");
0224   TPaveText* plotlabel3 = new TPaveText(0.23, 0.75, 0.43, 0.80, "NDC");
0225   plotlabel3->SetTextColor(kBlack);
0226   plotlabel3->SetFillColor(kWhite);
0227   plotlabel3->SetBorderSize(0);
0228   plotlabel3->SetTextAlign(12);
0229   plotlabel3->SetTextSize(0.03);
0230   char temp[100];
0231   sprintf(temp, "%.4f", LUMINOSITY);
0232   plotlabel3->AddText((string("#int#font[12]{L}dt = ") + temp + string(" pb^{ -1}")).c_str());
0233   TPaveText* plotlabel4 = new TPaveText(0.6, 0.82, 0.8, 0.87, "NDC");
0234   plotlabel4->SetTextColor(kBlack);
0235   plotlabel4->SetFillColor(kWhite);
0236   plotlabel4->SetBorderSize(0);
0237   plotlabel4->SetTextAlign(12);
0238   plotlabel4->SetTextSize(0.03);
0239   double nsig = ParNumSignal->getVal();
0240   double nErr = ParNumSignal->getError();
0241   double e = ParEfficiency->getVal();
0242   double eErr = ParEfficiency->getError();
0243   double corr = fitResult->correlation(*ParEfficiency, *ParNumSignal);
0244   double err = ErrorInProduct(nsig, nErr, e, eErr, corr);
0245   sprintf(temp, "Signal = %.2f #pm %.2f", NumSignalPass->getVal(), err);
0246   plotlabel4->AddText(temp);
0247   TPaveText* plotlabel5 = new TPaveText(0.6, 0.77, 0.8, 0.82, "NDC");
0248   plotlabel5->SetTextColor(kBlack);
0249   plotlabel5->SetFillColor(kWhite);
0250   plotlabel5->SetBorderSize(0);
0251   plotlabel5->SetTextAlign(12);
0252   plotlabel5->SetTextSize(0.03);
0253   sprintf(temp, "Bkg = %.2f #pm %.2f", ParNumBkgPass->getVal(), ParNumBkgPass->getError());
0254   plotlabel5->AddText(temp);
0255   TPaveText* plotlabel6 = new TPaveText(0.6, 0.87, 0.8, 0.92, "NDC");
0256   plotlabel6->SetTextColor(kBlack);
0257   plotlabel6->SetFillColor(kWhite);
0258   plotlabel6->SetBorderSize(0);
0259   plotlabel6->SetTextAlign(12);
0260   plotlabel6->SetTextSize(0.03);
0261   plotlabel6->AddText("Passing probes");
0262   TPaveText* plotlabel7 = new TPaveText(0.6, 0.72, 0.8, 0.77, "NDC");
0263   plotlabel7->SetTextColor(kBlack);
0264   plotlabel7->SetFillColor(kWhite);
0265   plotlabel7->SetBorderSize(0);
0266   plotlabel7->SetTextAlign(12);
0267   plotlabel7->SetTextSize(0.03);
0268   sprintf(temp, "Eff = %.3f #pm %.3f", ParEfficiency->getVal(), ParEfficiency->getErrorHi());
0269   plotlabel7->AddText(temp);
0270 
0271   plotlabel4->Draw();
0272   plotlabel5->Draw();
0273   plotlabel6->Draw();
0274   plotlabel7->Draw();
0275 
0276   //   c->SaveAs( cname + TString(".eps"));
0277   c->SaveAs(cname + TString(".gif"));
0278   delete c;
0279 
0280   cname = TString("fit_Fail");
0281   TCanvas* c2 = new TCanvas(cname, cname, 500, 500);
0282   RooPlot* frame2 = mass.frame();
0283   frame2->SetMinimum(0);
0284   dataFail->plotOn(frame2, RooFit::DataError(errorType));
0285   pdfFail.plotOn(frame2, RooFit::ProjWData(*dataFail), RooFit::Components(*bkgFailPdf), RooFit::LineColor(kRed));
0286   pdfFail.plotOn(frame2, RooFit::ProjWData(*dataFail));
0287   frame2->Draw("e0");
0288 
0289   plotlabel = new TPaveText(0.23, 0.87, 0.43, 0.92, "NDC");
0290   plotlabel->SetTextColor(kBlack);
0291   plotlabel->SetFillColor(kWhite);
0292   plotlabel->SetBorderSize(0);
0293   plotlabel->SetTextAlign(12);
0294   plotlabel->SetTextSize(0.03);
0295   plotlabel->AddText("CMS Preliminary 2010");
0296   plotlabel2 = new TPaveText(0.23, 0.82, 0.43, 0.87, "NDC");
0297   plotlabel2->SetTextColor(kBlack);
0298   plotlabel2->SetFillColor(kWhite);
0299   plotlabel2->SetBorderSize(0);
0300   plotlabel2->SetTextAlign(12);
0301   plotlabel2->SetTextSize(0.03);
0302   plotlabel2->AddText("#sqrt{s} = 7 TeV");
0303   plotlabel3 = new TPaveText(0.23, 0.75, 0.43, 0.80, "NDC");
0304   plotlabel3->SetTextColor(kBlack);
0305   plotlabel3->SetFillColor(kWhite);
0306   plotlabel3->SetBorderSize(0);
0307   plotlabel3->SetTextAlign(12);
0308   plotlabel3->SetTextSize(0.03);
0309   sprintf(temp, "%.4f", LUMINOSITY);
0310   plotlabel3->AddText((string("#int#font[12]{L}dt = ") + temp + string(" pb^{ -1}")).c_str());
0311   plotlabel4 = new TPaveText(0.6, 0.82, 0.8, 0.87, "NDC");
0312   plotlabel4->SetTextColor(kBlack);
0313   plotlabel4->SetFillColor(kWhite);
0314   plotlabel4->SetBorderSize(0);
0315   plotlabel4->SetTextAlign(12);
0316   plotlabel4->SetTextSize(0.03);
0317   err = ErrorInProduct(nsig, nErr, 1.0 - e, eErr, corr);
0318   sprintf(temp, "Signal = %.2f #pm %.2f", NumSignalFail->getVal(), err);
0319   plotlabel4->AddText(temp);
0320   plotlabel5 = new TPaveText(0.6, 0.77, 0.8, 0.82, "NDC");
0321   plotlabel5->SetTextColor(kBlack);
0322   plotlabel5->SetFillColor(kWhite);
0323   plotlabel5->SetBorderSize(0);
0324   plotlabel5->SetTextAlign(12);
0325   plotlabel5->SetTextSize(0.03);
0326   sprintf(temp, "Bkg = %.2f #pm %.2f", ParNumBkgFail->getVal(), ParNumBkgFail->getError());
0327   plotlabel5->AddText(temp);
0328   plotlabel6 = new TPaveText(0.6, 0.87, 0.8, 0.92, "NDC");
0329   plotlabel6->SetTextColor(kBlack);
0330   plotlabel6->SetFillColor(kWhite);
0331   plotlabel6->SetBorderSize(0);
0332   plotlabel6->SetTextAlign(12);
0333   plotlabel6->SetTextSize(0.03);
0334   plotlabel6->AddText("Failing probes");
0335   plotlabel7 = new TPaveText(0.6, 0.72, 0.8, 0.77, "NDC");
0336   plotlabel7->SetTextColor(kBlack);
0337   plotlabel7->SetFillColor(kWhite);
0338   plotlabel7->SetBorderSize(0);
0339   plotlabel7->SetTextAlign(12);
0340   plotlabel7->SetTextSize(0.03);
0341   sprintf(
0342       temp, "Eff = %.3f #pm %.3f", ParEfficiency->getVal(), ParEfficiency->getErrorHi(), ParEfficiency->getErrorLo());
0343   plotlabel7->AddText(temp);
0344 
0345   //   plotlabel->Draw();
0346   //   plotlabel2->Draw();
0347   //   plotlabel3->Draw();
0348   plotlabel4->Draw();
0349   plotlabel5->Draw();
0350   plotlabel6->Draw();
0351   plotlabel7->Draw();
0352 
0353   c2->SaveAs(cname + TString(".eps"));
0354   c2->SaveAs(cname + TString(".gif"));
0355   c2->SaveAs(cname + TString(".root"));
0356   delete c2;
0357 
0358   gBenchmark->Show("fitZCat");
0359 }
0360 
0361 void fitZCat() {
0362   performFit("EfficiencyFitter/input/Data_36_09122010_TP_TagWP80/",
0363              "Mass_TagPlusRecoPassVBTF95IdIso_EB_Pt20ToInf_Data_36",
0364              "Mass_TagPlusRecoFailVBTF95IdIso_EB_Pt20ToInf_Data_36",
0365              "Mass_TagPlusRecoPassVBTF95IdIso_EB_Pt20ToInf",
0366              "Mass_TagPlusRecoFailVBTF95IdIso_EB_Pt20ToInf");
0367 }