File indexing completed on 2024-04-06 12:24:14
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
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
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
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
0061
0062
0063 const Double_t mlow = 60;
0064 const Double_t mhigh = 120;
0065 const Int_t nbins = 20;
0066
0067
0068
0069 Bool_t cBgA1MuTrk = kFALSE;
0070 Bool_t cBgA2MuTrk = kTRUE;
0071 Bool_t cBgAlphaMuTrk = kFALSE;
0072
0073 Bool_t cBgA1MuMuNoIso = kFALSE;
0074 Bool_t cBgAlphaMuMuNoIso = kFALSE;
0075
0076 Bool_t cBgA1MuSta = kFALSE;
0077 Bool_t cBgAlphaMuSta = kFALSE;
0078
0079
0080
0081
0082 RooRealVar mass("mass", "mass", mlow, mhigh);
0083 mass.setBins(nbins);
0084
0085
0086
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
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
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);
0114 RooRealVar* ParSSBackgroundExpCoefficient = new RooRealVar("ParSSBackgroundExpCoefficient",
0115 "ParSSBackgroundExpCoefficient",
0116 -0.01,
0117 -1.0,
0118 1.0);
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);
0125 RooRealVar* ParSSSignalResolution = new RooRealVar(
0126 "ParSSSignalResolution", "ParSSSignalResolution", 1.0, 0.0, 10.0);
0127
0128
0129
0130
0131
0132
0133
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
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
0163
0164
0165 RooExponential* bkgOSPdf = new RooExponential("bkgOSPdf", "bkgOSPdf", mass, *ParOSBackgroundExpCoefficient);
0166 RooExponential* bkgSSPdf = new RooExponential("bkgSSPdf", "bkgSSPdf", mass, *ParSSBackgroundExpCoefficient);
0167
0168
0169
0170
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
0179 RooSimultaneous pdfTotal("pdfTotal", "pdfTotal", sample);
0180 pdfTotal.addPdf(pdfOS, "pdfOS");
0181 pdfTotal.addPdf(pdfSS, "pdfSS");
0182 pdfTotal.Print();
0183
0184
0185
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
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
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
0311
0312
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 }
