Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:58:34

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