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"

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

using namespace std;
using namespace RooFit;

TCanvas *can;

double get_fake_weight(TString uncert, SimpleLepton fake_lep, FR_struct fr);
void makeHTML(FOFlags &ctrl);
map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str="");
TVector3 get_corrected_met(double met, double metphi, TLorentzVector lepton);
TString entries_str(TH1D *hist);
TString integral_str(TH1D *hist);
double integrateHist(TH1D *hist) { return hist->Integral(0,hist->GetNbinsX()+1); }

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

//----------------------------------------------------------------------------------------
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 fr = initFRs(ctrl.mufakedir+"/fr.root",ctrl.elefakedir+"/fr.root");
  FR_struct fr = initFRs(ctrl.mufakedir,ctrl.elefakedir);

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

  vector<CSample*> samplev;

  ifstream instream;
  instream.open(ctrl.config); assert(instream.is_open());
  TString ntupledir;
  string line;
  while(getline(instream,line)) {
    if(line[0]=='#') continue;

    stringstream ss(line);
    if(line[0]=='^') {
      TString dummy;
      if(TString(line).Contains("^ntupdir")) ss >> dummy >> ntupledir;
      // if(TString(line).Contains("^skim")) ss >> dummy >> ntupledir;
      continue; 
    }

    if(line[0]=='$') {
      samplev.push_back(new CSample());
      stringstream ss(line);
      string chr;
      TString sname;
      Int_t color;
      ss >> chr >> sname >> color;
      string legend = line.substr(line.find('@')+1);
      samplev.back()->name  = sname;
      samplev.back()->legend = legend;
      samplev.back()->color = color;
      samplev.back()->hists = init_hists(ctrl,sname);
      continue;
    }
    
    TString dset;
    bool isdata;
    double xsec;
    ss >> dset >> isdata >> xsec;
    assert(ctrl.muSele==ctrl.eleSele);
    TString fname = ntupledir+"/"+ctrl.muSele+"/"+dset+"/merged.root";
    (samplev.back()->fsv).push_back(new filestuff(dset,fname,dset,isdata));
  }
  instream.close();


  for(unsigned ics=0; ics<samplev.size(); ics++) {
    CSample *cs = samplev[ics];
    
    for(unsigned ifs=0; ifs<(cs->fsv).size(); ifs++) {
      
      filestuff *fs = (cs->fsv)[ifs];
      for(unsigned ientry=0; ientry<fs->getentries(); ientry++) {
        fs->getentry(ientry,"info");

        // veto extra leptons
        unsigned npass = fs->passingL->size();
        unsigned nfail = fs->failingL->size();
        if(npass+nfail != 2) continue;

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

          double wgt = fs->isdata_ ? 1 : lumi*xstab.Get(fs->dataset_)/fs->total_entries_;

          // selection on the W lepton
          // if(w_lep.isoPF04/w_lep.vec.Pt() > 0.025)   continue;
          if(!(w_lep.tightCutsApplied))                 continue;
          if((abs(w_lep.type)==11) && (ctrl.faketype=="mu")) {
            if(w_lep.vec.Pt() < 20)                     continue;
          } else if((abs(w_lep.type)==13) && (ctrl.faketype=="ele")) {
            if(w_lep.vec.Pt() < 25)                     continue;
          } else continue;

          // require an extra OF SS lepton
          bool has_ssof=false;
          vector<SimpleLepton> ssof_leps;
          vector<TString> types;
          for(unsigned iextra=0; iextra<fs->passingL->size(); iextra++) {
            SimpleLepton ex_lep = (*fs->passingL)[iextra];

            if(iextra == iw)                              continue; // unnecessary
            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<fs->failingL->size(); 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);
          if(fs->info->met < 25)  continue;
          (*(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);
        
          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(dr(w_lep,ex_lep) < 0.05) continue;
            //????????????????????????????????????????????????????????????????????????????????????????

            if(mass < 1) {
              cout << "low mass: " << mass << endl
                   << "\tW lepton:     " << w_lep.type << setw(12) << w_lep.vec.Eta() << setw(12) << w_lep.vec.Phi()
                   << "\tf'ble lepton: " << ex_lep.type << setw(12) << ex_lep.vec.Eta() << setw(12) << ex_lep.vec.Phi()
                   << endl;
            }
            
            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);

              (*(cs->hists)["pred"])["w_pt"]->Fill(        w_lep.vec.Pt(),                      fake_wgt*wgt);
              (*(cs->hists)["pred"])["l2_pt"]->Fill(       ex_lep.vec.Pt(),                     fake_wgt*wgt);
              (*(cs->hists)["pred_lo"])["l2_pt"]->Fill(    ex_lep.vec.Pt(),                     fake_wgt_lo*wgt);
              (*(cs->hists)["pred_hi"])["l2_pt"]->Fill(    ex_lep.vec.Pt(),                     fake_wgt_hi*wgt);
              (*(cs->hists)["pred"])["w_eta"]->Fill(       w_lep.vec.Eta(),                     fake_wgt*wgt);
              (*(cs->hists)["pred"])["l2_eta"]->Fill(      ex_lep.vec.Eta(),                    fake_wgt*wgt);
              (*(cs->hists)["pred_lo"])["l2_eta"]->Fill(   ex_lep.vec.Eta(),                    fake_wgt_lo*wgt);
              (*(cs->hists)["pred_hi"])["l2_eta"]->Fill(   ex_lep.vec.Eta(),                    fake_wgt_hi*wgt);
              (*(cs->hists)["pred"])["w_isoPF04"]->Fill(   w_lep.isoPF04/w_lep.vec.Pt(),        fake_wgt*wgt);
              (*(cs->hists)["pred"])["l2_isoPF04"]->Fill(  ex_lep.isoPF04/ex_lep.vec.Pt(),      fake_wgt*wgt);
              (*(cs->hists)["pred"])["mass"]->Fill(        mass,                                fake_wgt*wgt);
              (*(cs->hists)["pred_lo"])["mass"]->Fill(     mass,                                fake_wgt_lo*wgt);
              (*(cs->hists)["pred_hi"])["mass"]->Fill(     mass,                                fake_wgt_hi*wgt);
              (*(cs->hists)["pred"])["met"]->Fill(         fs->info->met,                       fake_wgt*wgt);
              (*(cs->hists)["pred"])["mt"]->Fill(          mt,                                  fake_wgt*wgt);
              (*(cs->hists)["pred"])["fr"]->Fill(          fake_wgt                             );
            } else if(types[iextra] == "pass") {

              // fill the observation
            
              (*(cs->hists)["obs"])["w_pt"]->Fill(        w_lep.vec.Pt(),                 wgt);
              (*(cs->hists)["obs"])["l2_pt"]->Fill(       ex_lep.vec.Pt(),                wgt);
              (*(cs->hists)["obs"])["w_eta"]->Fill(       w_lep.vec.Eta(),                wgt);
              (*(cs->hists)["obs"])["l2_eta"]->Fill(      ex_lep.vec.Eta(),               wgt);
              (*(cs->hists)["obs"])["w_isoPF04"]->Fill(   w_lep.isoPF04/w_lep.vec.Pt(),   wgt);
              (*(cs->hists)["obs"])["l2_isoPF04"]->Fill(  ex_lep.isoPF04/ex_lep.vec.Pt(), wgt);
              (*(cs->hists)["obs"])["mass"]->Fill(        mass,                           wgt);
              (*(cs->hists)["obs"])["met"]->Fill(         fs->info->met,                  wgt);
              (*(cs->hists)["obs"])["mt"]->Fill(          mt,                             wgt);
            } else { cout << "ERROR: types: " << types[iextra] << endl; assert(0); }
          }
        }
      }
    }
  }
  
  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);

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

  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
    if(var!="mt") {
      (*(cs_mc->hists)["both"])[var]->Scale(total_mc_scale);
      (*(cs_zmc->hists)["both"])[var]->Scale(total_mc_scale);
    }
    CPlot cplot_both("both_"+var,"",(*(cs->hists)["both"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/plots");
    cplot_both.AddHist1D((*(cs->hists)["both"])[var],"tight+fk'ble: "+integral_str((*(cs->hists)["both"])[var]),"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*(*(cs->hists)["both"])[var]->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");

    // 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+"/plots");
    cplot.AddHist1D((*(cs->hists)["obs"])[var],"Obs: "+integral_str((*(cs->hists)["obs"])[var]),"E");
    cplot.AddHist1D((*(cs->hists)["pred"])[var],"FR predic: "+integral_str((*(cs->hists)["pred"])[var]),"hist",kRed);
    if(var=="mass" || var=="l2_pt" || var=="l2_eta") {
      cplot.AddHist1D((*(cs->hists)["pred_lo"])[var],"FR lo: "+integral_str((*(cs->hists)["pred_lo"])[var]),"hist",kRed,kDashed);
      cplot.AddHist1D((*(cs->hists)["pred_hi"])[var],"FR hi: "+integral_str((*(cs->hists)["pred_hi"])[var]),"hist",kRed,kDashed);
    }
    
    // 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")) cplot.SetLogy();
    cplot.Draw(can,true,"png");
  }
  
  makeHTML(ctrl);
}
//----------------------------------------------------------------------------------------
double get_fake_weight(TString uncert, SimpleLepton fake_lep, FR_struct fr)
{
  double rate=0;
  assert(uncert=="" || uncert=="lo" || uncert=="hi");
  if(abs(fake_lep.type) == 13) { // fake muons:
    double x = fake_lep.vec.Eta();
    if(fr.absEtaMu) x = fabs(x);
    double y = fake_lep.vec.Pt();
    rate = fr.mufr.getEff(x,y,true);
    if(uncert == "lo") rate -= fr.mufr.getErrLow(x,y,true);
    if(uncert == "hi") rate += fr.mufr.getErrHigh(x,y,true);
  } else if(abs(fake_lep.type) == 11) {
    double x = fake_lep.vec.Eta();
    if(fr.absEtaEle) x = fabs(x);
    double y = fake_lep.vec.Pt();
    rate = fr.elefr.getEff(x,y,true);
    if(uncert == "lo") rate -= fr.elefr.getErrLow(x,y,true);
    if(uncert == "hi") rate += fr.elefr.getErrHigh(x,y,true);
  } else assert(0);
  
  return rate/(1-rate);
}
//----------------------------------------------------------------------------------------
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))["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))["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();
  }
  
  return hists;
}
//--------------------------------------------------------------------------------------------------
void makeHTML(FOFlags &ctrl)
{
  // TString title(ctrl.inputdir);
  // title = title(title.Last('/')+1,title.Length());
  TString title("EMU closure test: Fake "+ctrl.faketype);
  ofstream htmlfile;
  char htmlfname[100];
  sprintf(htmlfname,"%s/"+ctrl.faketype+".html",ctrl.outdir.Data());
  htmlfile.open(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></a></td>" << endl;
  htmlfile << "<td width=\"33.3%\"><a></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=\"33.3%\"><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=\"33.3%\"><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=\"33.3%\"><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 << "</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/l2_pt.png\"><img src=\"plots/l2_pt.png\" alt=\"plots/l2_pt.png\" width=\"100%\"></a></td>" << endl;
  htmlfile << "<td width=\"33.3%\"><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=\"33.3%\"><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 << "</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(); 
}
//----------------------------------------------------------------------------------------
TVector3 get_corrected_met(double met, double metphi, TLorentzVector lepton)
{
  TLorentzVector corr_met;
  corr_met.SetPtEtaPhiM(met,0,metphi,0);
  TVector3 w_lep_trans(lepton.Px(),lepton.Py(),0);
  return corr_met.Vect() + w_lep_trans;
}
//----------------------------------------------------------------------------------------
TString entries_str(TH1D *hist)
{
  int entries = hist->GetEntries();
  stringstream ss;
  ss << entries;
  TString str;
  ss >> str;
  return str;
}
//----------------------------------------------------------------------------------------
TString integral_str(TH1D *hist)
{
  stringstream ss;
  ss << fixed << setprecision(1) << integrateHist(hist);
  TString str;
  ss >> str;
  return str;
}
//----------------------------------------------------------------------------------------
RooFitResult *fitMt(FOFlags &ctrl, TH1D *hmcZ, TH1D *hmcW, TH1D *hdata, double &total_mc_scale)
{
  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+"/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;
}
Revision:	1.5
Committed:	Tue Jun 19 17:11:28 2012 UTC (12 years, 11 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+41 -20 lines
Log Message:	* empty log message *