test/utilities/TwoCategorySimZFitter.C

0001 /*******************************************************************
0002  * Project: CMS detector at the CERN
0003  *
0004  * Package: Presently in the user code
0005  *
0006  *
0007  * Authors:
0008  *
0009  *   Kalanand Mishra, Fermilab - kalanand@fnal.gov
0010  *
0011  * Description:
0012  *   A standalone macro to perform single bin simultaneos fit of
0013  *   dilepton mass under Z peak.
0014  *   Works with  unbinned data.
0015  *
0016  * Implementation details:
0017  *  Uses RooFit classes.
0018  *
0019  * History:
0020  *
0021  *
0022  ********************************************************************/
0023
0024 //// Following are the 2 irreducible free parameters of the fit:
0025 ////        Z cross section and single electron efficiency.
0026 ////  Additionally, the following 2 nuisance parameters are floating:
0027 ////        nBkgTF, bkgGamma.
0028
0029
0030 // ROOT
0031 #include <TROOT.h>
0032 #include <TChain.h>
0033 #include <TFile.h>
0034 #include <TCanvas.h>
0035 #include <TH1F.h>
0036 #include <TH2F.h>
0037 #include <TTree.h>
0038 #include <TGraph.h>
0039 #include "tdrstyle.C"
0040
0041 using namespace RooFit;
0042
0043 void makeSignalPdf();
0044 void makeBkgPdf();
0045
0046 // The signal & background Pdf
0047 RooRealVar *rooMass_;
0048 RooAbsPdf* signalShapePdfTT_;
0049 RooAbsPdf* signalShapePdfTF_;
0050 RooAbsPdf *bkgShapePdfTF_;
0051
0052
0053 TFile* Zeelineshape_file;
0054 TCanvas *c;
0055
0056 // Acceptance for Gen-->SC
0057 const float A_BB= 0.2257;
0058 const float A_BE= 0.1612;
0059 const float A_EE= 0.0476;
0060
0061 // Specify the electron selection: "WP95" or "WP80"
0062 const char* selection = "WP95";
0063
0064
0065 // Integrated luminosity
0066 const double intLumi = 2.88;
0067
0068
0069 // In case we want to fix the nuisance parameters ...
0070 const bool FIX_NUIS_PARS=false;
0071
0072
0073
0074
0075 void TwoCategorySimZFitter()
0076 {
0077   // The fit variable - lepton invariant mass
0078   rooMass_ = new RooRealVar("mass","m_{ee}",60.0, 120.0, "GeV/c^{2}");
0079   rooMass_->setBins(120.0);
0080   RooRealVar Mass = *rooMass_;
0081
0082   // Make the category variable that defines the two fits,
0083   // namely whether the probe passes or fails the eff criteria.
0084   RooCategory sample("sample","") ;
0085   sample.defineType("pass", 1) ;
0086   sample.defineType("fail", 0) ;
0087
0088
0089   gROOT->cd();
0090   char inFile[50];
0091   sprintf(inFile, "ZeeEvents_%s.root", selection);
0092   TFile fin(inFile, "read");
0093   TTree* treeTT = (TTree*) fin.Get("tree");
0094   sprintf(inFile, "ZeeEvents_%s_TF.root", selection);
0095   TFile fin2(inFile, "read");
0096   TTree* treeTF = (TTree*) fin2.Get("tree");
0097
0098
0099   ///////// convert Histograms into RooDataHists
0100   RooDataSet* data_TT = new RooDataSet("data_TT","data_TT", treeTT, Mass);
0101   RooDataSet* data_TF = new RooDataSet("data_TF_BB","data_TF",treeTF, Mass);
0102
0103   RooDataSet* data = new RooDataSet( "data","data",
0104                        RooArgList(Mass),Index(sample),
0105                        Import("pass",*data_TT),
0106                        Import("fail",*data_TF) );
0107
0108   data->get()->Print();
0109   cout << "Made dataset" << endl;
0110
0111
0112
0113
0114   // ********** Construct signal & bkg shape PDF ********** //
0115   makeSignalPdf();
0116   cout << "Made signal pdf" << endl;
0117   makeBkgPdf();
0118   cout << "Made bkg pdf" << endl;
0119
0120   // Now supply integrated luminosity in inverse picobarn
0121   RooRealVar lumi("lumi","lumi", intLumi);
0122
0123
0124   // Now define Z production cross section variable (in pb)
0125   RooRealVar xsec("xsec","xsec", 912.187, 200.0, 2000.0);
0126
0127
0128   // Define efficiency variable: single bin efficiency, ignore trigger for now
0129   RooRealVar eff("eff","eff", 0.93, 0.0, 1.0);
0130
0131
0132   // Now define acceptance variables-->we get these numbers from MC
0133   RooRealVar acc("acc","acc", A_BB+A_BE+A_EE);
0134
0135
0136
0137   // Define background yield variables: they are not related to each other
0138   float numBkgHighPurity=11.8;
0139   if(!strcmp(selection,"WP80")) numBkgHighPurity=3.0;
0140   RooRealVar nBkgTT("nBkgTT","nBkgTT", numBkgHighPurity);
0141   RooRealVar nBkgTF("nBkgTF","nBkgTF", 58.0,     0.0, 500.);
0142   if(!strcmp(selection,"WP80")) {
0143      nBkgTF.setVal(0.0);
0144      nBkgTF.setConstant(true);
0145   }
0146    if(FIX_NUIS_PARS) nBkgTF.setConstant(true);
0147
0148   ////////////////////////////////////////////////////////////////////////////
0149   ////////////////////////////////////////////////////////////////////////////
0150  //  Define signal yield variables.
0151   // They are linked together by the total cross section:  e.g.
0152   //          Nbb = sigma*L*Abb*effB
0153
0154   const char* formula = 0;
0155   RooArgList* args;
0156   formula="lumi*xsec*acc*eff*eff+nBkgTT";
0157   args = new RooArgList(lumi,xsec,acc,eff,nBkgTT);
0158   RooFormulaVar nSigTT("nSigTT", formula, *args);
0159   delete args;
0160
0161   formula="lumi*xsec*acc*eff*(1.0-eff)";
0162   args = new RooArgList(lumi,xsec,acc,eff);
0163   RooFormulaVar nSigTF("nSigTF",formula, *args);
0164   delete args;
0165
0166   /////////////////////////////////////////////////////////////////////////
0167   /////////////////////////////////////////////////////////////////////////
0168
0169    RooArgList componentsTF(*signalShapePdfTF_,*bkgShapePdfTF_);
0170    RooArgList yieldsTF(nSigTF, nBkgTF);
0171    RooExtendPdf pdfTT("pdfTT","extended sum pdf", *signalShapePdfTT_, nSigTT);
0172    RooAddPdf pdfTF("pdfTF","extended sum pdf",componentsTF, yieldsTF);
0173
0174
0175    // The total simultaneous fit ...
0176    RooSimultaneous totalPdf("totalPdf","totalPdf", sample);
0177    totalPdf.addPdf(pdfTT,"pass");
0178    totalPdf.Print();
0179    totalPdf.addPdf(pdfTF,"fail");
0180    totalPdf.Print();
0181
0182
0183   // ********* Do the Actual Fit ********** //
0184    RooFitResult *fitResult = totalPdf.fitTo(*data, Save(true),
0185    RooFit::Extended(true),
0186    RooFit::Minos(true),
0187    PrintEvalErrors(-1),
0188    Warnings(false)
0189       );
0190   fitResult->Print("v");
0191
0192   // ********** Make and save Canvas for the plots ********** //
0193   gROOT->ProcessLine(".L ~/tdrstyle.C");
0194   setTDRStyle();
0195   tdrStyle->SetErrorX(0.5);
0196   tdrStyle->SetPadLeftMargin(0.19);
0197   tdrStyle->SetPadRightMargin(0.10);
0198   tdrStyle->SetPadBottomMargin(0.15);
0199   tdrStyle->SetLegendBorderSize(0);
0200   tdrStyle->SetTitleYOffset(1.5);
0201   RooAbsData::ErrorType errorType = RooAbsData::SumW2;
0202
0203
0204   TString cname = Form("Zmass_TT_%dnb", (int)(1000*intLumi) );
0205   c = new TCanvas(cname,cname,500,500);
0206   RooPlot* frame1 = Mass.frame(60., 120., 60);
0207   data_TT->plotOn(frame1,RooFit::DataError(errorType));
0208   pdfTT.plotOn(frame1,ProjWData(*data_TT));
0209   frame1->SetMinimum(0);
0210   frame1->Draw("e0");
0211   TPaveText *plotlabel = new TPaveText(0.23,0.87,0.43,0.92,"NDC");
0212    plotlabel->SetTextColor(kBlack);
0213    plotlabel->SetFillColor(kWhite);
0214    plotlabel->SetBorderSize(0);
0215    plotlabel->SetTextAlign(12);
0216    plotlabel->SetTextSize(0.03);
0217    plotlabel->AddText("CMS Preliminary 2010");
0218   TPaveText *plotlabel2 = new TPaveText(0.23,0.82,0.43,0.87,"NDC");
0219    plotlabel2->SetTextColor(kBlack);
0220    plotlabel2->SetFillColor(kWhite);
0221    plotlabel2->SetBorderSize(0);
0222    plotlabel2->SetTextAlign(12);
0223    plotlabel2->SetTextSize(0.03);
0224    plotlabel2->AddText("#sqrt{s} = 7 TeV");
0225   TPaveText *plotlabel3 = new TPaveText(0.23,0.75,0.43,0.80,"NDC");
0226    plotlabel3->SetTextColor(kBlack);
0227    plotlabel3->SetFillColor(kWhite);
0228    plotlabel3->SetBorderSize(0);
0229    plotlabel3->SetTextAlign(12);
0230    plotlabel3->SetTextSize(0.03);
0231   char temp[100];
0232   sprintf(temp, "%.1f", intLumi);
0233   plotlabel3->AddText((string("#int#font[12]{L}dt = ") +
0234   temp + string(" pb^{ -1}")).c_str());
0235   TPaveText *plotlabel4 = new TPaveText(0.6,0.87,0.8,0.92,"NDC");
0236    plotlabel4->SetTextColor(kBlack);
0237    plotlabel4->SetFillColor(kWhite);
0238    plotlabel4->SetBorderSize(0);
0239    plotlabel4->SetTextAlign(12);
0240    plotlabel4->SetTextSize(0.03);
0241    double nsig = nSigTT.getVal();
0242    double nsigerr = nSigTT.getPropagatedError(*fitResult) ;
0243    sprintf(temp, "Signal = %.1f #pm %.1f", nsig, nsigerr);
0244    plotlabel4->AddText(temp);
0245    TPaveText *plotlabel5 = new TPaveText(0.6,0.82,0.8,0.87,"NDC");
0246    plotlabel5->SetTextColor(kBlack);
0247    plotlabel5->SetFillColor(kWhite);
0248    plotlabel5->SetBorderSize(0);
0249    plotlabel5->SetTextAlign(12);
0250    plotlabel5->SetTextSize(0.03);
0251    sprintf(temp, "#epsilon = %.3f #pm %.3f", eff.getVal(), eff.getError() );
0252    plotlabel5->AddText(temp);
0253   TPaveText *plotlabel7 = new TPaveText(0.6,0.77,0.8,0.82,"NDC");
0254    plotlabel7->SetTextColor(kBlack);
0255    plotlabel7->SetFillColor(kWhite);
0256    plotlabel7->SetBorderSize(0);
0257    plotlabel7->SetTextAlign(12);
0258    plotlabel7->SetTextSize(0.03);
0259    sprintf(temp, "#sigma = %.1f #pm %.1f pb", xsec.getVal(), xsec.getError());
0260    plotlabel7->AddText(temp);
0261   plotlabel->Draw();
0262   plotlabel2->Draw();
0263   plotlabel3->Draw();
0264   plotlabel4->Draw();
0265   plotlabel5->Draw();
0266   plotlabel7->Draw();
0267 //   c->SaveAs( cname + TString(".eps"));
0268 //   c->SaveAs( cname + TString(".gif"));
0269 //   c->SaveAs( cname + TString(".root"));
0270 //   c->SaveAs( cname + TString(".png"));
0271 //   c->SaveAs( cname + TString(".C"));
0272
0273
0274   cname = Form("Zmass_TF%dnb", (int)(1000*intLumi) );
0275   c = new TCanvas(cname,cname,500,500);
0276   RooPlot* frame2 = Mass.frame(60., 120., 12);
0277   data_TF->plotOn(frame2,RooFit::DataError(errorType));
0278   pdfTF.plotOn(frame2,ProjWData(*data_TF),
0279   Components(*signalShapePdfTF_),DrawOption("LF"),FillStyle(1001),FillColor(kOrange),VLines());
0280   pdfTF.plotOn(frame2,ProjWData(*data_TF));
0281   pdfTF.plotOn(frame2,ProjWData(*data_TF),
0282   Components(*bkgShapePdfTF_),LineColor(kRed));
0283   data_TF->plotOn(frame2,RooFit::DataError(errorType));
0284   frame2->SetMinimum(0);
0285   frame2->Draw("e0");
0286   frame2->GetYaxis()->SetNdivisions(505);
0287   plotlabel = new TPaveText(0.23,0.87,0.43,0.92,"NDC");
0288    plotlabel->SetTextColor(kBlack);
0289    plotlabel->SetFillColor(kWhite);
0290    plotlabel->SetBorderSize(0);
0291    plotlabel->SetTextAlign(12);
0292    plotlabel->SetTextSize(0.03);
0293    plotlabel->AddText("CMS Preliminary 2010");
0294   plotlabel2 = new TPaveText(0.23,0.82,0.43,0.87,"NDC");
0295    plotlabel2->SetTextColor(kBlack);
0296    plotlabel2->SetFillColor(kWhite);
0297    plotlabel2->SetBorderSize(0);
0298    plotlabel2->SetTextAlign(12);
0299    plotlabel2->SetTextSize(0.03);
0300    plotlabel2->AddText("#sqrt{s} = 7 TeV");
0301   plotlabel3 = new TPaveText(0.23,0.75,0.43,0.80,"NDC");
0302    plotlabel3->SetTextColor(kBlack);
0303    plotlabel3->SetFillColor(kWhite);
0304    plotlabel3->SetBorderSize(0);
0305    plotlabel3->SetTextAlign(12);
0306    plotlabel3->SetTextSize(0.03);
0307   char temp2[100];
0308   sprintf(temp2, "%.1f", intLumi);
0309   plotlabel3->AddText((string("#int#font[12]{L}dt = ") +
0310   temp2 + string(" pb^{ -1}")).c_str());
0311   plotlabel4 = new TPaveText(0.6,0.87,0.8,0.92,"NDC");
0312    plotlabel4->SetTextColor(kBlack);
0313    plotlabel4->SetFillColor(kWhite);
0314    plotlabel4->SetBorderSize(0);
0315    plotlabel4->SetTextAlign(12);
0316    plotlabel4->SetTextSize(0.03);
0317    nsig = nSigTF.getVal();
0318    nsigerr = nSigTF.getPropagatedError(*fitResult) ;
0319    sprintf(temp2, "Signal = %.2f #pm %.2f", nsig, nsigerr);
0320    plotlabel4->AddText(temp2);
0321   plotlabel5 = new TPaveText(0.6,0.82,0.8,0.87,"NDC");
0322    plotlabel5->SetTextColor(kBlack);
0323    plotlabel5->SetFillColor(kWhite);
0324    plotlabel5->SetBorderSize(0);
0325    plotlabel5->SetTextAlign(12);
0326    plotlabel5->SetTextSize(0.03);
0327    sprintf(temp2, "Bkg = %.2f #pm %.2f", nBkgTF.getVal(), nBkgTF.getError());
0328    plotlabel5->AddText(temp2);
0329   TPaveText* plotlabel6 = new TPaveText(0.6,0.77,0.8,0.82,"NDC");
0330    plotlabel6->SetTextColor(kBlack);
0331    plotlabel6->SetFillColor(kWhite);
0332    plotlabel6->SetBorderSize(0);
0333    plotlabel6->SetTextAlign(12);
0334    plotlabel6->SetTextSize(0.03);
0335    sprintf(temp2, "#epsilon = %.3f #pm %.3f", eff.getVal(), eff.getError() );
0336    plotlabel6->AddText(temp2);
0337   plotlabel->Draw();
0338   plotlabel2->Draw();
0339   plotlabel3->Draw();
0340   plotlabel4->Draw();
0341   plotlabel5->Draw();
0342   plotlabel6->Draw();
0343   plotlabel7->Draw();
0344
0345 //   c->SaveAs( cname + TString(".eps"));
0346 //   c->SaveAs( cname + TString(".gif"));
0347 //   c->SaveAs( cname + TString(".root"));
0348 //   c->SaveAs( cname + TString(".png"));
0349 //   c->SaveAs( cname + TString(".C"));
0350
0351
0352   //    if(data) delete data;
0353   //    if(c) delete c;
0354 }
0355
0356
0357
0358
0359
0360 // // ***** Function to return the signal Pdf *** //
0361  void makeSignalPdf() {
0362
0363  // Tag+Tag selection pdf
0364   Zeelineshape_file =  new TFile("../Zlineshapes.root", "READ");
0365
0366   // Tag+Tag selection pdf
0367   TH1* histbbpass = (TH1D*) Zeelineshape_file->Get("pass_BB");
0368   TH1* histebpass = (TH1D*) Zeelineshape_file->Get("pass_BE");
0369   TH1* histeepass = (TH1D*) Zeelineshape_file->Get("pass_EE");
0370
0371   TH1D* th1_TT = (TH1D*) histbbpass->Clone("th1_TT");
0372   th1_TT->Add(histebpass);
0373   th1_TT->Add(histeepass);
0374
0375   RooRealVar* zero  = new RooRealVar("zero","", 0.0);
0376   RooDataHist* rdh_TT = new RooDataHist("rdh_TT","", *rooMass_, th1_TT);
0377   RooAbsPdf* signalModelTT_ = new RooHistPdf("signalModelTT", "", *rooMass_, *rdh_TT);
0378   RooRealVar* resoTT_  = new RooRealVar("resoTT","resoTT",1.88, 0.0, 5.);
0379   if(FIX_NUIS_PARS) resoTT_->setConstant(true);
0380   RooGaussModel* resModelTT_ = new RooGaussModel("resModelTT","",
0381                 *rooMass_, *zero, *resoTT_);
0382   signalShapePdfTT_  = new RooFFTConvPdf("sigModel","final signal shape",
0383                        *rooMass_, *signalModelTT_, *resModelTT_);
0384
0385
0386  // Tag+Fail selection pdf
0387   TH1* histbb = (TH1D*) Zeelineshape_file->Get("fail_BB");
0388   TH1* histeb = (TH1D*) Zeelineshape_file->Get("fail_BE");
0389   TH1* histee = (TH1D*) Zeelineshape_file->Get("fail_EE");
0390
0391   TH1D* th1_TF = (TH1D*) histbb->Clone("th1_TF");
0392   th1_TF->Add(histeb);
0393   th1_TF->Add(histee);
0394   RooDataHist* rdh_TF_ = new RooDataHist("rdh_TF_BB","", *rooMass_, th1_TF);
0395   RooAbsPdf* signalModelTF_ = new RooHistPdf("signalModelTF",
0396   "",*rooMass_,*rdh_TF_);
0397   RooRealVar* resoTF_  = new RooRealVar("resoTF","resoTF", 2.76762, 0.0, 10.);
0398   if(FIX_NUIS_PARS) resoTF_->setConstant(true);
0399   RooGaussModel* resModelTF_ = new RooGaussModel("resModelTF",
0400   "gaussian resolution model", *rooMass_, *zero, *resoTF_);
0401   signalShapePdfTF_  = new RooFFTConvPdf("signalShapePdfTF","",
0402                        *rooMass_, *signalModelTF_, *resModelTF_);
0403 }
0404
0405
0406
0407
0408 // ***** Function to return the background Pdf **** //
0409 void makeBkgPdf()
0410 {
0411   // Background PDF variables
0412    RooRealVar* bkgGammaFailTF_ = new RooRealVar("bkgGammaFailTF",
0413    "",-0.028394, -10., 10.);
0414    if(FIX_NUIS_PARS) bkgGammaFailTF_->setConstant(true);
0415    bkgShapePdfTF_ = new RooExponential("bkgShapePdfTF",
0416    "",*rooMass_, *bkgGammaFailTF_);
0417 }
0418
0419
0420
0421
0422
0423
0424
0425
0426