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 #include "RooFFTConvPdf.h"
0031 
0032 #include "TFile.h"
0033 #include "TH1D.h"
0034 #include "TPaveText.h"
0035 #include "RooPlot.h"
0036 
0037 #endif
0038 
0039 #define LUMINOSITY 2.88  //(in pb^-1)
0040 #define NBINSPASS 24
0041 #define NBINSFAIL 24
0042 
0043 using namespace RooFit;
0044 
0045 double ErrorInProduct(double x, double errx, double y, double erry, double corr) {
0046   double xFrErr = errx / x;
0047   double yFrErr = erry / y;
0048   return sqrt(pow(xFrErr, 2) + pow(yFrErr, 2) + 2.0 * corr * xFrErr * yFrErr) * x * y;
0049 }
0050 
0051 //=== MAIN MACRO =================================================================================================
0052 
0053 void performFit(string inputDir,
0054                 string OSInputDataFilename,
0055                 string SSInputDataFilename,
0056                 string OSSignalTemplateHistName,
0057                 string SSSignalTemplateHistName) {
0058   gBenchmark->Start("fitZCat");
0059 
0060   //--------------------------------------------------------------------------------------------------------------
0061   // Settings
0062   //==============================================================================================================
0063 
0064   const Double_t mlow = 60;
0065   const Double_t mhigh = 120;
0066   const Int_t nbins = 20;
0067 
0068   // Which fit parameters to set constant
0069   // MuTrk
0070   Bool_t cBgA1MuTrk = kFALSE;
0071   Bool_t cBgA2MuTrk = kTRUE;
0072   Bool_t cBgAlphaMuTrk = kFALSE;
0073   // MuMuNoIso
0074   Bool_t cBgA1MuMuNoIso = kFALSE;
0075   Bool_t cBgAlphaMuMuNoIso = kFALSE;
0076   // MuSta
0077   Bool_t cBgA1MuSta = kFALSE;
0078   Bool_t cBgAlphaMuSta = kFALSE;
0079 
0080   //--------------------------------------------------------------------------------------------------------------
0081   // Main analysis code
0082   //==============================================================================================================
0083   RooRealVar mass("mass", "mass", mlow, mhigh);
0084   mass.setBins(nbins);
0085 
0086   //
0087   // Prepare data for the fits
0088   //
0089   char datname[100];
0090   RooDataSet* dataOS = RooDataSet::read((inputDir + OSInputDataFilename).c_str(), RooArgList(mass));
0091   RooDataSet* dataSS = RooDataSet::read((inputDir + SSInputDataFilename).c_str(), RooArgList(mass));
0092 
0093   //Define categories
0094   RooCategory sample("sample", "");
0095   sample.defineType("OS", 1);
0096   sample.defineType("SS", 2);
0097 
0098   RooDataSet* dataCombined = new RooDataSet(
0099       "dataCombined", "dataCombined", RooArgList(mass), Index(sample), Import("OS", *dataOS), Import("SS", *dataSS));
0100 
0101   //*********************************************************************************************
0102   //Define Free Parameters
0103   //*********************************************************************************************
0104 
0105   RooRealVar* ParNumSignalOS = new RooRealVar("ParNumSignalOS", "ParNumSignalOS", 1000.0, 0.0, 10000000.0);
0106   RooRealVar* ParNumSignalSS = new RooRealVar("ParNumSignalSS", "ParNumSignalSS", 100.0, 0.0, 10000000.0);
0107   RooRealVar* ParNumBkgOS = new RooRealVar("ParNumBkgOS", "ParNumBkgOS", 10.0, 0.0, 10000000.0);
0108   RooRealVar* ParNumBkgSS = new RooRealVar("ParNumBkgSS", "ParNumBkgSS", 1.0, 0.0, 10000000.0);
0109 
0110   RooRealVar* ParOSBackgroundExpCoefficient = new RooRealVar("ParOSBackgroundExpCoefficient",
0111                                                              "ParOSBackgroundExpCoefficient",
0112                                                              -0.01,
0113                                                              -1.0,
0114                                                              1.0);  //ParOSBackgroundExpCoefficient.setConstant(kTRUE);
0115   RooRealVar* ParSSBackgroundExpCoefficient = new RooRealVar("ParSSBackgroundExpCoefficient",
0116                                                              "ParSSBackgroundExpCoefficient",
0117                                                              -0.01,
0118                                                              -1.0,
0119                                                              1.0);  //ParSSBackgroundExpCoefficient.setConstant(kTRUE);
0120   RooRealVar* ParOSSignalMassShift = new RooRealVar("ParOSSignalMassShift", "ParOSSignalMassShift", 0.0, -10.0, 10.0);
0121   ParOSSignalMassShift->setConstant(kTRUE);
0122   RooRealVar* ParSSSignalMassShift = new RooRealVar("ParSSSignalMassShift", "ParSSSignalMassShift", 0.0, -10.0, 10.0);
0123   ParSSSignalMassShift->setConstant(kTRUE);
0124   RooRealVar* ParOSSignalResolution = new RooRealVar(
0125       "ParOSSignalResolution", "ParOSSignalResolution", 1.0, 0.0, 10.0);  //ParOSSignalResolution.setConstant(kTRUE);
0126   RooRealVar* ParSSSignalResolution = new RooRealVar(
0127       "ParSSSignalResolution", "ParSSSignalResolution", 1.0, 0.0, 10.0);  //ParSSSignalResolution.setConstant(kTRUE);
0128 
0129   //*********************************************************************************************
0130   //
0131   //Load Signal PDFs
0132   //
0133   //*********************************************************************************************
0134   //Load histogram templates
0135   TFile* Zeelineshape_file = new TFile("ZMassLineshape_SSOS.root", "READ");
0136   TH1* histTemplateOS = (TH1D*)Zeelineshape_file->Get(OSSignalTemplateHistName.c_str());
0137   TH1* histTemplateSS = (TH1D*)Zeelineshape_file->Get(SSSignalTemplateHistName.c_str());
0138 
0139   //Introduce mass shift coordinate transformation
0140   RooFormulaVar OSShiftedMass("OSShiftedMass", "@0-@1", RooArgList(mass, *ParOSSignalMassShift));
0141   RooFormulaVar SSShiftedMass("SSShiftedMass", "@0-@1", RooArgList(mass, *ParSSSignalMassShift));
0142 
0143   RooDataHist* dataHistOS = new RooDataHist("dataHistOS", "dataHistOS", RooArgSet(mass), histTemplateOS);
0144   RooDataHist* dataHistSS = new RooDataHist("dataHistSS", "dataHistSS", RooArgSet(mass), histTemplateSS);
0145 
0146   RooGaussian* OSSignalResolution =
0147       new RooGaussian("OSSignalResolution", "OSSignalResolution", mass, *ParOSSignalMassShift, *ParOSSignalResolution);
0148   RooGaussian* SSSignalResolution =
0149       new RooGaussian("SSSignalResolution", "SSSignalResolution", mass, *ParSSSignalMassShift, *ParSSSignalResolution);
0150 
0151   RooHistPdf* OSSignalShapeTemplatePdf =
0152       new RooHistPdf("OSSignalShapeTemplatePdf", "OSSignalShapeTemplatePdf", RooArgSet(mass), *dataHistOS, 1);
0153   RooHistPdf* SSSignalShapeTemplatePdf =
0154       new RooHistPdf("SSSignalShapeTemplatePdf", "SSSignalShapeTemplatePdf", RooArgSet(mass), *dataHistSS, 1);
0155 
0156   RooFFTConvPdf* signalShapeOSPdf = new RooFFTConvPdf(
0157       "signalShapeOSPdf", "signalShapeOSPdf", mass, *OSSignalShapeTemplatePdf, *OSSignalResolution, 2);
0158   RooFFTConvPdf* signalShapeSSPdf = new RooFFTConvPdf(
0159       "signalShapeSSPdf", "signalShapeSSPdf", mass, *SSSignalShapeTemplatePdf, *SSSignalResolution, 2);
0160 
0161   //*********************************************************************************************
0162   //
0163   // Create Background PDFs
0164   //
0165   //*********************************************************************************************
0166   RooExponential* bkgOSPdf = new RooExponential("bkgOSPdf", "bkgOSPdf", mass, *ParOSBackgroundExpCoefficient);
0167   RooExponential* bkgSSPdf = new RooExponential("bkgSSPdf", "bkgSSPdf", mass, *ParSSBackgroundExpCoefficient);
0168 
0169   //*********************************************************************************************
0170   //
0171   // Create Total PDFs
0172   //
0173   //*********************************************************************************************
0174   RooAddPdf pdfOS(
0175       "pdfOS", "pdfOS", RooArgList(*signalShapeOSPdf, *bkgOSPdf), RooArgList(*ParNumSignalOS, *ParNumBkgOS));
0176   RooAddPdf pdfSS(
0177       "pdfSS", "pdfSS", RooArgList(*signalShapeSSPdf, *bkgSSPdf), RooArgList(*ParNumSignalSS, *ParNumBkgSS));
0178 
0179   // PDF for simultaneous fit
0180   RooSimultaneous pdfTotal("pdfTotal", "pdfTotal", sample);
0181   pdfTotal.addPdf(pdfOS, "pdfOS");
0182   pdfTotal.addPdf(pdfSS, "pdfSS");
0183   pdfTotal.Print();
0184 
0185   //
0186   // Define likelihood, add constraints, and run the fit
0187   //
0188   RooFitResult* fitResult =
0189       pdfTotal.fitTo(*dataCombined, RooFit::Save(true), RooFit::Extended(true), RooFit::PrintLevel(-1));
0190   fitResult->Print("v");
0191 
0192   double nSignalOS = ParNumSignalOS->getVal();
0193   double nSignalSS = ParNumSignalSS->getVal();
0194   double nBkgOS = ParNumBkgOS->getVal();
0195   double nBkgSS = ParNumBkgSS->getVal();
0196 
0197   printf("\nFit results:\n");
0198   cout << "Signal OS: " << nSignalOS << endl;
0199   cout << "Bkg OS: " << nBkgOS << endl;
0200   cout << "Signal SS: " << nSignalSS << endl;
0201   cout << "Bkg SS: " << nBkgSS << endl;
0202 
0203   //--------------------------------------------------------------------------------------------------------------
0204   // Make plots
0205   //==============================================================================================================
0206   RooAbsData::ErrorType errorType = RooAbsData::Poisson;
0207 
0208   TString cname = TString("fit_OS");
0209   TCanvas* c = new TCanvas(cname, cname, 800, 600);
0210   RooPlot* frame1 = mass.frame();
0211   frame1->SetMinimum(0);
0212   dataOS->plotOn(frame1, RooFit::DataError(errorType));
0213   pdfOS.plotOn(frame1, RooFit::ProjWData(*dataOS), RooFit::Components(*bkgOSPdf), RooFit::LineColor(kRed));
0214   pdfOS.plotOn(frame1, RooFit::ProjWData(*dataOS));
0215   frame1->Draw("e0");
0216 
0217   TPaveText* plotlabel = new TPaveText(0.23, 0.87, 0.43, 0.92, "NDC");
0218   plotlabel->SetTextColor(kBlack);
0219   plotlabel->SetFillColor(kWhite);
0220   plotlabel->SetBorderSize(0);
0221   plotlabel->SetTextAlign(12);
0222   plotlabel->SetTextSize(0.03);
0223   plotlabel->AddText("CMS Preliminary 2010");
0224   TPaveText* plotlabel2 = new TPaveText(0.23, 0.82, 0.43, 0.87, "NDC");
0225   plotlabel2->SetTextColor(kBlack);
0226   plotlabel2->SetFillColor(kWhite);
0227   plotlabel2->SetBorderSize(0);
0228   plotlabel2->SetTextAlign(12);
0229   plotlabel2->SetTextSize(0.03);
0230   plotlabel2->AddText("#sqrt{s} = 7 TeV");
0231   TPaveText* plotlabel3 = new TPaveText(0.23, 0.75, 0.43, 0.80, "NDC");
0232   plotlabel3->SetTextColor(kBlack);
0233   plotlabel3->SetFillColor(kWhite);
0234   plotlabel3->SetBorderSize(0);
0235   plotlabel3->SetTextAlign(12);
0236   plotlabel3->SetTextSize(0.03);
0237   char temp[100];
0238   sprintf(temp, "%.4f", LUMINOSITY);
0239   plotlabel3->AddText((string("#int#font[12]{L}dt = ") + temp + string(" pb^{ -1}")).c_str());
0240   TPaveText* plotlabel4 = new TPaveText(0.6, 0.82, 0.8, 0.87, "NDC");
0241   plotlabel4->SetTextColor(kBlack);
0242   plotlabel4->SetFillColor(kWhite);
0243   plotlabel4->SetBorderSize(0);
0244   plotlabel4->SetTextAlign(12);
0245   plotlabel4->SetTextSize(0.03);
0246   sprintf(temp, "Signal = %.2f #pm %.2f", ParNumSignalOS->getVal(), ParNumSignalOS->getError());
0247   plotlabel4->AddText(temp);
0248   TPaveText* plotlabel5 = new TPaveText(0.6, 0.77, 0.8, 0.82, "NDC");
0249   plotlabel5->SetTextColor(kBlack);
0250   plotlabel5->SetFillColor(kWhite);
0251   plotlabel5->SetBorderSize(0);
0252   plotlabel5->SetTextAlign(12);
0253   plotlabel5->SetTextSize(0.03);
0254   sprintf(temp, "Bkg = %.2f #pm %.2f", ParNumBkgOS->getVal(), ParNumBkgOS->getError());
0255   plotlabel5->AddText(temp);
0256   plotlabel4->Draw();
0257   plotlabel5->Draw();
0258 
0259   //   c->SaveAs( cname + TString(".eps"));
0260   c->SaveAs(cname + TString(".gif"));
0261   delete c;
0262 
0263   cname = TString("fit_SS");
0264   TCanvas* c2 = new TCanvas(cname, cname, 500, 500);
0265   RooPlot* frame2 = mass.frame();
0266   frame2->SetMinimum(0);
0267   dataSS->plotOn(frame2, RooFit::DataError(errorType));
0268   pdfSS.plotOn(frame2, RooFit::ProjWData(*dataSS), RooFit::Components(*bkgSSPdf), RooFit::LineColor(kRed));
0269   pdfSS.plotOn(frame2, RooFit::ProjWData(*dataSS));
0270   frame2->Draw("e0");
0271 
0272   plotlabel = new TPaveText(0.23, 0.87, 0.43, 0.92, "NDC");
0273   plotlabel->SetTextColor(kBlack);
0274   plotlabel->SetFillColor(kWhite);
0275   plotlabel->SetBorderSize(0);
0276   plotlabel->SetTextAlign(12);
0277   plotlabel->SetTextSize(0.03);
0278   plotlabel->AddText("CMS Preliminary 2010");
0279   plotlabel2 = new TPaveText(0.23, 0.82, 0.43, 0.87, "NDC");
0280   plotlabel2->SetTextColor(kBlack);
0281   plotlabel2->SetFillColor(kWhite);
0282   plotlabel2->SetBorderSize(0);
0283   plotlabel2->SetTextAlign(12);
0284   plotlabel2->SetTextSize(0.03);
0285   plotlabel2->AddText("#sqrt{s} = 7 TeV");
0286   plotlabel3 = new TPaveText(0.23, 0.75, 0.43, 0.80, "NDC");
0287   plotlabel3->SetTextColor(kBlack);
0288   plotlabel3->SetFillColor(kWhite);
0289   plotlabel3->SetBorderSize(0);
0290   plotlabel3->SetTextAlign(12);
0291   plotlabel3->SetTextSize(0.03);
0292   sprintf(temp, "%.4f", LUMINOSITY);
0293   plotlabel3->AddText((string("#int#font[12]{L}dt = ") + temp + string(" pb^{ -1}")).c_str());
0294   plotlabel4 = new TPaveText(0.6, 0.82, 0.8, 0.87, "NDC");
0295   plotlabel4->SetTextColor(kBlack);
0296   plotlabel4->SetFillColor(kWhite);
0297   plotlabel4->SetBorderSize(0);
0298   plotlabel4->SetTextAlign(12);
0299   plotlabel4->SetTextSize(0.03);
0300   sprintf(temp, "Signal = %.2f #pm %.2f", ParNumSignalSS->getVal(), ParNumSignalSS->getError());
0301   plotlabel4->AddText(temp);
0302   plotlabel5 = new TPaveText(0.6, 0.77, 0.8, 0.82, "NDC");
0303   plotlabel5->SetTextColor(kBlack);
0304   plotlabel5->SetFillColor(kWhite);
0305   plotlabel5->SetBorderSize(0);
0306   plotlabel5->SetTextAlign(12);
0307   plotlabel5->SetTextSize(0.03);
0308   sprintf(temp, "Bkg = %.2f #pm %.2f", ParNumBkgSS->getVal(), ParNumBkgSS->getError());
0309   plotlabel5->AddText(temp);
0310 
0311   //   plotlabel->Draw();
0312   //   plotlabel2->Draw();
0313   //   plotlabel3->Draw();
0314   plotlabel4->Draw();
0315   plotlabel5->Draw();
0316 
0317   c2->SaveAs(cname + TString(".eps"));
0318   c2->SaveAs(cname + TString(".gif"));
0319   c2->SaveAs(cname + TString(".root"));
0320   delete c2;
0321 
0322   gBenchmark->Show("fitZCat");
0323 }
0324 
0325 void ZeeBkgFit() {
0326   performFit("MitFitter/inputs/BkgFit/",
0327              "M_VBTF95PlusVBTF95_OS",
0328              "M_VBTF95PlusVBTF95_SS",
0329              "Mass_WP95WP95_OppositeSign",
0330              "Mass_WP95WP95_SameSign");
0331 }