Plotting/Modules/WZStudy.C

#include <iostream>
#include <vector>
#include <sys/stat.h>

#include <TCut.h>
#include <TROOT.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TColor.h>
#include <TPaveText.h>
#include <TRandom.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TSQLResult.h>
#include <TProfile.h>
#include <TLegendEntry.h>

using namespace std;

namespace TriLeptons {
    float Rsfof_CorrectionFactor=-1;
    float Rsfof_CorrectionFactor_Error=-1;
}

string IdentifierForUnderlinedEntry(string Entry) {
  if(Contains(Entry,"t_bar_t__")) return "ttbar";
  if(Contains(Entry,"Single_top__")) return "SingleTop";
  if(Contains(Entry,"W_Jets__")) return "W+Jets";
  if(Contains(Entry,"WZ__")) return "WZ";
  if(Contains(Entry,"Z_Jets")) return "Z+Jets";
  if(Contains(Entry,"DATA")) return "Data";
  if(Contains(Entry,"2l2q_")) return "ZZ->2l2q";
  if(Contains(Entry,"4l_")) return "ZZ->4l";
  return Entry;
}

void WZ_kin_plot(string variable, int nbins, float min, float high, string xlabel, TCut cut, string saveas) {
  TCanvas *can = new TCanvas("can","can");
  TPad *pad = new TPad("pad","pad",0,0,1,1);
  pad->cd();

  THStack stack = allsamples.DrawStack("stack",variable,nbins,min,high,xlabel, "events", cut,mc, luminosity);
  TH1F *hdata = SingleLeptonData.Draw("hdata",variable,nbins,min,high,xlabel, "events", cut,data, luminosity);

  float StackMax=CollapseStack(stack)->GetMaximum();
  if(StackMax>hdata->GetMaximum()) hdata->SetMaximum(1.2*StackMax);
  hdata->Draw("e1");
  stack.Draw("histo,same");
  hdata->Draw("e1,same");

  TLegend *kinleg = allsamples.allbglegend();
  kinleg->Draw();
  DrawPrelim();
  
  float total=0;
  float wz=0;
//  cout << "Numbers with variable " << variable << " and cut " << cut << ": " << endl;
  TIter nextSF = TIter(stack.GetHists());
  TH1F *h;
  stringstream spurity;
  int bracecounter=0;


  while ( h = (TH1F*)nextSF() ) {
//    cout << IdentifierForUnderlinedEntry(h->GetName()) << " : " << h->Integral() << endl;
    if(IdentifierForUnderlinedEntry(h->GetName())=="WZ") {
      wz+=h->Integral();
      total+=h->Integral();
    } else {
      if(IdentifierForUnderlinedEntry(h->GetName())!="Data") total+=h->Integral();
    }
    bracecounter++;
    spurity << "#splitline{" << IdentifierForUnderlinedEntry(h->GetName()) << " : " << h->Integral() << "}{";
  }
  
  spurity << "#splitline{Purity: " << wz/total << "}{N(events): " << total << "}";
  for(int i=0;i<bracecounter;i++) spurity << "}";
  
  TPaveText *PurityPaveText = new TPaveText(0.60, 0.40,0.89, 0.6,"blNDC");
  PurityPaveText->SetFillStyle(4000);
  PurityPaveText->SetBorderSize(0);
  PurityPaveText->SetFillColor(kWhite);
  PurityPaveText->SetTextFont(42);
  PurityPaveText->SetTextSize(0.021);
  PurityPaveText->AddText(spurity.str().c_str());
  //PurityPaveText->Draw();

  
//  TText *purity = new TText(0.8,0.6,spurity.str().c_str());
//  purity->Draw();
//  cout << "DATA : " << hdata->Integral() << endl;
  
  Save_With_Ratio(hdata,stack,pad->cd(),"WZ_study/"+saveas);
  
  float sumData=0;
  float sumMC=0;
  TH1F *cstack = CollapseStack(stack);
  
  for(int i=1;i<=hdata->GetNbinsX();i++) {
    if(hdata->GetBinCenter(i)>=20 && hdata->GetBinCenter(i)<=70) {
      sumData+=hdata->GetBinContent(i);
      sumMC+=cstack->GetBinContent(i);
    }
  }
  delete cstack;
  cout << " [20,70] for " << saveas << " : " << sumData << " (data) vs " << sumMC << " (mc) " << endl;
  
  CleanLegends();
  delete hdata;
  delete can;
}

void SF_vs_OF_plots(string variable, int nbins, float min, float high, string xlabel, TCut cut, TCut sf, TCut of, int mcordata, string saveas) {
    TCanvas *can = new TCanvas("can","can");
    TPad *pad = new TPad("pad","pad",0,0,1,1);
    pad->cd();
    
    TLegend *kinleg = make_legend();
    if(mcordata==data) kinleg->SetHeader("Data");
    else kinleg->SetHeader("MC");
    
    TH1F *sfdata;
    TH1F *ofdata;
    
    if(mcordata==data) {
      ofdata = SingleLeptonData.Draw("ofdata",variable,nbins,min,high,xlabel, "events", cut&&of,mcordata, luminosity);
      sfdata = SingleLeptonData.Draw("sfdata",variable,nbins,min,high,xlabel, "events", cut&&sf,mcordata, luminosity);
    } else {
      ofdata = allsamples.Draw("ofdata",variable,nbins,min,high,xlabel, "events", cut&&of,mcordata, luminosity);
      sfdata = allsamples.Draw("sfdata",variable,nbins,min,high,xlabel, "events", cut&&sf,mcordata, luminosity);
    }
    
    float OFmax=ofdata->GetMaximum();
    if(OFmax>sfdata->GetMaximum()) sfdata->SetMaximum(1.2*OFmax);
    ofdata->SetLineColor(kRed);
    sfdata->Draw("e1");
    ofdata->Draw("histo,same");
    sfdata->Draw("e1,same");
    
    kinleg->AddEntry(sfdata,"Same flavor","p");
    kinleg->AddEntry(ofdata,"Opposite flavor","l");
    kinleg->Draw();
    if(mcordata==data) DrawPrelim();
    else DrawMCPrelim();
    
    Save_With_Ratio(sfdata,ofdata,pad->cd(),"WZ_study/SF_vs_OF_plots/"+saveas);
    delete sfdata;
    delete ofdata;
    delete can;
}

void write_report(TCut BasicCut, TCut cut, ofstream &report) {
  THStack stack = allsamples.DrawStack("stack","tri_mll",10,0,1000,"m_{ll}", "events", BasicCut&&cut,mc, luminosity);
  float total=0;
  float wz=0;
  TIter nextSF = TIter(stack.GetHists());
  TH1F *h;
  while ( h = (TH1F*)nextSF() ) {
    cout << IdentifierForUnderlinedEntry(h->GetName()) << " : " << h->Integral() << endl;
    if(IdentifierForUnderlinedEntry(h->GetName())=="WZ") {
      wz+=h->Integral();
      total+=h->Integral();
    } else {
      if(IdentifierForUnderlinedEntry(h->GetName())!="Data") total+=h->Integral();
    }
    report << h->Integral() << ";";
  }
  report << wz/total << ";\n";
  report << flush;

}

vector<string> producecuts(string variablename, float cut1, float cut2, float cut3=-1, float cut4=-1, float cut5=-1, float cut6=-1, float cut7=-1, float cut8=-1, float cut9=-1, float cut10=-1, float cut11=-1 ) {
  float allcuts[12] = {cut1,cut2,cut3,cut4,cut5,cut6,cut7,cut8,cut9,cut10,cut11};
  vector<string> cuts;
  for(int i = 0; i < 11; ++i ) {
    if(allcuts[i]<0) continue;
    stringstream nowcut;
    nowcut << variablename << allcuts[i];
    cuts.push_back(nowcut.str());
  }
  return cuts;
}
    

void OptimizeSelection(TCut BasicCut) {
  ofstream report;
  report.open ("report3.txt");
  
  vector<string> mt   = producecuts("tri_mT>",0,10,15,20,25,30,35,40,45,50);
  vector<string> et   = producecuts("met[4]>",0,10,20,30,40,50);
  vector<string> mll  = producecuts("abs(tri_mll-91.2)<",10000,5,10,15,20);
  vector<string> mlll = producecuts("abs(tri_mlll-91)>",0,5,10,15,20,25,30,40,50);
  
  for(int imt=0;imt<mt.size();imt++) {
    for(int iet=0;iet<et.size();iet++) {
      for(int imll=0;imll<mll.size();imll++) {
        for(int imlll=0;imlll<mlll.size();imlll++) {
          stringstream cut;
          cut << mt[imt] << " && " << et[iet] << " && " << mll[imll] << " && " << mlll[imlll];
          report << mt[imt] << ";" << et[iet] << ";" << mll[imll] << ";" << mlll[imlll] <<";";
//        cout << "Cut is : " << cut.str() << endl;
          write_report(BasicCut, TCut(cut.str().c_str()),report);
        }
      }
    }
  }
  report.close();
}

void SF_OF_comparison(TCut BasicCut, TCut Selection) {
  int nbins=38;
  float mll_low=10;
  float mll_hi=200;
  
  // all MC samples
  TH1F *sfmc   = allsamples.Draw("sfmc"  ,"tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepSF,mc  , luminosity);
  TH1F *ofmc   = allsamples.Draw("ofmc"  ,"tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepOF,mc  , luminosity);
  
  TCanvas *can = new TCanvas("can","can");
  TPad *pad = new TPad("pad","pad",0,0,1,1);
  pad->cd();
  if(ofmc->GetMaximum()>sfmc->GetMaximum()) sfmc->SetMaximum(1.3*ofmc->GetMaximum());
  ofmc->SetLineColor(kRed);
  TLegend *leg = make_legend();
  leg->AddEntry(sfmc,"same flavor","P");
  leg->AddEntry(ofmc,"opposite flavor","L");
  leg->SetHeader("MC:");
  sfmc->Draw("e1");
  ofmc->Draw("histo,same");
  sfmc->Draw("e1,same");
  leg->Draw();
  DrawMCPrelim();

  Save_With_Ratio(sfmc,ofmc,pad->cd(),"WZ_study/OFSFcomparison_MC");
  
  // WZ only
  TH1F *sfwz   = allsamples.Draw("sfwz"  ,"tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepSF,mc  , luminosity,allsamples.FindSample("WZJetsTo3"));
  TH1F *ofwz   = allsamples.Draw("ofwz"  ,"tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepOF,mc  , luminosity,allsamples.FindSample("WZJetsTo3"));
  
  TPad *pad2 = new TPad("pad","pad",0,0,1,1);
  pad2->cd();
  if(ofwz->GetMaximum()>sfwz->GetMaximum()) sfwz->SetMaximum(1.3*ofwz->GetMaximum());
  ofwz->SetLineColor(kRed);
  sfwz->Draw("e1");
  ofwz->Draw("histo,same");
  sfwz->Draw("e1,same");
  leg->SetHeader("WZ MC:");
  leg->Draw();
  DrawMCPrelim();
  
  Save_With_Ratio(sfwz,ofwz,pad2->cd(),"WZ_study/OFSFcomparison_WZ");
  
  // data
  TH1F *sfdata = SingleLeptonData.Draw("sfdata","tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepSF,data, luminosity);
  TH1F *ofdata = SingleLeptonData.Draw("ofdata","tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepOF,data, luminosity);
  
  TPad *pad3 = new TPad("pad","pad",0,0,1,1);
  pad3->cd();
  if(ofdata->GetMaximum()>sfdata->GetMaximum()) sfdata->SetMaximum(1.3*ofdata->GetMaximum());
  ofdata->SetLineColor(kRed);
  sfdata->Draw("e1");
  ofdata->Draw("histo,same");
  sfdata->Draw("e1,same");
  leg->SetHeader("Data:");
  leg->Draw();
  DrawPrelim();
  
  Save_With_Ratio(sfdata,ofdata,pad3->cd(),"WZ_study/OFSFcomparison_data");

  delete can;
  
  dout << "SF events: " << sfdata->Integral() << "   (mc : " << sfmc->Integral() << " , wz : " << sfwz->Integral() << " )   WZ purity: " << 100*sfwz->Integral()/sfmc->Integral() << " %" << endl;
  dout << "OF events: " << ofdata->Integral() << "   (mc : " << ofmc->Integral() << " , wz : " << ofwz->Integral() << " )   WZ purity: " << 100*ofwz->Integral()/ofmc->Integral() << " %" << endl;
  
  dout << "RATIO : SF / OF = " << sfdata->Integral() / ofdata->Integral() << " +/- " << (sfdata->Integral()/ofdata->Integral()) * TMath::Sqrt(1/sfdata->Integral()+1/ofdata->Integral()) << " (stat) " << endl;
    
  TriLeptons::Rsfof_CorrectionFactor=ofwz->Integral() / sfwz->Integral();
  TriLeptons::Rsfof_CorrectionFactor_Error=(ofwz->Integral()/sfwz->Integral()) * TMath::Sqrt(1/sfwz->GetEntries()+1/ofwz->GetEntries());
  delete sfdata;
  delete ofdata;
  delete sfwz;
  delete ofwz;
  
}

void DrawMatchGenDistribution(string variable, int nbins, float low, float high, bool uselog, string xlabel, string SaveAs, TCut basecut) {
  TLegend *leg = make_legend();
  leg->SetY1(0.7);
  TCanvas *can = new TCanvas("can","can");
  TPad *pad = new TPad("pad","pad",0,0,1,1);
  pad->cd();
  pad->SetLogy(uselog);
  
  TCut IsMatch("(gentri_GoodWMatch&&gentri_GoodZMatch)");
  TCut NoMatch("!(gentri_GoodWMatch&&gentri_GoodZMatch)");
  
  TH1F *goodMatch = allsamples.Draw("goodMatch",variable,nbins,low,high,xlabel,"events",basecut&&IsMatch,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  TH1F *noMatch =   allsamples.Draw("noMatch",  variable,nbins,low,high,xlabel,"events",basecut&&NoMatch,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  
  noMatch->SetLineColor(TColor::GetColor("#FF0000"));//an evil red
  goodMatch->SetLineColor(TColor::GetColor("#04B404")); // a nice green
  
//   noMatch->Scale(1.0/noMatch->Integral());
//   goodMatch->Scale(1.0/goodMatch->Integral());
  
  if(noMatch->GetMaximum()>goodMatch->GetMaximum()) goodMatch->SetMaximum(1.3*noMatch->GetMaximum());
  
  goodMatch->Draw("histo");
  noMatch->Draw("histo,same");
  
  leg->AddEntry(goodMatch,"correct match","l");
  leg->AddEntry(noMatch,"incorrect match","l");
  leg->Draw();
  
  Save_With_Ratio(goodMatch,noMatch,pad->cd(),"WZ_study/GenMatchComparison/"+SaveAs);
  
  delete goodMatch;
  delete noMatch;
  delete can;
}

void DrawMatchRecoDistribution(string variable, int nbins, float low, float high, bool uselog, string xlabel, string SaveAs, TCut basecut) {
  TLegend *leg = make_legend();
  leg->SetY1(0.7);
  TCanvas *can = new TCanvas("can","can");
  TPad *pad = new TPad("pad","pad",0,0,1,1);
  pad->cd();
  pad->SetLogy(uselog);
  
  TCut IsMatch("(tri_GoodWMatch&&tri_GoodZMatch)");
  TCut NoMatch("!(tri_GoodWMatch&&tri_GoodZMatch)");
  
  TH1F *goodMatch = allsamples.Draw("goodMatch",variable,nbins,low,high,xlabel,"events",basecut&&IsMatch,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  TH1F *noMatch =   allsamples.Draw("noMatch",  variable,nbins,low,high,xlabel,"events",basecut&&NoMatch,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  
  noMatch->SetLineColor(TColor::GetColor("#FF0000"));//an evil red
  goodMatch->SetLineColor(TColor::GetColor("#04B404")); // a nice green
  
//   noMatch->Scale(1.0/noMatch->Integral());
//   goodMatch->Scale(1.0/goodMatch->Integral());
  
  if(noMatch->GetMaximum()>goodMatch->GetMaximum()) goodMatch->SetMaximum(1.3*noMatch->GetMaximum());
  
  goodMatch->Draw("histo");
  noMatch->Draw("histo,same");
  
  leg->AddEntry(goodMatch,"correct match","l");
  leg->AddEntry(noMatch,"incorrect match","l");
  leg->Draw();
  DrawPrelim();
  
  Save_With_Ratio(goodMatch,noMatch,pad->cd(),"WZ_study/RecoMatchComparison/"+SaveAs);
  
  delete goodMatch;
  delete noMatch;
  delete can;
}

void GenSFOFCalculator(string title, TCut Cut) {
  TH1F *histoSF = allsamples.Draw("histoSF","gentri_mll",1,0,1000,"processing","events",Cut&&TCut("abs(gentri_id2)==abs(gentri_id3)"),mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  TH1F *histoOF = allsamples.Draw("histoOF","gentri_mll",1,0,1000,"processing","events",Cut&&TCut("abs(gentri_id2)!=abs(gentri_id3)"),mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  cout << title << endl;
  cout << "     " << (const char*) Cut << endl;
  cout << "     SF:    \t" << histoSF->Integral() << " +/- " << sqrt(histoSF->GetEntries())*(histoSF->Integral()/histoSF->GetEntries()) << endl;
  cout << "     OF:    \t" << histoOF->Integral() << " +/- " << sqrt(histoOF->Integral())*(histoOF->Integral()/histoOF->GetEntries()) << endl;
  cout << "     SF/OF: \t" << histoSF->Integral() / histoOF->Integral() << " +/- " << histoSF->GetEntries() / histoOF->GetEntries() * sqrt(1.0/histoSF->GetEntries() + 1.0 / histoOF->GetEntries()) << endl;
  
  delete histoSF;
  delete histoOF;
}
  
void reco_SF_OF_comparison(TCut Base, TCut Selection) {
  
  bool dolog=true;
  bool nolog=false;
  
  DrawMatchRecoDistribution("tri_mlll",60,0,300,dolog,"m_{lll}","TrileptonMass",Base);
  DrawMatchRecoDistribution("tri_mll",60,0,300,dolog,"m_{l_{1}l_{2}}","LeadingDileptonMass",Base);
  DrawMatchRecoDistribution("tri_mT",50,0,100,nolog,"m_{T}","MT",Base);
  DrawMatchRecoDistribution("tri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","SubleadingDileptonMass",Base);
  DrawMatchRecoDistribution("tri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","SubleadingDileptonMass_SF",Base&&TCut("abs(tri_id2)==abs(tri_id3)"));
  DrawMatchRecoDistribution("tri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","SubleadingDileptonMass_OF",Base&&TCut("abs(tri_id2)!=abs(tri_id3)"));
  DrawMatchRecoDistribution("genMET",60,0,300,dolog,"gen MET","MET",Base);
  DrawMatchRecoDistribution("abs(tri_id1)",50,-0.5,49.5,nolog,"id(l_{1})","genID1",Base);
  DrawMatchRecoDistribution("abs(tri_id2)",50,-0.5,49.5,nolog,"id(l_{2})","genID2",Base);
  DrawMatchRecoDistribution("abs(tri_id3)",50,-0.5,49.5,nolog,"id(l_{3})","genID3",Base);
  DrawMatchRecoDistribution("abs(tri_genMID1)",50,-0.5,49.5,nolog,"mother id(l_{1})","genMID1",Base);
  DrawMatchRecoDistribution("abs(tri_genMID2)",50,-0.5,49.5,nolog,"mother id(l_{2})","genMID2",Base);
  DrawMatchRecoDistribution("abs(tri_genMID3)",50,-0.5,49.5,nolog,"mother id(l_{3})","genMID3",Base);
  
  DrawMatchRecoDistribution("tri_mlll",60,0,300,dolog,"m_{lll}","FullSelection/TrileptonMass",Base&&Selection);
  DrawMatchRecoDistribution("tri_mll",60,0,300,dolog,"m_{l_{1}l_{2}}","FullSelection/LeadingDileptonMass",Base&&Selection);
  DrawMatchRecoDistribution("tri_mT",50,0,100,nolog,"m_{T}","FullSelection/MT",Base&&Selection);
  DrawMatchRecoDistribution("tri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","FullSelection/SubleadingDileptonMass",Base&&Selection);
  DrawMatchRecoDistribution("tri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","FullSelection/SubleadingDileptonMass_SF",Base&&TCut("abs(tri_id2)==abs(tri_id3)")&&Selection);
  DrawMatchRecoDistribution("tri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","FullSelection/SubleadingDileptonMass_OF",Base&&TCut("abs(tri_id2)!=abs(tri_id3)")&&Selection);
  DrawMatchRecoDistribution("genMET",60,0,300,dolog,"gen MET","FullSelection/MET",Base&&Selection);
  DrawMatchRecoDistribution("abs(tri_id1)",50,-0.5,49.5,nolog,"id(l_{1})","FullSelection/genID1",Base&&Selection);
  DrawMatchRecoDistribution("abs(tri_id2)",50,-0.5,49.5,nolog,"id(l_{2})","FullSelection/genID2",Base&&Selection);
  DrawMatchRecoDistribution("abs(tri_id3)",50,-0.5,49.5,nolog,"id(l_{3})","FullSelection/genID3",Base&&Selection);
  DrawMatchRecoDistribution("abs(tri_genMID1)",50,-0.5,49.5,nolog,"mother id(l_{1})","FullSelection/genMID1",Base&&Selection);
  DrawMatchRecoDistribution("abs(tri_genMID2)",50,-0.5,49.5,nolog,"mother id(l_{2})","FullSelection/genMID2",Base&&Selection);
  DrawMatchRecoDistribution("abs(tri_genMID3)",50,-0.5,49.5,nolog,"mother id(l_{3})","FullSelection/genMID3",Base&&Selection);
}
  
void generator_SF_OF_comparison() {
  /*we want to analyze the cuts at generator level (and their impact), plus understand how bad matching can affect the selection. 
    we want to know the distributions for correctly and incorrectly matched W/Z for the following variables: 
      - M_{lll}
      - M_{ll}
      - M_{T}
      - M_{l_{2}l_{3}}
      - MET
  */
  TCut genBase("genNleptons>=3 && gentri_pt1>30 && gentri_pt2>20 && gentri_pt3>20");
  TCut genParticleBase("genNleptons>=3 && pgentri_pt1>30 && pgentri_pt2>20 && pgentri_pt3>20");
  
  TCut Selection("abs(gentri_mlll-91)>30 && genMET>30 && abs(gentri_mll-91)<10 && gentri_mT > 40");
  TCut genParticleSelection("abs(pgentri_mlll-91)>30 && genMET>30 && abs(pgentri_mll-91)<10 && abs(pgentri_submll-91)>20 && pgentri_mT > 40");
  
  bool dolog=true;
  bool nolog=false;
  
  DrawMatchGenDistribution("gentri_mlll",60,0,300,dolog,"m_{lll}","TrileptonMass",genBase);
  DrawMatchGenDistribution("gentri_mll",60,0,300,dolog,"m_{l_{1}l_{2}}","LeadingDileptonMass",genBase);
  DrawMatchGenDistribution("gentri_mT",50,0,100,nolog,"m_{T}","MT",genBase);
  DrawMatchGenDistribution("gentri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","SubleadingDileptonMass",genBase);
  DrawMatchGenDistribution("gentri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","SubleadingDileptonMass_SF",genBase&&TCut("abs(gentri_id2)==abs(gentri_id3)"));
  DrawMatchGenDistribution("gentri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","SubleadingDileptonMass_OF",genBase&&TCut("abs(gentri_id2)!=abs(gentri_id3)"));
  DrawMatchGenDistribution("genMET",60,0,300,dolog,"gen MET","MET",genBase);
  DrawMatchGenDistribution("abs(gentri_id1)",50,-0.5,49.5,nolog,"id(l_{1})","genID1",genBase);
  DrawMatchGenDistribution("abs(gentri_id2)",50,-0.5,49.5,nolog,"id(l_{2})","genID2",genBase);
  DrawMatchGenDistribution("abs(gentri_id3)",50,-0.5,49.5,nolog,"id(l_{3})","genID3",genBase);
  DrawMatchGenDistribution("abs(gentri_genMID1)",50,-0.5,49.5,nolog,"mother id(l_{1})","genMID1",genBase);
  DrawMatchGenDistribution("abs(gentri_genMID2)",50,-0.5,49.5,nolog,"mother id(l_{2})","genMID2",genBase);
  DrawMatchGenDistribution("abs(gentri_genMID3)",50,-0.5,49.5,nolog,"mother id(l_{3})","genMID3",genBase);
  
  DrawMatchGenDistribution("gentri_mlll",60,0,300,dolog,"m_{lll}","FullSelection/TrileptonMass",genBase&&Selection);
  DrawMatchGenDistribution("gentri_mll",60,0,300,dolog,"m_{l_{1}l_{2}}","FullSelection/LeadingDileptonMass",genBase&&Selection);
  DrawMatchGenDistribution("gentri_mT",50,0,100,nolog,"m_{T}","FullSelection/MT",genBase&&Selection);
  DrawMatchGenDistribution("gentri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","FullSelection/SubleadingDileptonMass",genBase&&Selection);
  DrawMatchGenDistribution("gentri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","FullSelection/SubleadingDileptonMass_SF",genBase&&TCut("abs(gentri_id2)==abs(gentri_id3)")&&Selection);
  DrawMatchGenDistribution("gentri_submll",60,0,300,nolog,"m_{l_{2}l_{3}}","FullSelection/SubleadingDileptonMass_OF",genBase&&TCut("abs(gentri_id2)!=abs(gentri_id3)")&&Selection);
  DrawMatchGenDistribution("genMET",60,0,300,dolog,"gen MET","FullSelection/MET",genBase&&Selection);
  DrawMatchGenDistribution("abs(gentri_id1)",50,-0.5,49.5,nolog,"id(l_{1})","FullSelection/genID1",genBase&&Selection);
  DrawMatchGenDistribution("abs(gentri_id2)",50,-0.5,49.5,nolog,"id(l_{2})","FullSelection/genID2",genBase&&Selection);
  DrawMatchGenDistribution("abs(gentri_id3)",50,-0.5,49.5,nolog,"id(l_{3})","FullSelection/genID3",genBase&&Selection);
  DrawMatchGenDistribution("abs(gentri_genMID1)",50,-0.5,49.5,nolog,"mother id(l_{1})","FullSelection/genMID1",genBase&&Selection);
  DrawMatchGenDistribution("abs(gentri_genMID2)",50,-0.5,49.5,nolog,"mother id(l_{2})","FullSelection/genMID2",genBase&&Selection);
  DrawMatchGenDistribution("abs(gentri_genMID3)",50,-0.5,49.5,nolog,"mother id(l_{3})","FullSelection/genMID3",genBase&&Selection);
  
return;


  GenSFOFCalculator("Basic selection",genBase);
  GenSFOFCalculator("Basic selection, correct Z",genBase&&TCut("gentri_GoodZMatch"));
  GenSFOFCalculator("Basic selection, wrong Z",genBase&&TCut("!gentri_GoodZMatch"));
  GenSFOFCalculator("Outside |m_{23}-Z|<20",genBase&&TCut("abs(gentri_submll-91)>20"));
  GenSFOFCalculator("Outside |m_{23}-Z|<20, Correct Z", genBase&&TCut("abs(gentri_submll-91)>20&&gentri_GoodZMatch"));
  GenSFOFCalculator("Outside |m_{23}-Z|<20, Wrong Z", genBase&&TCut("abs(gentri_submll-91)>20&&!gentri_GoodZMatch"));
  GenSFOFCalculator("Outside |m_{23}-Z|<20, Correct W", genBase&&TCut("abs(gentri_submll-91)>20&&gentri_GoodWMatch"));
  GenSFOFCalculator("Outside |m_{23}-Z|<20, Wrong W", genBase&&TCut("abs(gentri_submll-91)>20&&!gentri_GoodWMatch"));
  GenSFOFCalculator("Outside |m_{23}-Z|<20, Correct Z, correct W", genBase&&TCut("abs(gentri_submll-91)>20&&gentri_GoodZMatch&&gentri_GoodWMatch"));
  GenSFOFCalculator("Outside |m_{23}-Z|<20, Wrong Z or W", genBase&&TCut("abs(gentri_submll-91)>20&&!(gentri_GoodZMatch&&gentri_GoodWMatch)"));
  GenSFOFCalculator("Full selection, Correct Z and W", genBase&&Selection&&TCut("abs(gentri_submll-91)>20&&(gentri_GoodZMatch&&gentri_GoodWMatch)"));
  GenSFOFCalculator("Full selection, Wrong Z or W", genBase&&Selection&&TCut("abs(gentri_submll-91)>20&&!(gentri_GoodZMatch&&gentri_GoodWMatch)"));
//  GenSFOFCalculator("Full selection", genBase&&Selection);
//  GenSFOFCalculator("Full selection with genParticles", genBase&&genParticleSelection&&TCut("abs(gentri_submll-91)>20"));
  
  
}


void AlgorithmTest(TCut BasicCut, TCut Selection) {
  TCut Pure("tri_GoodWMatch && tri_GoodZMatch");
  
  int nbins=38;
  float mll_low=10;
  float mll_hi=200;
  
  TCanvas *can = new TCanvas("can","can");
  TH1F *sfwz   = allsamples.Draw("sfwz"  ,"tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepSF,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  TH1F *ofwz   = allsamples.Draw("ofwz"  ,"tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepOF,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  
  TH1F *sfwzpure = allsamples.Draw("sfwzpure","tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepSF&&Pure,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  TH1F *ofwzpure = allsamples.Draw("ofwzpure","tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepOF&&Pure,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  
  delete can;
  
  dout << "SF WZ events: " << sfwz->Integral() << "   (correctly combined WZ: " << sfwzpure->Integral() << ")" << endl;
  dout << "OF WZ events: " << ofwz->Integral() << "   (correctly combined WZ: " << ofwzpure->Integral() << ")" << endl;
  
  dout << "RATIO : SF / OF = " << sfwz->Integral() / ofwz->Integral() << " +/- " << sfwz->Integral()/ofwz->Integral() * sqrt(1/sfwz->GetEntries() + 1/ofwz->GetEntries()) << endl;
  dout << "  PURE: SF / OF = " << sfwzpure->Integral() / ofwzpure->Integral() << " +/- " << sfwzpure->Integral()/ofwzpure->Integral() * sqrt(1/sfwzpure->GetEntries() + 1/ofwzpure->GetEntries()) << endl;
  dout << "PURITY: " << endl;
  dout << "   SF : " << 100.0*sfwzpure->Integral()/sfwz->Integral() << " +/- " << 100.0*sfwzpure->Integral()/sfwz->Integral() * sqrt(1/sfwzpure->GetEntries() + 1/sfwz->GetEntries()) <<endl;
  dout << "   OF : " << 100.0*ofwzpure->Integral()/ofwz->Integral() << " +/- " << 100.0*ofwzpure->Integral()/ofwz->Integral() * sqrt(1/ofwzpure->GetEntries() + 1/ofwz->GetEntries()) <<endl;
  

  delete sfwz;
  delete ofwz;
  delete sfwzpure;
  delete ofwzpure;
}

void GetResult(TCut BasicCut, TCut Selection) {
    if(TriLeptons::Rsfof_CorrectionFactor<0 || TriLeptons::Rsfof_CorrectionFactor_Error<0) SF_OF_comparison(BasicCut,Selection);
    int nbins=38;
    float mll_low=10;
    float mll_hi=200;
    
    TCanvas *can = new TCanvas("can","can");
    TH1F *sfwz   = SingleLeptonData.Draw("sfwz"  ,"tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepSF,data,luminosity);
    TH1F *ofwz   = SingleLeptonData.Draw("ofwz"  ,"tri_submll",nbins,mll_low,mll_hi,"m_{l_{2}l_{3}} [GeV]", "events", BasicCut&&Selection&&TriLepOF,data,luminosity);
    
    delete can;
    
    dout << " ****** FINAL RESULT ****** " << endl;
    dout << "SF events: " << sfwz->Integral() << " +/- " << sqrt(sfwz->Integral()) << endl;
    dout << "OF events: " << ofwz->Integral() << " +/- " << sqrt(ofwz->Integral()) << endl;
    
    dout << "Correction factor : " << TriLeptons::Rsfof_CorrectionFactor << " +/- " << TriLeptons::Rsfof_CorrectionFactor_Error << endl;
    
    float rawstatuncert=sfwz->Integral()/ofwz->Integral() * sqrt(1/sfwz->GetEntries() + 1/ofwz->GetEntries());
    float rawresult=sfwz->Integral() / ofwz->Integral();
    dout << "UNCORRECTED RATIO : SF / OF = " << rawresult << " +/- " << rawstatuncert << endl;
    dout << "CORRECTED RATIO : SF / OF = " << TriLeptons::Rsfof_CorrectionFactor * rawresult << " +/- " << (TriLeptons::Rsfof_CorrectionFactor * rawresult) * sqrt(rawstatuncert*rawstatuncert/(rawresult*rawresult) + ( TriLeptons::Rsfof_CorrectionFactor_Error * TriLeptons::Rsfof_CorrectionFactor_Error )/(TriLeptons::Rsfof_CorrectionFactor * TriLeptons::Rsfof_CorrectionFactor));
    
    
    delete sfwz;
    delete ofwz;
//    delete can;
}


void PredictionFailurePlot(string truthvariable, string predictionvariable, int nbins, float min, float max, TCut PredictionCut, TCut TruthCut, string saveas, string xlabel) {
  TLegend *leg = make_legend();
  leg->SetY1(0.7);
  TCanvas *can = new TCanvas("can","can");
  TPad *pad = new TPad("pad","pad",0,0,1,1);
  pad->cd();
//    pad->SetLogy(1);
  
  TH1F *Prediction = allsamples.Draw("Prediction",predictionvariable,nbins,min,max,xlabel,"events",PredictionCut,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  TH1F *MCTruth    = allsamples.Draw("MCTruth",   truthvariable,     nbins,min,max,xlabel,"events",TruthCut     ,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  
  Prediction->SetLineColor(TColor::GetColor("#FF0000"));//an evil red
  MCTruth->SetLineColor(TColor::GetColor("#04B404")); // a nice green
  
//  if(Prediction->GetMaximum()>MCTruth->GetMaximum()) MCTruth->SetMaximum(1.3*Prediction->GetMaximum());
  
  if(MCTruth->GetMaximum()>Prediction->GetMaximum()) {
    MCTruth->SetMaximum(MCTruth->GetMaximum()*1.2);
    MCTruth->DrawNormalized("histo");
    Prediction->DrawNormalized("histo,same");
  } else {
    Prediction->SetMaximum(Prediction->GetMaximum()*1.2);
    Prediction->DrawNormalized("histo");
    MCTruth->DrawNormalized("histo,same");
  }
    
  
  leg->AddEntry(Prediction,"prediction (norm.)","l");
  leg->AddEntry(MCTruth,"MC truth (norm.)","l");
  leg->Draw();
  DrawMCPrelim();
  Save_With_Ratio(MCTruth,Prediction,pad->cd(),"WZ_study/PredictionMethodFailure/"+saveas);
  
  delete Prediction;
  delete MCTruth;
  delete can;
  
  
}

void WhyDoesThePredictionMethodFail(TCut PredictionCut, TCut TruthCut) {
  PredictionFailurePlot("tri_submll","tri_badsubmll",50,0,200,PredictionCut,TruthCut,"SubMll","m_{23}");
  PredictionFailurePlot("tri_mT","tri_badmT",50,0,200,PredictionCut,TruthCut,"MT","m_{T}");
  PredictionFailurePlot("tri_mll","tri_badmll",25,0,200,PredictionCut,TruthCut,"Mll","m_{12}");
}
  
  
void ComputeNewMiscombinationRate(TCut origSelection, TCut sBasicCut) {
  TCanvas *can = new TCanvas("can","can");
  
  string sSelection= (const char*)origSelection;
  sSelection= ReplaceAll(sSelection,"tri_","tri_bad");
  sSelection= ReplaceAll(sSelection,"tri_badmlll","tri_mlll");
  cout << "modified " << (const char*)origSelection << " to " << sSelection << endl;
  
  TCut FlavorBlindSelection((sSelection+(string)"&&(abs(tri_id2)!=abs(tri_id3))").c_str()); // making sure we select OF events
  TCut CorrectSelection    (origSelection && TCut("(abs(tri_id2)!=abs(tri_id3))"));        // making sure we select OF events
  TCut OFZ                 ("tri_mll!=tri_badmll");
  
  TH1F *WZ_OF_OFZ = allsamples.Draw("WZ_OFZ"    ,"tri_badsubmll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", FlavorBlindSelection&&OFZ ,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  TH1F *WZ_OF_any = allsamples.Draw("WZ_OF_any" ,"tri_badsubmll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", FlavorBlindSelection      ,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
  
  TH1F *data_OF_OFZ = SingleLeptonData.Draw("data_OFZ"    ,"tri_badsubmll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", FlavorBlindSelection&&OFZ ,data,luminosity);
  TH1F *data_OF_any = SingleLeptonData.Draw("data_OF_any" ,"tri_badsubmll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", FlavorBlindSelection      ,data,luminosity);
  
  dout << "Miscombination rate in WZ sample with selection in place : " << 100 * WZ_OF_OFZ->Integral() / WZ_OF_any->Integral() << " +/- " << 100*(WZ_OF_OFZ->Integral() / WZ_OF_any->Integral())* sqrt( 1.0/WZ_OF_OFZ->GetEntries() + 1.0/WZ_OF_any->GetEntries())<< " %    (" << WZ_OF_OFZ->Integral() << " / " << WZ_OF_any->Integral() << ")" << endl;
  dout << "Miscombination rate in DATA with selection in place : " << 100 * data_OF_OFZ->Integral() / data_OF_any->Integral() << " +/- " << 100*(data_OF_OFZ->Integral() / data_OF_any->Integral()) * sqrt( 1.0/data_OF_OFZ->Integral() + 1.0/data_OF_any->Integral() ) << " %    (" << data_OF_OFZ->Integral() << " / " << data_OF_any->Integral() << ")" << endl;
  
  delete WZ_OF_any;
  delete WZ_OF_OFZ;
  delete can;
}

void ComputeMiscombinationRate(TCut origSelection, TCut sBasicCut) {
/*  TCanvas *can = new TCanvas("can","can");
  //idea: use the flavor-blind version of the algorithm to obtain the yields in our signal region.
  string sSelection= (const char*)origSelection;
  sSelection= ReplaceAll(sSelection,"tri_","tri_bad");
  sSelection= ReplaceAll(sSelection,"tri_badmlll","tri_mlll");
  cout << "modified " << (const char*)origSelection << " to " << sSelection << endl;
  
  TCut BadSelection((sSelection+(string)"&&abs(tri_id2)!=abs(tri_id3)").c_str()); // making sure we select OF events
  TCut GoodSelection(origSelection && TCut("abs(tri_id2)!=abs(tri_id3)")); // making sure we select OF events
  
  TCut BasicCut(sBasicCut&&TCut("tri_id2!=tri_id3"));
  
  TCut TruthBad("!(tri_GoodZMatch&&tri_GoodZMatch)");
  TCut MethodBad("!(tri_id2==tri_badid2&&tri_id3==tri_badid3)");
  
  for(int icut=0;icut<2;icut++) {
    TCut Goodused;
    TCut Badused;
    if(icut==0) {
      Goodused=BasicCut;
      Badused =BasicCut;
      cout << "Preselection only: " << endl;
      cout << "      used : " << (const char*) Goodused << endl;
    } else {
      Goodused=BasicCut&&GoodSelection;
      Badused =BasicCut&&BadSelection;
      cout << "Full selection: " << endl;
      cout << "      Good used : " << (const char*) Goodused << endl;
      cout << "      Bad used : " << (const char*) Goodused << endl;
    }
    
    TH1F *Truth_badcomb = allsamples.Draw("Truth_badcomb"    ,"tri_submll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", Goodused           ,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
    TH1F *Truth_all     = allsamples.Draw("Truth_all"        ,"tri_submll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", Goodused           ,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
    
    TH1F *Events_WZ    = allsamples.Draw("Events_WZ"    ,"tri_submll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", Goodused           ,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
    TH1F *bad_WZ       = allsamples.Draw("bad_WZ"       ,"tri_submll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", Goodused&&TruthBad ,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
    TH1F *MethodBad_WZ = allsamples.Draw("MethodBad_WZ" ,"tri_submll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", Badused&&MethodBad,mc,luminosity,allsamples.FindSample("WZJetsTo3"));
    TH1F *MethodBad_WZg= allsamples.Draw("MethodBad_WZg","tri_submll",1,0.0,1000.0,"m_{l_{2}l_{3}} [GeV]", "events", Badused,mc,luminosity,allsamples.FindSample("WZJetsTo3"));

    cout << "    ********** WZ **********" << endl;
    cout << "    N events :          " << Events_WZ->Integral() << "\t (mc events: " << Events_WZ->GetEntries() << ")" << endl;
    cout << "    Bad combinations:   " << bad_WZ->Integral() << "\t (mc events: " << bad_WZ->GetEntries() << ")" << endl;
    cout << "    Predicted bad com.: " << MethodBad_WZ->Integral() << " / " << MethodBad_WZg->Integral() << "\t (mc events: " << MethodBad_WZ->GetEntries() << ")" << endl;


    delete Events_WZ;
    delete bad_WZ;
    delete MethodBad_WZ;
    delete MethodBad_WZg;
  }
  delete can;
  
  WhyDoesThePredictionMethodFail(BasicCut&&MethodBad,BasicCut&&TruthBad);
*/  
}
  
  
void CarryOutWZStudy() {
  
  WZ_kin_plot("mll",40,0,200,"m_{Z} [GeV]",TCut("id1==id2 && pfJetGoodNum40>=3 && abs(mll-91.2)<20 && pt1>30 && pt2>20 && abs(eta1)<2.4 && abs(eta2)<2.4"),"Crosscheck_Lumi");
  WZ_kin_plot("mll",40,0,200,"m_{Z} [GeV]",TCut("id1==id2 && id1==0 && pfJetGoodNum40>=3 && abs(mll-91.2)<20 && pt1>30 && pt2>20 && abs(eta1)<2.4 && abs(eta2)<2.4"),"Crosscheck_Lumi_ee");
  WZ_kin_plot("mll",40,0,200,"m_{Z} [GeV]",TCut("id1==id2 && id1==1 && pfJetGoodNum40>=3 && abs(mll-91.2)<20 && pt1>30 && pt2>20 && abs(eta1)<2.4 && abs(eta2)<2.4"),"Crosscheck_Lumi_mm");
  
  TCut BasicCut("tri_MatchFound&&leptonNum>2&&tri_pt1>30&&tri_pt2>20&&tri_pt3>20");
  TCut Selection("abs(tri_mlll-91)>30 && tri_mT>40 && met[4] > 30 && abs(tri_mll-91)<10 && tri_submll > 10 && abs(tri_submll-91)>5");

  //OptimizeSelection(BasicCut);
  
  //ComputeMiscombinationRate(Selection, BasicCut);
  
//   ComputeNewMiscombinationRate(Selection, BasicCut);
  
  WZ_kin_plot("tri_mll",40,0,200,"m_{Z} [GeV]",BasicCut,"Inclusive/Zmass");
  WZ_kin_plot("tri_mlll",40,0,400,"m_{ll} [GeV]",BasicCut,"Inclusive/TrileptonMass");
  WZ_kin_plot("tri_mT",40,0,200,"m_{T} [GeV]",BasicCut,"Inclusive/mT");
  WZ_kin_plot("met[4]", 40,0,400,"raw PFMET [GeV]",BasicCut,"Inclusive/MET");
  WZ_kin_plot("tri_pt1",40,0,400,"p_{T}^{1} [GeV]",BasicCut,"Inclusive/Pt1");
  WZ_kin_plot("tri_pt2",40,0,400,"p_{T}^{2} [GeV]",BasicCut,"Inclusive/Pt2");
  WZ_kin_plot("tri_pt3",40,0,400,"p_{T}^{3} [GeV]",BasicCut,"Inclusive/Pt3");
/*
  WZ_kin_plot("tri_mll",40,0,200,"m_{Z} [GeV]",BasicCut&&TCut("abs(tri_mlll-91)>30"),"Mlll_cut/Zmass");
  WZ_kin_plot("tri_mlll",40,0,400,"m_{ll} [GeV]",BasicCut&&TCut("abs(tri_mlll-91)>30"),"Mlll_cut/TrileptonMass");
  WZ_kin_plot("tri_mT",40,0,200,"m_{T} [GeV]",BasicCut&&TCut("abs(tri_mlll-91)>30"),"Mlll_cut/mT");
  WZ_kin_plot("met[4]", 40,0,400,"raw PFMET [GeV]",BasicCut&&TCut("abs(tri_mlll-91)>30"),"Mlll_cut/MET");
  
  WZ_kin_plot("tri_mll",40,0,200,"m_{Z} [GeV]",BasicCut&&TCut("abs(tri_mll-91)<10"),"Zwindow/Zmass");
  WZ_kin_plot("tri_mlll",40,0,400,"m_{ll} [GeV]",BasicCut&&TCut("abs(tri_mll-91)<10"),"Zwindow/TrileptonMass");
  WZ_kin_plot("tri_mT",40,0,200,"m_{T} [GeV]",BasicCut&&TCut("abs(tri_mll-91)<10"),"Zwindow/mT");
  WZ_kin_plot("met[4]", 40,0,400,"raw PFMET [GeV]",BasicCut&&TCut("abs(tri_mll-91)<10"),"Zwindow/MET");
 
  WZ_kin_plot("tri_mll",40,0,200,"m_{Z} [GeV]",BasicCut,"CutFlow/Zmass__0_CutOnNothing");
  WZ_kin_plot("tri_mll",40,0,200,"m_{Z} [GeV]",BasicCut&&TCut("tri_mT>40"),"CutFlow/Zmass__1_CutOn_MT");
  WZ_kin_plot("tri_mll",40,0,200,"m_{Z} [GeV]",BasicCut&&TCut("tri_mT>40&&abs(tri_mlll-91)>30"),"CutFlow/Zmass__2_CutOn_MT_n_Mll");
  WZ_kin_plot("tri_mll",40,0,200,"m_{Z} [GeV]",BasicCut&&TCut("tri_mT>40&&abs(tri_mlll-91)>30&&met[4]>30"),"CutFlow/Zmass__3_CutOn_MT_n_Mll_n_Met");
  WZ_kin_plot("tri_mll",40,0,200,"m_{Z} [GeV]",BasicCut&&TCut("tri_mT>40&&abs(tri_mlll-91)>30&&met[4]>30&&abs(tri_mll-91)<10"),"CutFlow/Zmass__4_CutOn_MT_n_Mll_n_Met_MllWindow");
   
  WZ_kin_plot("tri_mT",40,0,200,"m_{T} [GeV]",BasicCut,"CutFlow/MT__0_CutOnNothing");
  WZ_kin_plot("tri_mT",40,0,200,"m_{T} [GeV]",BasicCut&&TCut("tri_mT>40"),"CutFlow/MT__1_CutOn_MT");
  WZ_kin_plot("tri_mT",40,0,200,"m_{T} [GeV]",BasicCut&&TCut("tri_mT>40&&abs(tri_mlll-91)>30"),"CutFlow/MT__2_CutOn_MT_n_Mll");
  WZ_kin_plot("tri_mT",40,0,200,"m_{T} [GeV]",BasicCut&&TCut("tri_mT>40&&abs(tri_mlll-91)>30&&met[4]>30"),"CutFlow/MT__3_CutOn_MT_n_Mll_n_Met");
  WZ_kin_plot("tri_mT",40,0,200,"m_{T} [GeV]",BasicCut&&TCut("tri_mT>40&&abs(tri_mlll-91)>30&&met[4]>30&&abs(tri_mll-91)<10"),"CutFlow/MT__4_CutOn_MT_n_Mll_n_Met_MllWindow");
   
  WZ_kin_plot("met[4]",40,0,200,"MET [GeV]",BasicCut,"CutFlow/MET__0_CutOnNothing");
  WZ_kin_plot("met[4]",40,0,200,"MET [GeV]",BasicCut&&TCut("tri_mT>40"),"CutFlow/MET__1_CutOn_MT");
  WZ_kin_plot("met[4]",40,0,200,"MET [GeV]",BasicCut&&TCut("tri_mT>40&&abs(tri_mlll-91)>30"),"CutFlow/MET__2_CutOn_MT_n_Mll");
  WZ_kin_plot("met[4]",40,0,200,"MET [GeV]",BasicCut&&TCut("tri_mT>40&&abs(tri_mlll-91)>30&&met[4]>30"),"CutFlow/MET__3_CutOn_MT_n_Mll_n_Met");
  WZ_kin_plot("met[4]",40,0,200,"MET [GeV]",BasicCut&&TCut("tri_mT>40&&abs(tri_mlll-91)>30&&met[4]>30&&abs(tri_mll-91)<10"),"CutFlow/MET__4_CutOn_MT_n_Mll_n_Met_MllWindow");

  */
  WZ_kin_plot("tri_mll",40,0,200,"m_{Z} [GeV]",BasicCut&&Selection,"Zmass");
  WZ_kin_plot("tri_mlll",40,0,400,"m_{lll} [GeV]",BasicCut&&Selection,"TrileptonMass");
  WZ_kin_plot("tri_mT",40,0,200,"m_{T} [GeV]",BasicCut&&Selection,"mT");
  WZ_kin_plot("met[4]", 40,0,400,"raw PFMET [GeV]",BasicCut&&Selection,"MET");
  WZ_kin_plot("tri_pt1",40,0,400,"p_{T}^{1} [GeV]",BasicCut&&Selection,"Pt1");
  WZ_kin_plot("tri_pt2",40,0,400,"p_{T}^{2} [GeV]",BasicCut&&Selection,"Pt2");
  WZ_kin_plot("tri_pt3",40,0,400,"p_{T}^{3} [GeV]",BasicCut&&Selection,"Pt3");
  

  WZ_kin_plot("tri_submll",40,0,200,"m_{l_{2}l_{3}} [GeV]",BasicCut&&Selection&&TriLepSF,"mllsub_SF");
  WZ_kin_plot("tri_submll",40,0,200,"m_{l_{2}l_{3}} [GeV]",BasicCut&&Selection&&TriLepOF,"mllsub_OF");
  WZ_kin_plot("tri_submll",40,0,200,"m_{l_{2}l_{3}} [GeV]",BasicCut&&Selection&&TriLepSF&&TCut("tri_id2==0"),"mllsub_SF_ee");
  WZ_kin_plot("tri_submll",40,0,200,"m_{l_{2}l_{3}} [GeV]",BasicCut&&Selection&&TriLepSF&&TCut("tri_id2==1"),"mllsub_SF_mm");
  WZ_kin_plot("tri_submll",40,0,200,"m_{l_{2}l_{3}} [GeV]",BasicCut&&Selection&&TriLepOF&&TCut("tri_id2==0"),"mllsub_OF_em");
  WZ_kin_plot("tri_submll",40,0,200,"m_{l_{2}l_{3}} [GeV]",BasicCut&&Selection&&TriLepOF&&TCut("tri_id2==1"),"mllsub_OF_me");

  
  SF_OF_comparison(BasicCut,Selection);  
//  generator_SF_OF_comparison();
//  reco_SF_OF_comparison(BasicCut,Selection);

  SF_vs_OF_plots("tri_submll",40,0,200,"m_{l_{2}l_{3}} [GeV]",BasicCut&&Selection,TriLepSF,TriLepOF,mc,"mllsub_MC");
  SF_vs_OF_plots("tri_submll",40,0,200,"m_{l_{2}l_{3}} [GeV]",BasicCut&&Selection,TriLepSF,TriLepOF,data,"mllsub_data");
  AlgorithmTest(BasicCut,Selection);
  GetResult(BasicCut,Selection);
  

}

void Make_Chi2_Resolution_Minimization(TGraph *gr) {
  TCanvas *can2 = new TCanvas("can2","can2");
  can2->cd();
  gStyle->SetOptFit(0);
  gr->Draw("AP*");
  gr->GetXaxis()->SetTitle("MET scale factor");
  gr->GetYaxis()->SetTitle("#chi^{2}");
  gr->GetXaxis()->CenterTitle();
  gr->GetYaxis()->CenterTitle();
  gr->Draw("AP*");
  gr->Fit("pol3");
  TF1 *fit = (TF1*)gr->GetFunction("pol3");
  fit->SetLineWidth(2);
  fit->SetLineColor(kBlue);
  gr->Draw("AP*");
  DrawPrelim();
  CompleteSave(can2,"WZ_study/DYResolution/Result_Chi2");       
  delete can2;
}

void ExtractDYMissingResolution() {
  
  TCanvas *can = new TCanvas("can","can");
  
  TH1F *histoSF = SingleLeptonData.Draw("histoSF","met[4]",50,0,100,"MET [GeV]","events",TCut("id1==id2&&mll>0&&pfJetGoodNum40>=2"),data,luminosity);
  TH1F *histoOF = SingleLeptonData.Draw("histoOF","met[4]",50,0,100,"MET [GeV]","events",TCut("id1!=id2&&mll>0&&pfJetGoodNum40>=2"),data,luminosity);
  THStack histoMC = allsamples.DrawStack("histoMC","met[4]",50,0,100,"MET [GeV]","events",TCut("id1==id2&&mll>0&&pfJetGoodNum40>=2"),mc,luminosity);
  TFile *f = new TFile("file.root","RECREATE");
  
  TH1F *MC_wo_DY = CollapseStack(histoMC);
  TIter nextSF = TIter(histoMC.GetHists());
  TH1F *h;
  while ( h = (TH1F*)nextSF() ) {
    if(Contains(h->GetName(),"Z_Jets")) {
      MC_wo_DY->Add(h,-1);
    }
  }
  MC_wo_DY->Write();
    
  
  int ntests=150;
  TH1F *DYres[ntests];
  
  float min_chi2=-1;
  int best_index=-1;
  stringstream allresults;
  TGraph *gr = new TGraph();
  for(int i=0;i<ntests;i++) {
    
    float smearby=0.15*i/((float)ntests);
    
    stringstream drawvariable;
    drawvariable << "((" << 1+smearby << ")*met[4])";
    cout << "Going to use the variable " << drawvariable.str() << endl;
    stringstream name;
    name << "histo_DY_" << i;
    DYres[i] = allsamples.Draw(name.str().c_str(),drawvariable.str(),50,0,100,"MET [GeV]","events",TCut("id1==id2&&mll>0&&pfJetGoodNum40>=2"),mc,luminosity, allsamples.FindSample("DYJets"));
    DYres[i]->Write();
    
    TPad *pad = new TPad("pad","pad",0,0,1,1);
    pad->cd();
    
    DrawPrelim();
    MC_wo_DY->SetFillColor(kWhite);
    MC_wo_DY->SetLineColor(kBlue);
    TLegend *leg = make_legend();
    leg->SetY1(0.8);
    leg->AddEntry(histoSF,"data","pl");
    leg->AddEntry(MC_wo_DY,"MC","l");
    
    TH1F *sum = (TH1F*)MC_wo_DY->Clone("sum");
    sum->Add(DYres[i]);
    sum->Draw("histo");
    histoSF->Draw("e1,same");
    stringstream Result;
    
    Double_t chi2;
    Int_t ndf,igood;
    Double_t res=0.0;
    Double_t chi2prob = MarcosChi2TestX(histoSF,sum,chi2,ndf,igood,"UW");
    Result << "#splitline{Smeared by : " << std::setprecision(3) << 100*smearby << " %}{#splitline{#chi^{2} / NDF : " << std::setprecision(3) << chi2 << " / " << ndf << "}{ #chi^{2} prob : " << std::setprecision(3)  << chi2prob << "}}";
    
    allresults << "\n" << smearby << ";" << chi2;
    gr->SetPoint(i,smearby,chi2);
    
    if(min_chi2<0||chi2<min_chi2) {
      min_chi2=chi2;
      best_index=i;
    }
    
    TPaveText *PurityPaveText = new TPaveText(0.60, 0.40,0.89, 0.6,"blNDC");
    PurityPaveText->SetFillStyle(4000);
    PurityPaveText->SetBorderSize(0);
    PurityPaveText->SetFillColor(kWhite);
    PurityPaveText->SetTextFont(42);
    PurityPaveText->SetTextSize(0.042);
    PurityPaveText->AddText(Result.str().c_str());
    PurityPaveText->Draw();
    DrawPrelim();
    leg->Draw();
    
    
    stringstream savethis;
    savethis << "WZ_study/DYResolution/Response" << 100*smearby << "Percent";
    
    Save_With_Ratio(histoSF,sum,pad->cd(),savethis.str());
  }
  
  gr->SetTitle("Chi2Diagram");
  gr->SetName("Chi2Diagram");
  
  
  dout << "Result summary: " << endl;
  dout << allresults.str() << endl;
  
  nextSF = TIter(histoMC.GetHists());
  THStack stack;
  while ( h = (TH1F*)nextSF() ) {
    if(!Contains(h->GetName(),"Z_Jets")) {
      stack.Add(h);
    }
  }
  
  
  TPad *pad = new TPad("pad","pad",0,0,1,1);
  pad->cd();

  DYres[best_index]->SetFillColor(kYellow);
  DYres[best_index]->SetLineColor(kYellow);
  
  stack.Add(DYres[best_index]);
  stack.Draw();
  histoSF->Draw("e1,same");
  DrawPrelim();
  Save_With_Ratio(histoSF,CollapseStack(stack),pad->cd(),"WZ_study/DYResolution/Result");
  
  TFile *fr = new TFile("rumba.root","RECREATE");
  gr->Write();
  fr->Close();
  
  delete can;
  
  Make_Chi2_Resolution_Minimization(gr);

}

void WZstudy() {
//  ExtractDYMissingResolution();
  
  cutWeight=TCut("(weight*(weight<1000)*(is_data*tri_MatchFound+(!is_data)*(0.96*(tri_id1==1)+0.90*(tri_id1==0))))");
  switch_overunderflow(true);
  cout << "Essential cut is " << (const char*) essentialcut << endl;
  cout << "Going to set essential cut to new triggers" << endl;
  TCut essential_bkp = essentialcut;
  essentialcut=TCut("mll>5||mll<6");
  cout << "Essential cut is now " << (const char*) essentialcut << endl;
  
  CarryOutWZStudy();
  
  
  cout << "Essential cut is now " << (const char*) essentialcut << endl;
  essentialcut = essential_bkp;
  cout << "Essential cut is now " << (const char*) essentialcut << "   (restored)" << endl;
  switch_overunderflow(false);
  
}

Revision:	1.10
Committed:	Fri Apr 26 06:56:58 2013 UTC (12 years ago) by buchmann
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.9:	+30 -17 lines
Log Message:	- Using dedicated sample collection (SingleLeptonData) to avoid confusion - Added luminosity cross-checks (to make sure SingleLeptonData is complete) - Adapted mT cut from 50 to 40 GeV - Enforce Trigger matching (tri_MatchFound), and use weights (96% for SingleMu, 90% for SingleElectron)