NonMCBackground/src/plot_emu_fakes.cc

#include <iostream>
#include <algorithm>
#include <iomanip>

#include "RooHistPdf.h"
#include "RooGlobalFunc.h"
#include "RooPlot.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooGenericPdf.h"
#include "RooDataSet.h"
#include "RooAddPdf.h"
#include "RooFormulaVar.h"
#include "RooFitResult.h"
#include "RooCurve.h"
#include "RooBinning.h"
#include "RooExponential.h"
#include "RooExtendPdf.h"

#include "TCanvas.h"
#include "TChain.h"
#include "TString.h"
#include "TStyle.h"
#include "TH1D.h"

#include "CPlot.h"
#include "FOArgs.h"

#include "fake_defs.h"
#include "Various.h"
#include "PlotHeaders.h"

#ifndef CMSSW_BASE
#define CMSSW_BASE "../../"
#endif

using namespace std;
using namespace RooFit;
using namespace mithep;

TCanvas *can;

void makeHTML(FOFlags &ctrl, TString plotLabel);
map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str="");
void fillHist(CSample *cs, TString type, TString var, double val, double wgt, double wgt_lo, double wgt_hi);

RooFitResult *fitMt(FOFlags &ctrl, TH1D *hmcZ, TH1D *hmcW, TH1D *hSSmet, double &total_mc_scale, TString plotLabel);

//----------------------------------------------------------------------------------------
int main(int argc, char** argv) 
{
  double lumi = 1600;

  FOFlags ctrl;
  parse_foargs( argc, argv, ctrl );
  ctrl.dump();

  can = new TCanvas("can","can");

  FR_struct fr2011 = initFRs(ctrl.mufakefile2011,ctrl.elefakefile2011);
  FR_struct fr2012 = initFRs(ctrl.mufakefile2012,ctrl.elefakefile2012);

  TString cmsswpath(CMSSW_BASE + TString("/src"));
  string xspath = (string(cmsswpath)+"/MitPhysics/data/xs.dat");
  SimpleTable xstab(xspath.c_str());

  vector<CSample*> samplev;
  TString ntupledir(""),label(""),plotLabel(""),jsonFile("");
  readConfigFile(ctrl.config, ntupledir, label, plotLabel, jsonFile, samplev, &ctrl, init_hists);

  vector<pair<unsigned,unsigned> > runEvtv; // vector to veto duplicate events
  UInt_t minRun=999999999,maxRun=0;
  for(unsigned ics=0; ics<samplev.size(); ics++) {
    CSample *cs = samplev[ics];
    cout << cs->name << ": " << endl;
    for(unsigned ifs=0; ifs<(cs->fsv).size(); ifs++) {
      filestuff *fs = (cs->fsv)[ifs];
      cout << "\t" << fs->fname_ << endl;
      FR_struct *fr = (fs->era_==2011) ? &fr2011 : &fr2012;
      unsigned nDuplSkipped=0;
      for(unsigned ientry=0; ientry<fs->getentries("FOtree"); ientry++) {
        fs->getentry(ientry,"info");

        if(fs->isdata_) {
          setMinMaxRun(fs->info->run, minRun, maxRun);
          // bool dupl = takeCareOfDuplicateEvents(fs->info->run, fs->info->evt, runEvtv, nDuplSkipped);
          // if(dupl) continue;
          pair<unsigned,unsigned> rl(fs->info->run, fs->info->lumi);
          if(!fs->rlrm_.HasRunLumi(rl)) continue;
        }

        unsigned npass = fs->passingL->size();
        unsigned nfail = fs->failingL->size();

        // z veto
        bool zVeto=false;
        for(unsigned ilep1=0; ilep1<npass; ilep1++) {
          SimpleLepton lep1 = (*fs->passingL)[ilep1];
          for(unsigned ilep2=ilep1+1; ilep2<npass; ilep2++) {
            SimpleLepton lep2 = (*fs->passingL)[ilep2];
            // if(abs(lep1.type) == abs(lep2.type)) zVeto = true;
            // if(lep1.charge != lep2.charge)       zVeto = true;
            if(abs(lep1.type) == abs(lep2.type) &&
               lep1.charge != lep2.charge)       zVeto = true;
          }
        }
        if(zVeto) continue;

        for(unsigned iw=0; iw<npass; iw++) {
          SimpleLepton w_lep = (*fs->passingL)[iw];

          double wgt=1;
          if(!fs->isdata_) {
            double xsWgt = lumi*xstab.Get(fs->dataset_)/fs->total_entries_;
            double puWgt = weightTrue2012(fs->info->npu);
            wgt = xsWgt*puWgt;
          }

          // make fake ele and fake mu plots separately
          assert(ctrl.faketype=="mu" || ctrl.faketype=="ele");
          assert(abs(w_lep.type)==11 || abs(w_lep.type)==13);
          if( abs(w_lep.type)==11 && ctrl.faketype=="ele") continue;
          if( abs(w_lep.type)==13 && ctrl.faketype=="mu" ) continue;
          
          // selection on the W lepton
          // if(w_lep.isoPF04/w_lep.vec.Pt() > 0.025)   continue;
          if(fs->info->met < 25)                        continue;
          if(!(w_lep.tightCutsApplied))                 continue;
          if(w_lep.vec.Pt() < 25)                       continue;
          if(fabs(w_lep.ip3dSig) > 4)                   continue;
          // if( fs->isdata_)
          //   if( !(w_lep.bdtfail & 2) )                    continue;
          // if(!w_lep.isTight)                            continue; // this bit stores whether the w lepton is trigger matched

          // require an extra OF SS lepton
          bool has_ssof=false;
          vector<SimpleLepton> ssof_leps;
          vector<TString> types;
          for(unsigned iextra=iw+1; iextra<npass; iextra++) {
            SimpleLepton ex_lep = (*fs->passingL)[iextra];
            if(abs(ex_lep.type) == abs(w_lep.type))     continue;
            if(ex_lep.charge != w_lep.charge)           continue;
            has_ssof = true;
            ssof_leps.push_back(ex_lep);
            types.push_back("pass");
          }
          for(unsigned iextra=0; iextra<nfail; iextra++) {
            SimpleLepton ex_lep = (*fs->failingL)[iextra];
            if(abs(ex_lep.type) == abs(w_lep.type))     continue;
            if(ex_lep.charge != w_lep.charge)           continue;
            has_ssof = true;
            ssof_leps.push_back(ex_lep);
            types.push_back("fail");
          }
          if(!has_ssof) continue;

          // fill hists with W + (fake or pass)
          double mt = get_mt(w_lep.vec.Phi(),fs->info->metphi,w_lep.vec.Pt(),fs->info->met);
          (*(cs->hists)["both"])["npv"]->Fill(        fs->info->npv,                    wgt);
          (*(cs->hists)["both"])["w_pt"]->Fill(       w_lep.vec.Pt(),                   wgt);
          (*(cs->hists)["both"])["w_eta"]->Fill(      w_lep.vec.Eta(),                  wgt);
          (*(cs->hists)["both"])["w_isoPF04"]->Fill(  w_lep.isoPF04/w_lep.vec.Pt(),     wgt);
          (*(cs->hists)["both"])["met"]->Fill(        fs->info->met,                    wgt);
          (*(cs->hists)["both"])["mt"]->Fill(         mt,                               wgt);
          (*(cs->hists)["both"])["npass"]->Fill(      npass,                            wgt);
          (*(cs->hists)["both"])["nfail"]->Fill(      nfail,                            wgt);
          (*(cs->hists)["both"])["run"]->Fill(        fs->info->run,                    wgt);
        
          for(unsigned iextra=0; iextra<ssof_leps.size(); iextra++) {
            SimpleLepton ex_lep = ssof_leps[iextra];
            double mass = (w_lep.vec + ex_lep.vec).M();

            if(fabs(ex_lep.ip3dSig) > 4)                  continue;
            if(types[iextra] == "fail") {
              
              // fill fake rate prediction

              double fake_wgt    = get_fake_weight("",ex_lep,*fr);
              double fake_wgt_lo = get_fake_weight("lo",ex_lep,*fr);
              double fake_wgt_hi = get_fake_weight("hi",ex_lep,*fr);

              double centr=fake_wgt*wgt,lo=fake_wgt_lo*wgt,hi=fake_wgt_hi*wgt;
              fillHist( cs, "pred", "w_pt",             w_lep.vec.Pt(),                 centr, lo, hi);
              fillHist( cs, "pred", "l2_pt",            ex_lep.vec.Pt(),                centr, lo, hi);
              fillHist( cs, "pred", "w_eta",            w_lep.vec.Eta(),                centr, lo, hi);
              fillHist( cs, "pred", "l2_eta",           ex_lep.vec.Eta(),               centr, lo, hi);
              fillHist( cs, "pred", "w_isoPF04",        w_lep.isoPF04/w_lep.vec.Pt(),   centr, lo, hi);
              fillHist( cs, "pred", "l2_isoPF04",       ex_lep.isoPF04/ex_lep.vec.Pt(), centr, lo, hi);
              fillHist( cs, "pred", "w_ip3ds",          w_lep.ip3dSig,                  centr, lo, hi);
              fillHist( cs, "pred", "l2_ip3ds",         ex_lep.ip3dSig,                 centr, lo, hi);
              fillHist( cs, "pred", "mass",             mass,                           centr, lo, hi);
              fillHist( cs, "pred", "met",              fs->info->met,                  centr, lo, hi);
              fillHist( cs, "pred", "mt",               mt,                             centr, lo, hi);
              fillHist( cs, "pred", "run",              fs->info->run,                  centr, lo, hi);

              (*(cs->hists)["pred"])["fr"]->Fill(          fake_wgt                             );

            } else if(types[iextra] == "pass") {

              // fill the observation
            
              fillHist( cs, "obs", "w_pt",             w_lep.vec.Pt(),                  wgt, 0, 0);
              fillHist( cs, "obs", "l2_pt",            ex_lep.vec.Pt(),                 wgt, 0, 0);
              fillHist( cs, "obs", "w_eta",            w_lep.vec.Eta(),                 wgt, 0, 0);
              fillHist( cs, "obs", "l2_eta",           ex_lep.vec.Eta(),                wgt, 0, 0);
              fillHist( cs, "obs", "w_isoPF04",        w_lep.isoPF04/w_lep.vec.Pt(),    wgt, 0, 0);
              fillHist( cs, "obs", "l2_isoPF04",       ex_lep.isoPF04/ex_lep.vec.Pt(),  wgt, 0, 0);
              fillHist( cs, "obs", "w_ip3ds",          w_lep.ip3dSig,                   wgt, 0, 0);
              fillHist( cs, "obs", "l2_ip3ds",         ex_lep.ip3dSig,                  wgt, 0, 0);
              fillHist( cs, "obs", "mass",             mass,                            wgt, 0, 0);
              fillHist( cs, "obs", "met",              fs->info->met,                   wgt, 0, 0);
              fillHist( cs, "obs", "mt",               mt,                              wgt, 0, 0);
              fillHist( cs, "obs", "run",              fs->info->run,                   wgt, 0, 0);
            } else { cout << "ERROR: types: " << types[iextra] << endl; assert(0); }
          }
        }
      }
      cout << "\t\tWARNING: not checking for duplicate events" << endl; //skipped " << nDuplSkipped << " duplicate events" << endl;
    }
  }
  cout << "run range: " << setw(12) << minRun << setw(12) << maxRun << endl;
  
  assert(samplev.size()==3);
  CSample *cs = samplev[0];
  CSample *cs_mc = samplev[1];
  CSample *cs_zmc = samplev[2];

  double total_mc_scale=-1;
  RooFitResult *fitres = fitMt(ctrl,(*(cs_zmc->hists)["both"])["mt"],(*(cs_mc->hists)["both"])["mt"],(*(cs->hists)["both"])["mt"],total_mc_scale,plotLabel);

  cout << "\n\n\ntotal mc scale: " << total_mc_scale << endl << endl << endl;

  // RooArgList float_pars = fitres->floatParsFinal();
  // RooRealVar *nWandZ = (RooRealVar*)float_pars.at(float_pars.index("nWandZ"));
  // cout << nWandZ->getVal() << endl;

  gSystem->mkdir(ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype,true);
  TFile runHistFile(ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/runs.root","recreate");
  map<TString,TH1D*>::iterator it_v;
  for(it_v=(*(cs->hists)["pred"]).begin(); it_v!=(*(cs->hists)["pred"]).end(); it_v++) {
    TString var((*it_v).first);

    // plot passing and failing fakes together
    if(var!="mt") {
      (*(cs_mc->hists)["both"])[var]->Scale(total_mc_scale);
      (*(cs_zmc->hists)["both"])[var]->Scale(total_mc_scale);
    }
    TH1D *hBoth = (*(cs->hists)["both"])[var];
    CPlot cplot_both("both_"+var,"",hBoth->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/plots");
    cplot_both.AddHist1D(hBoth,"tight+fk'ble: "+integral_str(hBoth),"E");
    cplot_both.AddToStack((*(cs_zmc->hists)["both"])[var],"Z+j MC: "+integral_str((*(cs_zmc->hists)["both"])[var]),920);
    cplot_both.AddToStack((*(cs_mc->hists)["both"])[var],"W+j MC: "+integral_str((*(cs_mc->hists)["both"])[var]),426);

    if(var.Contains("eta")) cplot_both.SetYRange(0,1.4*hBoth->GetMaximum());
    if(var.Contains("isoPF")) cplot_both.SetLogy();
    if(var=="met") cplot_both.AddLine(25,0,25,18050,843,kDashed);
    cplot_both.Draw(can,true,"png");
    if(cs->isdata && var=="run") {
      assert(hBoth->GetBinContent(0)==0 && hBoth->GetBinContent(hBoth->GetXaxis()->GetNbins()+1)==0);
      hBoth->Write("runs");
    }

    // plot observation vs. prediction
    (*(cs_mc->hists)["obs"])[var]->Scale(total_mc_scale);
    (*(cs_zmc->hists)["obs"])[var]->Scale(total_mc_scale);
    CPlot cplot(var,"",(*(cs->hists)["pred"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/plots");
    cplot.AddHist1D((*(cs->hists)["obs"])[var],cs->legend+": "+integral_str((*(cs->hists)["obs"])[var]),"E");
    cplot.AddHist1D((*(cs->hists)["pred"])[var],"FR predic: "+integral_str((*(cs->hists)["pred"])[var]),"hist",kRed);
    cplot.AddHist1D((*(cs->hists)["pred_lo"])[var],"stat lo: "+integral_str((*(cs->hists)["pred_lo"])[var]),"hist",kRed,kDashed);
    cplot.AddHist1D((*(cs->hists)["pred_hi"])[var],"stat hi: "+integral_str((*(cs->hists)["pred_hi"])[var]),"hist",kRed,kDashed);

    // stack on the mc
    // cplot.AddToStack((*(cs_zmc->hists)["obs"])[var],"Z+j MC: "+integral_str((*(cs_zmc->hists)["obs"])[var]),920);
    // cplot.AddToStack((*(cs_mc->hists)["obs"])[var],"W+j MC: "+integral_str((*(cs_mc->hists)["obs"])[var]),426);

    if(var.Contains("eta")) cplot.SetYRange(0,1.4*max((*(cs->hists)["obs"])[var]->GetMaximum(),(*(cs->hists)["pred"])[var]->GetMaximum()));
    if(var.Contains("isoPF") || var.Contains("ip3ds")) {
      cplot.SetLogy();
      TH1D *hobs  = (*(cs->hists)["obs"])[var];
      TH1D *hpred = (*(cs->hists)["pred"])[var];
      double isoMin = 0.008*min(hobs->GetBinContent(1), hpred->GetBinContent(1));
      double isoMax = 1.8*max(hobs->GetBinContent(1), hpred->GetBinContent(1));

      cplot.SetYRange(isoMin,isoMax);
      // cplot.SetYRange(0,0.8*min((*(cs->hists)["obs"])[var]->GetMinimum(),(*(cs->hists)["pred"])[var]->GetMinimum()));
    }
    cplot.Draw(can,true,"png");
  }
  runHistFile.Close();
  
  makeHTML(ctrl,plotLabel);
}
//----------------------------------------------------------------------------------------
map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str)
{
  map<TString,map<TString,TH1D*>* > hists;
  map<TString,TH1D*> *h_both = new map<TString,TH1D*>;          hists["both"]    = h_both;
  map<TString,TH1D*> *h_pred = new map<TString,TH1D*>;          hists["pred"]    = h_pred;
  map<TString,TH1D*> *h_pred_lo = new map<TString,TH1D*>;       hists["pred_lo"] = h_pred_lo;
  map<TString,TH1D*> *h_pred_hi = new map<TString,TH1D*>;       hists["pred_hi"] = h_pred_hi;
  map<TString,TH1D*> *h_obs  = new map<TString,TH1D*>;          hists["obs"]     = h_obs;
  map<TString,map<TString,TH1D*>* >::iterator it_h;
  double w_pt_max,l2_pt_max,mass_max;
  if(ctrl.faketype=="mu") {
    w_pt_max = 80;
    l2_pt_max = 35;
    mass_max = 120;
  } else if(ctrl.faketype=="ele"){
    w_pt_max = 90;
    l2_pt_max = 60;
    mass_max = 210;
  } else assert(0);

  for(it_h=hists.begin(); it_h!=hists.end(); it_h++) {
    (*((*it_h).second))["npv"]          = new TH1D(TString("npv") +"_"+(*it_h).first+str,";#bf{n PV};",                 50,-.5,49.5);        (*((*it_h).second))["npv"]->Sumw2();
    (*((*it_h).second))["w_pt"]         = new TH1D(TString("w_pt") +"_"+(*it_h).first+str,";#bf{real p_{T} [GeV]};",    37,15,w_pt_max);     (*((*it_h).second))["w_pt"]->Sumw2();
    (*((*it_h).second))["l2_pt"]        = new TH1D(TString("l2_pt") +"_"+(*it_h).first+str,";#bf{fake p_{T} [GeV]};",   25,0,l2_pt_max);     (*((*it_h).second))["l2_pt"]->Sumw2();
    (*((*it_h).second))["w_eta"]        = new TH1D(TString("w_eta")+"_"+(*it_h).first+str,";#bf{real #eta};",           25,-2.5,2.5);        (*((*it_h).second))["w_eta"]->Sumw2();
    (*((*it_h).second))["l2_eta"]       = new TH1D(TString("l2_eta")+"_"+(*it_h).first+str,";#bf{fake #eta};",          25,-2.5,2.5);        (*((*it_h).second))["l2_eta"]->Sumw2();
    (*((*it_h).second))["w_isoPF04"]    = new TH1D(TString("w_isoPF04")+"_"+(*it_h).first+str,";#bf{real PF Iso};",     25,0,1);             (*((*it_h).second))["w_isoPF04"]->Sumw2();
    (*((*it_h).second))["l2_isoPF04"]   = new TH1D(TString("l2_isoPF04")+"_"+(*it_h).first+str,";#bf{fake PF Iso};",    25,0,3);             (*((*it_h).second))["l2_isoPF04"]->Sumw2();
    (*((*it_h).second))["w_ip3ds"]      = new TH1D(TString("w_ip3ds")+"_"+(*it_h).first+str,";#bf{real ip3dSig};",      25,-10,10);            (*((*it_h).second))["w_ip3ds"]->Sumw2();
    (*((*it_h).second))["l2_ip3ds"]     = new TH1D(TString("l2_ip3ds")+"_"+(*it_h).first+str,";#bf{fake ip3dSig};",     25,-10,10);            (*((*it_h).second))["l2_ip3ds"]->Sumw2();
    (*((*it_h).second))["mass"]         = new TH1D(TString("mass")  +"_"+(*it_h).first+str,";#bf{m_{ll}};",             25,0,mass_max);      (*((*it_h).second))["mass"]->Sumw2();
    (*((*it_h).second))["met"]          = new TH1D(TString("met")  +"_"+(*it_h).first+str,";#bf{MET};",                 25,0,100);           (*((*it_h).second))["met"]->Sumw2();
    (*((*it_h).second))["mt"]           = new TH1D(TString("mt")   +"_"+(*it_h).first+str,";#bf{m_{T}};",               25,2,150);           (*((*it_h).second))["mt"]->Sumw2();
    (*((*it_h).second))["fr"]           = new TH1D(TString("fr")   +"_"+(*it_h).first+str,";#bf{FR};",                  125,-.01,.1);        (*((*it_h).second))["fr"]->Sumw2();
    (*((*it_h).second))["npass"]        = new TH1D(TString("npass")   +"_"+(*it_h).first+str,";#bf{N Pass};",           125,-.01,4);         (*((*it_h).second))["npass"]->Sumw2();
    (*((*it_h).second))["nfail"]        = new TH1D(TString("nfail")   +"_"+(*it_h).first+str,";#bf{N Fail};",           125,-.01,4);         (*((*it_h).second))["nfail"]->Sumw2();
    (*((*it_h).second))["run"]          = new TH1D(TString("run")   +"_"+(*it_h).first+str,";#bf{Run};",                50,160405,196550);   (*((*it_h).second))["run"]->Sumw2();
  }
  
  return hists;
}
//--------------------------------------------------------------------------------------------------
void makeHTML(FOFlags &ctrl, TString plotLabel)
{
  // TString title(ctrl.inputdir);
  // title = title(title.Last('/')+1,title.Length());
  TString title("EMU closure test: Fake "+ctrl.faketype);
  TString htmlfname(ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/plots.html");
  ofstream htmlfile(htmlfname);;

  htmlfile << "<!DOCTYPE html" << endl;
  htmlfile << "    PUBLIC \"-//W3C//DTD HTML 3.2//EN\">" << endl;
  htmlfile << "<html>" << endl;

  htmlfile << "<head><title>"+title+"</title></head>" << endl;
  htmlfile << "<body bgcolor=\"EEEEEE\">" << endl;
  htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">"+title+"</h3>" << endl;

  htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;  
  htmlfile << "<tr>" << endl;
  htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/mtfitW.png\"><img src=\"plots/mtfitW.png\" alt=\"plots/mtfitW.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/both_met.png\"><img src=\"plots/both_met.png\" alt=\"plots/both_met.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "</tr>" << endl;
  htmlfile << "</table>" << endl;

  htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;  
  htmlfile << "<tr>" << endl;
  htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/fr.png\"><img src=\"plots/fr.png\" alt=\"plots/fr.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/both_npv.png\"><img src=\"plots/both_npv.png\" alt=\"plots/both_npv.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/run.png\"><img src=\"plots/run.png\" alt=\"plots/run.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "</tr>" << endl;
  htmlfile << "</table>" << endl;

  htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;  
  htmlfile << "<tr>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/mass.png\"><img src=\"plots/mass.png\" alt=\"plots/mass.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/mt.png\"><img src=\"plots/mt.png\" alt=\"plots/mt.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/met.png\"><img src=\"plots/met.png\" alt=\"plots/met.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "</tr>" << endl;
  htmlfile << "</table>" << endl;
  
  htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;  
  htmlfile << "<tr>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/w_pt.png\"><img src=\"plots/w_pt.png\" alt=\"plots/w_pt.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/w_eta.png\"><img src=\"plots/w_eta.png\" alt=\"plots/w_eta.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/w_isoPF04.png\"><img src=\"plots/w_isoPF04.png\" alt=\"plots/w_isoPF04.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/w_ip3ds.png\"><img src=\"plots/w_ip3ds.png\" alt=\"plots/w_ip3ds.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "</tr>" << endl;
  htmlfile << "</table>" << endl;
  
  htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;  
  htmlfile << "<tr>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/l2_pt.png\"><img src=\"plots/l2_pt.png\" alt=\"plots/l2_pt.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/l2_eta.png\"><img src=\"plots/l2_eta.png\" alt=\"plots/l2_eta.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/l2_isoPF04.png\"><img src=\"plots/l2_isoPF04.png\" alt=\"plots/l2_isoPF04.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/l2_ip3ds.png\"><img src=\"plots/l2_ip3ds.png\" alt=\"plots/l2_ip3ds.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "</tr>" << endl;    
  htmlfile << "</table>" << endl;
  htmlfile << "<hr />" << endl;

  // htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;  
  // htmlfile << "<tr>" << endl;
  // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_w_isoPF04.png\"><img src=\"plots/both_w_isoPF04.png\" alt=\"plots/both_w_isoPF04.png\" width=\"100%\"></a></td>" << endl;
  // htmlfile << "</tr>" << endl;  
  // htmlfile << "</table>" << endl;

  htmlfile << "<hr />" << endl;

  htmlfile << "<hr />" << endl;
    
  htmlfile << "</body>" << endl;
  htmlfile << "</html>" << endl;
  htmlfile.close(); 
}
//----------------------------------------------------------------------------------------
RooFitResult *fitMt(FOFlags &ctrl, TH1D *hmcZ, TH1D *hmcW, TH1D *hdata, double &total_mc_scale, TString plotLabel)
{
  double original_w_integral = integrateHist(hmcW);

  RooRealVar mt("mt","MT [GeV]",hdata->GetXaxis()->GetXmin(),hdata->GetXaxis()->GetXmax());
  RooBinning mtbins(hdata->GetNbinsX(),hdata->GetXaxis()->GetXmin(),hdata->GetXaxis()->GetXmax());
  mt.setBinning(mtbins);

  // data
  RooDataHist obsMt("obsMt","obsMt",RooArgList(mt),hdata);
  // W mc pdf
  hmcW->Smooth(2);
  RooDataHist hW("hW","hW",RooArgList(mt),hmcW);
  RooHistPdf  pW("pW","pW",RooArgSet(mt),hW,3);
  // Z mc pdf
  RooDataHist hZ("hZ","hZ",RooArgList(mt),hmcZ);
  RooHistPdf  pZ("pZ","pZ",RooArgSet(mt),hZ,3);
  // rayleigh pdf
  RooRealVar  sig("sig","sig",400,15,585);
  RooRealVar  a1("a1","a1",-3,-10,10);
  // RooGenericPdf rayl("rayl","rayl","mt*exp(-(mt)*(mt)/(2*sig*sig))",RooArgList(mt,sig));
  RooGenericPdf rayl("rayl","rayl","mt*exp(-(mt)*(mt)/(sig*sig +2*sig*a1*(mt) + a1*a1*(mt)*(mt)))",RooArgList(mt,sig,a1));

  RooRealVar nObs("nObs","nObs",integrateHist(hdata));
  RooRealVar nWandZ("nWandZ","nWandZ",0.98*nObs.getVal(),.95*nObs.getVal(),1.1*nObs.getVal());
  RooRealVar zFrac("zFrac","zFrac",integrateHist(hmcZ) / (integrateHist(hmcZ) + integrateHist(hmcW)));
  RooAddPdf modelWandZ("modelWandZ","modelWandZ",RooArgList(pZ,pW),zFrac);
  RooRealVar nRayl("nRayl","nRayl",0.05*nObs.getVal(),0.0*nObs.getVal(),.05*nObs.getVal());

  RooExtendPdf extWandZ("extWandZ","extWandZ",modelWandZ,nWandZ);
  RooExtendPdf extRayl("extRayl","extRayl",rayl,nRayl);

  RooAddPdf *model = new RooAddPdf("model","model",RooArgList(extWandZ,extRayl));

  RooFitResult *fitres = model->fitTo(obsMt,Extended(kTRUE),Save(kTRUE));
  fitres->Print();

  TH1D *rayl_hist = (TH1D*)rayl.createHistogram("rayl_mt",mt,Binning(mtbins));
  rayl_hist->Scale(nRayl.getVal()/integrateHist(rayl_hist));
  TH1D *totalHist = (TH1D*)model->createHistogram("total_mt",mt,Binning(mtbins));
  totalHist->Scale( (nWandZ.getVal() + nRayl.getVal()) / integrateHist(totalHist) );

  // // why the fuck do I need this?
  // double tot = nWandZ.getVal() + nRayl.getVal();
  // if(use_exp) tot += nExp.getVal();
  // double hack_fac = integrateHist(hdata)/tot;
  // rayl_hist->Scale(hack_fac);
  // exp_hist->Scale(hack_fac);
  // hmcZ->Scale(hack_fac);
  // hmcW->Scale(hack_fac);

  // RooPlot *frame = mt.frame();
  // obsMt.plotOn(frame);
  // model->plotOn(frame);
  // model->plotOn(frame,Components("pW"),LineStyle(kDashed),LineColor(kRed));
  // model->plotOn(frame,Components("pZ"),LineColor(kRed));
  // model->plotOn(frame,Components("rayl"),LineColor(kGreen));
  // frame->Draw();
  // can->SaveAs("~/public_html/foo.png");
  
  hmcW->Scale( (nWandZ.getVal()*(1-zFrac.getVal())) / integrateHist(hmcW) );
  hmcZ->Scale( (nWandZ.getVal()*zFrac.getVal()) / integrateHist(hmcZ) );

  CPlot plotMtFit("mtfitW","","MT [GeV]","Events",ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/plots");
  plotMtFit.AddHist1D(hdata,"tight+fk'ble: " + integral_str(hdata),"E");
  plotMtFit.AddHist1D(totalHist,"total MC: " + integral_str(totalHist),"hist",kRed);
  plotMtFit.AddToStack(rayl_hist,"jet (Rayleigh): " + integral_str(rayl_hist),kBlue);
  plotMtFit.AddToStack(hmcZ,"Z+j MC: " + integral_str(hmcZ),kGray);
  plotMtFit.AddToStack(hmcW,"W+j MC: " + integral_str(hmcW),kCyan-6);
  plotMtFit.Draw(can,kTRUE,"png");

  total_mc_scale = integrateHist(hmcW) / original_w_integral;

  return fitres;
}
//----------------------------------------------------------------------------------------
void fillHist(CSample *cs, TString type, TString var, double val, double wgt, double wgt_lo, double wgt_hi)
{
  (*(cs->hists)[type])[var]->Fill(         val,   wgt);
  if(type=="pred") {
    (*(cs->hists)[type+"_lo"])[var]->Fill(   val,   wgt_lo);
    (*(cs->hists)[type+"_hi"])[var]->Fill(   val,   wgt_hi);
  }
}
Revision:	1.11
Committed:	Mon Jul 16 10:21:13 2012 UTC (12 years, 10 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
Changes since 1.10:	+7 -4 lines
Log Message:	Add denominator def with Dz and D0 but no IP, move some things from Various to PlotHeaders, clean up plotters