Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-01-08 23:50:50

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