File indexing completed on 2024-04-06 12:24:13
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
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
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
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
0060
0061
0062 const Double_t mlow = 60;
0063 const Double_t mhigh = 120;
0064 const Int_t nbins = 20;
0065
0066
0067
0068 Bool_t cBgA1MuTrk = kFALSE;
0069 Bool_t cBgA2MuTrk = kTRUE;
0070 Bool_t cBgAlphaMuTrk = kFALSE;
0071
0072 Bool_t cBgA1MuMuNoIso = kFALSE;
0073 Bool_t cBgAlphaMuMuNoIso = kFALSE;
0074
0075 Bool_t cBgA1MuSta = kFALSE;
0076 Bool_t cBgAlphaMuSta = kFALSE;
0077
0078
0079
0080
0081 RooRealVar mass("mass", "mass", mlow, mhigh);
0082 mass.setBins(nbins);
0083
0084
0085
0086
0087
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
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);
0117 RooRealVar* ParFailBackgroundExpCoefficient =
0118 new RooRealVar("ParFailBackgroundExpCoefficient",
0119 "ParFailBackgroundExpCoefficient",
0120 -0.2,
0121 -1.0,
0122 0.0);
0123 RooRealVar* ParPassSignalMassShift = new RooRealVar(
0124 "ParPassSignalMassShift", "ParPassSignalMassShift", -10.0, 10.0);
0125 RooRealVar* ParFailSignalMassShift = new RooRealVar(
0126 "ParFailSignalMassShift", "ParFailSignalMassShift", -10.0, 10.0);
0127 RooRealVar* ParPassSignalResolution = new RooRealVar(
0128 "ParPassSignalResolution", "ParPassSignalResolution", 0.0, 10.0);
0129 RooRealVar* ParFailSignalResolution = new RooRealVar(
0130 "ParFailSignalResolution", "ParFailSignalResolution", 0.0, 10.0);
0131
0132
0133
0134
0135
0136
0137
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
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
0159
0160
0161 RooExponential* bkgPassPdf = new RooExponential("bkgPassPdf", "bkgPassPdf", mass, *ParPassBackgroundExpCoefficient);
0162 RooExponential* bkgFailPdf = new RooExponential("bkgFailPdf", "bkgFailPdf", mass, *ParFailBackgroundExpCoefficient);
0163
0164
0165
0166
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
0175 RooSimultaneous pdfTotal("pdfTotal", "pdfTotal", sample);
0176 pdfTotal.addPdf(pdfPass, "pdfPass");
0177 pdfTotal.addPdf(pdfFail, "pdfFail");
0178 pdfTotal.Print();
0179
0180
0181
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
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
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
0345
0346
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 }
