SchieferD/jetcalib/mctopmass_x.cpp

////////////////////////////////////////////////////////////////////////////////
//
// mctopmass_x
// -----------
//
//            11/24/2007 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
////////////////////////////////////////////////////////////////////////////////

#include "utils/cmdline.h"
#include "jetcalib/jetcorr.h"

#include <TROOT.h>
#include <TStyle.h>
#include <TRint.h>
#include <TFile.h>
#include <TTree.h>
#include <TKey.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TLine.h>
#include <TText.h>
#include <TPaveText.h>
#include <TLatex.h>
#include <TH1F.h>
#include <THStack.h>
#include <TF1.h>
#include <TMultiGraph.h>
#include <TGraphErrors.h>
#include <TLorentzVector.h>

#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <vector>
#include <set>
#include <cmath>


//
// define collection of histograms recorded per algorithm
//
class histograms
{
public:
  histograms(const std::string& alg_name,float pt_min,float dr_max)
    : algname(alg_name)
    , ptmin(pt_min)
    , drmax(dr_max)
    , ndilep(0)
    , nlepjets(0)
    , nalljets(0)
    , nlj_njt4(0)
    , naj_njt6(0)
    , ntop_lj(0)
    , ntop_aj(0)
    , nttbar_aj(0)
  {
    TDirectory* currentDir = gDirectory;
    gROOT->cd();
    
    hmt[0]   =new TH1F(("mt_mc_"   +algname).c_str(),"",100,174.9,175.1);
    hmw[0]   =new TH1F(("mw_mc_"   +algname).c_str(),"",100, 80.3, 80.5);
    hmt_lj[0]=new TH1F(("mt_mc_lj_"+algname).c_str(),"",100,174.9,175.1);
    hmw_lj[0]=new TH1F(("mw_mc_lj_"+algname).c_str(),"",100, 80.3, 80.5);
    hmt_aj[0]=new TH1F(("mt_mc_aj_"+algname).c_str(),"",100,174.9,175.1);
    hmw_aj[0]=new TH1F(("mw_mc_aj_"+algname).c_str(),"",100, 80.3, 80.5);
    
    hmt[1]   =new TH1F(("mt_gen_"   +algname).c_str(),"",100,0,300);
    hmw[1]   =new TH1F(("mw_gen_"   +algname).c_str(),"",100,0,150);
    hmt_lj[1]=new TH1F(("mt_gen_lj_"+algname).c_str(),"",100,0,300);
    hmw_lj[1]=new TH1F(("mw_gen_lj_"+algname).c_str(),"",100,0,150);
    hmt_aj[1]=new TH1F(("mt_gen_aj_"+algname).c_str(),"",100,0,300);
    hmw_aj[1]=new TH1F(("mw_gen_aj_"+algname).c_str(),"",100,0,150);
    
    hmt[2]   =new TH1F(("mt_rec_"   +algname).c_str(),"",100,0,300);
    hmw[2]   =new TH1F(("mw_rec_"   +algname).c_str(),"",100,0,150);
    hmt_lj[2]=new TH1F(("mt_rec_lj_"+algname).c_str(),"",100,0,300);
    hmw_lj[2]=new TH1F(("mw_rec_lj_"+algname).c_str(),"",100,0,150);
    hmt_aj[2]=new TH1F(("mt_rec_aj_"+algname).c_str(),"",100,0,300);
    hmw_aj[2]=new TH1F(("mw_rec_aj_"+algname).c_str(),"",100,0,150);
    
    hmt[3]   =new TH1F(("mt_corr_"   +algname).c_str(),"",100,0,300);
    hmw[3]   =new TH1F(("mw_corr_"   +algname).c_str(),"",100,0,150);
    hmt_lj[3]=new TH1F(("mt_corr_lj_"+algname).c_str(),"",100,0,300);
    hmw_lj[3]=new TH1F(("mw_corr_lj_"+algname).c_str(),"",100,0,150);
    hmt_aj[3]=new TH1F(("mt_corr_aj_"+algname).c_str(),"",100,0,300);
    hmw_aj[3]=new TH1F(("mw_corr_aj_"+algname).c_str(),"",100,0,150);

    hmt[4]   =new TH1F(("mt_l5_"   +algname).c_str(),"",100,0,300);
    hmw[4]   =new TH1F(("mw_l5_"   +algname).c_str(),"",100,0,150);
    hmt_lj[4]=new TH1F(("mt_l5_lj_"+algname).c_str(),"",100,0,300);
    hmw_lj[4]=new TH1F(("mw_l5_lj_"+algname).c_str(),"",100,0,150);
    hmt_aj[4]=new TH1F(("mt_l5_aj_"+algname).c_str(),"",100,0,300);
    hmw_aj[4]=new TH1F(("mw_l5_aj_"+algname).c_str(),"",100,0,150);

    hjetmult    = new TH1F(("jetmult_"   +algname).c_str(),"",31,-0.5,30.5);
    hjetmult_lj = new TH1F(("jetmult_lj_"+algname).c_str(),"",31,-0.5,30.5);
    hjetmult_aj = new TH1F(("jetmult_aj_"+algname).c_str(),"",31,-0.5,30.5);

    hseleff_lj  = new TH1F(("seleff_lj_"+algname).c_str(),"",46,5,50);
    hseleff_aj  = new TH1F(("seleff_aj_"+algname).c_str(),"",46,5,50);

    hmw[0]   ->SetXTitle("m_{q#bar{q}'} [GeV]");
    hmw_lj[0]->SetXTitle("m_{q#bar{q}'} [GeV]");
    hmw_aj[0]->SetXTitle("m_{q#bar{q}'} [GeV]");
    hmt[0]   ->SetXTitle("m_{bq#bar{q}'} [GeV]");
    hmt_lj[0]->SetXTitle("m_{bq#bar{q}'} [GeV]");
    hmt_aj[0]->SetXTitle("m_{bq#bar{q}'} [GeV]");

    for (int i=1;i<5;i++) {
      hmw[i]   ->SetXTitle("m_{W} [GeV]");
      hmw_lj[i]->SetXTitle("m_{W} [GeV]");
      hmw_aj[i]->SetXTitle("m_{W} [GeV]");
      hmt[i]   ->SetXTitle("m_{top} [GeV]");
      hmt_lj[i]->SetXTitle("m_{top} [GeV]");
      hmt_aj[i]->SetXTitle("m_{top} [GeV]");
    }
    std::stringstream ssxjetmult;
    ssxjetmult<<"jet multiplicity (p_{T}>"<<ptmin<<")";
    hjetmult   ->SetXTitle(ssxjetmult.str().c_str());
    hjetmult_lj->SetXTitle(ssxjetmult.str().c_str());
    hjetmult_aj->SetXTitle(ssxjetmult.str().c_str());
    hseleff_lj->SetXTitle("p_{T}^{min} [GeV]");
    hseleff_aj->SetXTitle("p_{T}^{min} [GeV]");
    hseleff_lj->SetYTitle("#epsilon(Njet#geq4)");
    hseleff_aj->SetYTitle("#epsilon(Njet#geq6)");
    
    currentDir->cd();
  }
  
  histograms(TDirectory* d,float pt_min,float dr_max)
    : algname(d->GetName())
    , ptmin(pt_min)
    , drmax(dr_max)
  {
    TH1I* hcounts = (TH1I*)d->Get(("counts_"+algname).c_str());
    ndilep   =(unsigned int)hcounts->GetBinContent(1);
    nlepjets =(unsigned int)hcounts->GetBinContent(2);
    nalljets =(unsigned int)hcounts->GetBinContent(3);
    nlj_njt4 =(unsigned int)hcounts->GetBinContent(4);
    naj_njt6 =(unsigned int)hcounts->GetBinContent(5);
    ntop_lj  =(unsigned int)hcounts->GetBinContent(6);
    ntop_aj  =(unsigned int)hcounts->GetBinContent(7);
    nttbar_aj=(unsigned int)hcounts->GetBinContent(8);

    hmt[0]   =(TH1F*)d->Get(("mt_mc_"   +algname).c_str());
    hmw[0]   =(TH1F*)d->Get(("mw_mc_"   +algname).c_str());
    hmt_lj[0]=(TH1F*)d->Get(("mt_mc_lj_"+algname).c_str());
    hmw_lj[0]=(TH1F*)d->Get(("mw_mc_lj_"+algname).c_str());
    hmt_aj[0]=(TH1F*)d->Get(("mt_mc_aj_"+algname).c_str());
    hmw_aj[0]=(TH1F*)d->Get(("mw_mc_aj_"+algname).c_str());
    
    hmt[1]   =(TH1F*)d->Get(("mt_gen_"   +algname).c_str());
    hmw[1]   =(TH1F*)d->Get(("mw_gen_"   +algname).c_str());
    hmt_lj[1]=(TH1F*)d->Get(("mt_gen_lj_"+algname).c_str());
    hmw_lj[1]=(TH1F*)d->Get(("mw_gen_lj_"+algname).c_str());
    hmt_aj[1]=(TH1F*)d->Get(("mt_gen_aj_"+algname).c_str());
    hmw_aj[1]=(TH1F*)d->Get(("mw_gen_aj_"+algname).c_str());
    
    hmt[2]   =(TH1F*)d->Get(("mt_rec_"   +algname).c_str());
    hmw[2]   =(TH1F*)d->Get(("mw_rec_"   +algname).c_str());
    hmt_lj[2]=(TH1F*)d->Get(("mt_rec_lj_"+algname).c_str());
    hmw_lj[2]=(TH1F*)d->Get(("mw_rec_lj_"+algname).c_str());
    hmt_aj[2]=(TH1F*)d->Get(("mt_rec_aj_"+algname).c_str());
    hmw_aj[2]=(TH1F*)d->Get(("mw_rec_aj_"+algname).c_str());
    
    hmt[3]   =(TH1F*)d->Get(("mt_corr_"   +algname).c_str());
    hmw[3]   =(TH1F*)d->Get(("mw_corr_"   +algname).c_str());
    hmt_lj[3]=(TH1F*)d->Get(("mt_corr_lj_"+algname).c_str());
    hmw_lj[3]=(TH1F*)d->Get(("mw_corr_lj_"+algname).c_str());
    hmt_aj[3]=(TH1F*)d->Get(("mt_corr_aj_"+algname).c_str());
    hmw_aj[3]=(TH1F*)d->Get(("mw_corr_aj_"+algname).c_str());

    hmt[4]   =(TH1F*)d->Get(("mt_l5_"   +algname).c_str());
    hmw[4]   =(TH1F*)d->Get(("mw_l5_"   +algname).c_str());
    hmt_lj[4]=(TH1F*)d->Get(("mt_l5_lj_"+algname).c_str());
    hmw_lj[4]=(TH1F*)d->Get(("mw_l5_lj_"+algname).c_str());
    hmt_aj[4]=(TH1F*)d->Get(("mt_l5_aj_"+algname).c_str());
    hmw_aj[4]=(TH1F*)d->Get(("mw_l5_aj_"+algname).c_str());

    hjetmult   =(TH1F*)d->Get(("jetmult_"   +algname).c_str());
    hjetmult_lj=(TH1F*)d->Get(("jetmult_lj_"+algname).c_str());
    hjetmult_aj=(TH1F*)d->Get(("jetmult_aj_"+algname).c_str());
    
    hseleff_lj=(TH1F*)d->Get(("seleff_lj_"+algname).c_str());
    hseleff_aj=(TH1F*)d->Get(("seleff_aj_"+algname).c_str());

    set_directory(gROOT);
  }
  
  bool save_to_file()
  {
    std::stringstream ssfilename;
    ssfilename<<"topmass_pt"<<ptmin<<"_dr"<<drmax<<".root";
    TFile*      f = new TFile(ssfilename.str().c_str(),"UPDATE");
    TDirectory* d = (TDirectory*)f->GetDirectory(algname.c_str());
    if (d==0) d=f->mkdir(algname.c_str());
    d->DeleteAll();
    d->cd();
    TH1I* hcounts = new TH1I(("counts_"+algname).c_str(),"",8,0,7);
    hcounts->SetBinContent(1,ndilep);
    hcounts->SetBinContent(2,nlepjets);
    hcounts->SetBinContent(3,nalljets);
    hcounts->SetBinContent(4,nlj_njt4);
    hcounts->SetBinContent(5,naj_njt6);
    hcounts->SetBinContent(6,ntop_lj);
    hcounts->SetBinContent(7,ntop_aj);
    hcounts->SetBinContent(8,nttbar_aj);
    set_directory(d);
    f->Write();
    set_directory(gROOT);
    f->Close();
    delete f;
    gROOT->cd();
    return true;
  }
  
  void set_directory(TDirectory* d)
  {
    for (int i=0;i<5;i++) {
      hmw[i]   ->SetDirectory(d);
      hmt[i]   ->SetDirectory(d);
      hmw_lj[i]->SetDirectory(d);
      hmt_lj[i]->SetDirectory(d);
      hmw_aj[i]->SetDirectory(d);
      hmt_aj[i]->SetDirectory(d);
    }
    hjetmult   ->SetDirectory(d);
    hjetmult_lj->SetDirectory(d);
    hjetmult_aj->SetDirectory(d);
    hseleff_lj ->SetDirectory(d);
    hseleff_aj ->SetDirectory(d);
  }
  
  std::string  algname;
  float        ptmin;
  float        drmax;
  unsigned int ndilep;
  unsigned int nlepjets;
  unsigned int nalljets;
  unsigned int nlj_njt4;
  unsigned int naj_njt6;
  unsigned int ntop_lj;
  unsigned int ntop_aj;
  unsigned int nttbar_aj;
  
  TH1F* hmt[5];
  TH1F* hmw[5];
  TH1F* hmt_lj[5];
  TH1F* hmw_lj[5];
  TH1F* hmt_aj[5];
  TH1F* hmw_aj[5];
  
  TH1F* hjetmult;
  TH1F* hjetmult_lj;
  TH1F* hjetmult_aj;
  TH1F* hseleff_lj;
  TH1F* hseleff_aj;
};


//
// DEFINE GLOBAL FUNCTIONS
//
histograms* analyze_tree(const std::string& filename,float ptmin,float drmax);

void plot_mcmass(histograms* histos,const std::string& channel="");
void plot_topmass(histograms* histos,const std::string& channel="");
void plot_Wmass(histograms* histos,const std::string& channel="");
void plot_jetmult(histograms* histos,const std::string& channel="");

void table_seleff(std::vector<histograms*>vec_histos,const std::string& channel);
void plot_seleff(std::vector<histograms*>vec_histos,const std::string& channel);
void plot_seleff_vspt(std::vector<histograms*>vec_histos,const std::string& channel);

void plot_Wmass_mean(std::vector<histograms*> vec_histos,
                     const std::string& channel="");
void plot_topmass_mean(std::vector<histograms*> vec_histos,
                       const std::string& channel="");
void plot_Wmass_width(std::vector<histograms*> vec_histos,
                      const std::string& channel="");
void plot_topmass_width(std::vector<histograms*> vec_histos,
                        const std::string& channel="");

TGraphErrors* Wmass_mean(std::vector<histograms*> vec_histos,int level,
                         const std::string& channel="");
TGraphErrors* topmass_mean(std::vector<histograms*> vec_histos,int level,
                           const std::string& channel="");
TGraphErrors* Wmass_width(std::vector<histograms*> vec_histos,int level,
                          const std::string& channel="");
TGraphErrors* topmass_width(std::vector<histograms*> vec_histos,int level,
                            const std::string& channel="");

bool parse_filename(const std::string& filename,float& ptmin,float& drmax);

std::string alg_label(const std::string& alg_name);
Color_t     alg_color(const std::string& alg_name);
Style_t     alg_style(const std::string& alg_name);


using namespace std;


//
// MAIN
//

//______________________________________________________________________________
int main(int argc,char**argv)
{
  cmdline cl;
  cl.parse(argc,argv);
  
  vector<string> input   = cl.get_vector<string>("input");
  float          ptmin   = cl.get_value<float>  ("ptmin", 5.0);
  float          drmax   = cl.get_value<float>  ("drmax", 0.5);
  string         channel = cl.get_value<string> ("channel","");

  bool show_all           = cl.get_value<bool>("show_all",          false);
  bool show_mcmass        = cl.get_value<bool>("show_mcmass",       false);
  bool show_Wmass         = cl.get_value<bool>("show_Wmass",        false);
  bool show_topmass       = cl.get_value<bool>("show_topmass",      false);
  bool show_jetmult       = cl.get_value<bool>("show_jetmult",      false);
  bool show_seleff_table  = cl.get_value<bool>("show_seleff_table", false);
  bool show_seleff        = cl.get_value<bool>("show_seleff",       false);
  bool show_seleff_vspt   = cl.get_value<bool>("show_seleff_vspt",  false);
  bool show_Wmass_mean    = cl.get_value<bool>("show_Wmass_mean",   false);
  bool show_topmass_mean  = cl.get_value<bool>("show_topmass_mean", false);
  bool show_Wmass_width   = cl.get_value<bool>("show_Wmass_width",  false);
  bool show_topmass_width = cl.get_value<bool>("show_topmass_width",false);
  
  if (!channel.empty()&&channel!="lj"&&channel!="aj") {
    cout<<"invalid channel '"<<channel<<"'"<<endl;
    return 0;
  }
  
  if(!cl.check()) return 0;
  argc=1;
  TRint* app = new TRint("mctopmass_x",&argc,argv);
  cl.print();
  
  
  vector<histograms*> vec_histos;
  
  if (input[0].find("topmass")!=string::npos) {
    if (!parse_filename(input[0],ptmin,drmax)) return 0;
    set<string> jetalgs;
    for (unsigned int i=1;i<input.size();i++) {
      unsigned int pos;
      string algname(input[i]);
      while ((pos=algname.find('/')!=string::npos)) algname=algname.substr(pos+1);
      if ((pos=algname.find(".root"))!=string::npos) algname=algname.substr(0,pos);
      jetalgs.insert(algname);
    }
    TFile* f = new TFile(input[0].c_str(),"READ");
    TIter next(f->GetListOfKeys());
    TKey* key;
    while ((key = (TKey*)next())!=0) {
      string classname = key->GetClassName();
      if (classname!="TDirectoryFile") continue;
      TDirectoryFile* d = (TDirectoryFile*)key->ReadObj();
      string dirname = d->GetName();
      if (jetalgs.empty()||jetalgs.find(dirname)!=jetalgs.end())
        vec_histos.push_back(new histograms(d,ptmin,drmax));
    }
  }
  else {
    for (vector<string>::const_iterator it=input.begin();it!=input.end();++it) {
      histograms* histos = analyze_tree(*it,ptmin,drmax);
      if (0!=histos) {
        vec_histos.push_back(histos);
        histos->save_to_file();
      }
    }
  }    

  if (vec_histos.size()==0) return 0;
  
  gStyle->SetOptStat(0);

  vector<histograms*>::iterator it;
  for (it=vec_histos.begin();it!=vec_histos.end();++it) {
    if (show_mcmass ||show_all) plot_mcmass(*it,channel);
    if (show_Wmass  ||show_all) plot_Wmass(*it,channel);
    if (show_topmass||show_all) plot_topmass(*it,channel);
    if (show_jetmult||show_all) plot_jetmult(*it,channel);
  }
  if (show_seleff_table ||show_all) table_seleff(vec_histos,channel);
  if (show_seleff       ||show_all) plot_seleff(vec_histos,channel);
  if (show_seleff_vspt  ||show_all) plot_seleff_vspt(vec_histos,channel);
  if (show_Wmass_mean   ||show_all) plot_Wmass_mean(vec_histos,channel);
  if (show_topmass_mean ||show_all) plot_topmass_mean(vec_histos,channel);
  if (show_Wmass_width  ||show_all) plot_Wmass_width(vec_histos,channel);
  if (show_topmass_width||show_all) plot_topmass_width(vec_histos,channel);
  
  app->Run();
  
  return 0;
}


//
// IMPLEMENTATION OF GLOBAL FUNCTIONS
//

//______________________________________________________________________________
histograms* analyze_tree(const string& filename,float ptmin,float drmax)
{
  TFile* f = new TFile(filename.c_str(),"READ"); if (!f->IsOpen()) return 0;
  TTree* t = (TTree*)f->Get("t");                if (0==t)         return 0;
  
  char  njt;            t->SetBranchAddress("njt",         &njt);
    
  float jte[100];       t->SetBranchAddress("jte",          jte);
  float jtpt[100];      t->SetBranchAddress("jtpt",        jtpt);
  float jteta[100];     t->SetBranchAddress("jteta",      jteta);
  float jtphi[100];     t->SetBranchAddress("jtphi",      jtphi);
  float jtjes[100][1];  t->SetBranchAddress("jtjes",      jtjes);
    
  short jtgenflv[100];  t->SetBranchAddress("jtgenflv",jtgenflv);
  float jtgene[100];    t->SetBranchAddress("jtgene",    jtgene);
  float jtgenpt[100];   t->SetBranchAddress("jtgenpt",  jtgenpt);
  float jtgeneta[100];  t->SetBranchAddress("jtgeneta",jtgeneta);
  float jtgenphi[100];  t->SetBranchAddress("jtgenphi",jtgenphi);
  float jtgendr[100];   t->SetBranchAddress("jtgendr",  jtgendr);
    
  char  decay;          t->SetBranchAddress("decay",     &decay);
  float jtmce[100];     t->SetBranchAddress("jtmce",      jtmce);
  float jtmcpt[100];    t->SetBranchAddress("jtmcpt",    jtmcpt);
  float jtmceta[100];   t->SetBranchAddress("jtmceta",  jtmceta);
  float jtmcphi[100];   t->SetBranchAddress("jtmcphi",  jtmcphi);
  float jtmcdr[100];    t->SetBranchAddress("jtmcdr",    jtmcdr);

  jetcorr l5_b;
  jetcorr l5_c;
  jetcorr l5_uds;
  
  string path = getenv("CMSANA_DIR");
  if (!l5_b  .initialize(path+"/data/jetcorr/b.level5"))   return 0;
  if (!l5_c  .initialize(path+"/data/jetcorr/c.level5"))   return 0;
  if (!l5_uds.initialize(path+"/data/jetcorr/uds.level5")) return 0;
  
  string algname = filename;
  unsigned int pos;
  while ((pos = algname.find("/"))!=string::npos) algname = algname.substr(pos);
  algname = algname.substr(0,algname.find("."));
  
  histograms* histos = new histograms(algname,ptmin,drmax);
    
  int nevts=t->GetEntries();
    
  int ndl=0;
  int nlj=0;
  int naj=0;
    
  int nlj_njt4 =0;
  int naj_njt6 =0;
    
  int ntop_lj  =0;
  int ntop_aj  =0;
  int nttbar_aj=0;

  for (int ievt=0;ievt<nevts;ievt++) {
      
    t->GetEntry(ievt);
      
    if (decay==0) { naj++; }
    if (decay==1) { nlj++; }
    if (decay==2) { ndl++; continue; }
      
    if (decay==1&&njt>=4) {
      float pt = 5.0;
      while (pt<=50.0&&pt<=jtpt[3]) { histos->hseleff_lj->Fill(pt); pt+=1.0; }
    }
      
    if (decay==0&&njt>=6) {
      float pt = 5.0;
      while (pt<=50.0&&pt<=jtpt[5]) { histos->hseleff_aj->Fill(pt); pt+=1.0; }
    }
      
    TLorentzVector* b[5]    ={0,0,0,0,0};
    TLorentzVector* bbar[5] ={0,0,0,0,0};
    TLorentzVector* q1[5]   ={0,0,0,0,0};
    TLorentzVector* qbar1[5]={0,0,0,0,0};
    TLorentzVector* q2[5]   ={0,0,0,0,0};
    TLorentzVector* qbar2[5]={0,0,0,0,0};

    int jetmult(0);

    for (int ijt=0;ijt<njt;ijt++) {

      if (jtpt[ijt]<ptmin) continue;
      jetmult++;
        
      if (jtgenflv[ijt]==0||std::abs(jtgenflv[ijt])>5) continue;
      if (jtgendr[ijt]>drmax) continue;
        
      TLorentzVector* lvmc = new TLorentzVector();
      lvmc->SetPtEtaPhiE(jtmcpt[ijt],jtmceta[ijt],jtmcphi[ijt],jtmce[ijt]);
      TLorentzVector* lvgen = new TLorentzVector();
      lvgen->SetPtEtaPhiE(jtgenpt[ijt],jtgeneta[ijt],jtgenphi[ijt],jtgene[ijt]);
      TLorentzVector* lvrec = new TLorentzVector();
      lvrec->SetPtEtaPhiE(jtpt[ijt],jteta[ijt],jtphi[ijt],jte[ijt]);
      TLorentzVector* lvcorr = new TLorentzVector(*lvrec);
      *lvcorr *= jtjes[ijt][0];
      TLorentzVector* lvl5 = new TLorentzVector(*lvcorr);

      if (jtgenflv[ijt]==+5) {
        if (b[0]==0) {
          *lvl5 *= l5_b.eval(lvl5->Pt(),lvl5->Eta());
          b[0]=lvmc;
          b[1]=lvgen;
          b[2]=lvrec;
          b[3]=lvcorr;
          b[4]=lvl5;
        }
      }
      else if (jtgenflv[ijt]==-5) {
        if (bbar[0]==0) {
          *lvl5 *= l5_b.eval(lvl5->Pt(),lvl5->Eta());
          bbar[0]=lvmc;
          bbar[1]=lvgen;
          bbar[2]=lvrec;
          bbar[3]=lvcorr;
          bbar[4]=lvl5;
        }
      }
      else if (jtgenflv[ijt]==+1) {
        *lvl5 *= l5_uds.eval(lvl5->Pt(),lvl5->Eta());
        if (q2[0]==0) {
          q2[0]=lvmc;
          q2[1]=lvgen;
          q2[2]=lvrec;
          q2[3]=lvcorr;
          q2[4]=lvl5;
        }
      }
      else if (jtgenflv[ijt]==-1) {
        if (qbar1[0]==0) {
          *lvl5 *= l5_uds.eval(lvl5->Pt(),lvl5->Eta());
          qbar1[0]=lvmc;
          qbar1[1]=lvgen;
          qbar1[2]=lvrec;
          qbar1[3]=lvcorr;
          qbar1[4]=lvl5;
        }
      }
      else if (jtgenflv[ijt]==+2)  {
        if (q1[0]==0) {
          *lvl5 *= l5_uds.eval(lvl5->Pt(),lvl5->Eta());
          q1[0]=lvmc;
          q1[1]=lvgen;
          q1[2]=lvrec;
          q1[3]=lvcorr;
          q1[4]=lvl5;
        }
      }
      else if (jtgenflv[ijt]==-2) {
        if (qbar2[0]==0) {
          *lvl5 *= l5_uds.eval(lvl5->Pt(),lvl5->Eta());
          qbar2[0]=lvmc;
          qbar2[1]=lvgen;
          qbar2[2]=lvrec;
          qbar2[3]=lvcorr;
          qbar2[4]=lvl5;
        }
      }
      else if (jtgenflv[ijt]==+3) {
        if (q2[0]==0) {
          *lvl5 *= l5_uds.eval(lvl5->Pt(),lvl5->Eta());
          q2[0]=lvmc;
          q2[1]=lvgen;
          q2[2]=lvrec;
          q2[3]=lvcorr;
          q2[4]=lvl5;
        }
      }
      else if (jtgenflv[ijt]==-3) {
        if (qbar1[0]==0) {
          *lvl5 *= l5_uds.eval(lvl5->Pt(),lvl5->Eta());
          qbar1[0]=lvmc;
          qbar1[1]=lvgen;
          qbar1[2]=lvrec;
          qbar1[3]=lvcorr;
          qbar1[4]=lvl5;
        }
      }
      else if (jtgenflv[ijt]==+4) {
        if (q1[0]==0) {
          *lvl5 *= l5_c.eval(lvl5->Pt(),lvl5->Eta());
          q1[0]=lvmc;
          q1[1]=lvgen;
          q1[2]=lvrec;
          q1[3]=lvcorr;
          q1[4]=lvl5;
        }
      }
      else if (jtgenflv[ijt]==-4) {
        if (qbar2[0]==0) {
          *lvl5 *= l5_c.eval(lvl5->Pt(),lvl5->Eta());
          qbar2[0]=lvmc;
          qbar2[1]=lvgen;
          qbar2[2]=lvrec;
          qbar2[3]=lvcorr;
          qbar2[4]=lvl5;
        }
      }
      else assert(0);
    }
      
    bool has_t    = (b[0]   !=0&&q1[0]!=0&&qbar1[0]!=0);
    bool has_tbar = (bbar[0]!=0&&q2[0]!=0&&qbar2[0]!=0);
      
    if (decay==0&&has_t&&has_tbar) nttbar_aj++;

    assert(decay==0||(!has_t||!has_tbar));
      
    histos->hjetmult->Fill(jetmult);
    if (decay==1) histos->hjetmult_lj->Fill(jetmult);
    if (decay==0) histos->hjetmult_aj->Fill(jetmult);
      
    if (decay==1&&jetmult>=4) nlj_njt4++;
    if (decay==0&&jetmult>=6) naj_njt6++;

    if (has_t) {
      if (decay==1) ntop_lj++;
      if (decay==0) ntop_aj++;
      for (int i=0;i<5;i++) {
        TLorentzVector Wplus = *q1[i] + *qbar1[i];
        TLorentzVector t     = *b[i]  + Wplus;
        histos->hmw[i]->Fill(Wplus.M());
        histos->hmt[i]->Fill(t.M());
        if (decay==1) {
          histos->hmw_lj[i]->Fill(Wplus.M());
          histos->hmt_lj[i]->Fill(t.M());
        }
        else if (decay==0) {
          histos->hmw_aj[i]->Fill(Wplus.M());
          histos->hmt_aj[i]->Fill(t.M());
        }
      }
    }

    if (has_tbar) {
      if (decay==1) ntop_lj++;
      if (decay==0) ntop_aj++;
      for (int i=0;i<5;i++) {
        TLorentzVector Wminus = *q2[i] + *qbar2[i];
        TLorentzVector tbar   = *bbar[i]  + Wminus;
        histos->hmw[i]->Fill(Wminus.M());
        histos->hmt[i]->Fill(tbar.M());
        if (decay==1) {
          histos->hmw_lj[i]->Fill(Wminus.M());
          histos->hmt_lj[i]->Fill(tbar.M());
        }
        else if (decay==0) {
          histos->hmw_aj[i]->Fill(Wminus.M());
          histos->hmt_aj[i]->Fill(tbar.M());
        }
      }
    }
      
    if (b[0])    { for (int i=0;i<5;i++) delete b[i];    } 
    if (bbar[0]) { for (int i=0;i<5;i++) delete bbar[i]; }
    if (q1[0])   { for (int i=0;i<5;i++) delete q1[i];   }
    if (qbar1[0]){ for (int i=0;i<5;i++) delete qbar1[i];}
    if (q2[0])   { for (int i=0;i<5;i++) delete q2[i];   }
    if (qbar2[0]){ for (int i=0;i<5;i++) delete qbar2[i];}
      
  }

  for (int i=1;i<=histos->hseleff_lj->GetNbinsX();i++) {
    float nsel=histos->hseleff_lj->GetBinContent(i);
    float nall=nlj;
    float eff = (nall>0) ? nsel/nall : 1.0;
    float eeff = (nall>0) ? std::sqrt(nsel/nall/nall*(1-nsel/nall)):0.0;
    histos->hseleff_lj->SetBinContent(i,eff);
    histos->hseleff_lj->SetBinError(i,eeff);
  }
  
  for (int i=1;i<=histos->hseleff_aj->GetNbinsX();i++) {
    float nsel=histos->hseleff_aj->GetBinContent(i);
    float nall=naj;
    float eff = (nall>0) ? nsel/nall : 1.0;
    float eeff = (nall>0) ? std::sqrt(nsel/nall/nall*(1-nsel/nall)):0.0;
    histos->hseleff_aj->SetBinContent(i,eff);
    histos->hseleff_aj->SetBinError(i,eeff);
  }

  histos->ndilep    = ndl;
  histos->nlepjets  = nlj;
  histos->nalljets  = naj;
  
  histos->nlj_njt4  = nlj_njt4;
  histos->naj_njt6  = naj_njt6;
  histos->ntop_lj   = ntop_lj;
  histos->ntop_aj   = ntop_aj;
  histos->nttbar_aj = nttbar_aj;

  return histos;
}

//______________________________________________________________________________
void plot_mcmass(histograms* histos,const string& channel)
{
  stringstream sscname;
  string ch = (channel.empty()) ? channel : "_"+channel;
  sscname<<"mcmass"<<ch<<"_"<<histos->algname
         <<"_pt"<<histos->ptmin<<"_dr0"<<histos->drmax*10;
  TCanvas* cmcmass = new TCanvas(sscname.str().c_str(),
                                 sscname.str().c_str(),
                                 1000,500);
  cmcmass->Divide(2,1);
  
  TH1F* hmw =
    (channel.empty()) ? histos->hmw[0] :
    (channel=="lj") ? histos->hmw_lj[0] : histos->hmw_aj[0];
  TH1F* hmt =
    (channel.empty()) ? histos->hmt[0] :
    (channel=="lj") ? histos->hmt_lj[0] : histos->hmt_aj[0];
  
  cmcmass->cd(1);
  gPad->SetLogy();
  hmw->Draw();
  hmw->SetTitle("W mass (parton-level)");
  hmw->SetFillStyle(1001);
  hmw->SetFillColor(kBlue);
  stringstream meanW; meanW<<"Mean: "<<setprecision(4)<<hmw->GetMean();
  stringstream rmsW;  rmsW <<"RMS : "<<setprecision(1)<<hmw->GetRMS();
  TText txtW; txtW.SetTextColor(kBlue);
  txtW.DrawTextNDC(0.23,0.8, meanW.str().c_str());
  txtW.DrawTextNDC(0.23,0.74,rmsW.str().c_str());

  cmcmass->cd(2);
  gPad->SetLogy();
  hmt->Draw();
  hmt->SetTitle("top mass (parton-level)");
  hmt->SetFillStyle(1001);
  hmt->SetFillColor(kRed);
  stringstream meantop; meantop<<"Mean: "<<setprecision(3)<<hmt->GetMean();
  stringstream rmstop;  rmstop <<"RMS : "<<setprecision(1)<<hmt->GetRMS();
  TText txttop; txttop.SetTextColor(kRed);
  txttop.DrawTextNDC(0.23,0.8, meantop.str().c_str());
  txttop.DrawTextNDC(0.23,0.74,rmstop.str().c_str());
  
  cmcmass->Modified();
}


//______________________________________________________________________________
void plot_topmass(histograms* histos,const string& channel)
{
  stringstream sscname;
  string ch = (channel.empty()) ? channel : "_"+channel;
  sscname<<"topmass"<<ch<<"_"<<histos->algname
         <<"_pt"<<histos->ptmin<<"_dr0"<<histos->drmax*10;
  TCanvas* ctopmass = new TCanvas(sscname.str().c_str(),
                                  sscname.str().c_str(),
                                  700,700);
  ctopmass->SetLeftMargin(0.15);
  ctopmass->SetRightMargin(0.05);

  TH1F* hgen  =
    (channel.empty()) ? histos->hmt[1] :
    (channel=="lj") ? histos->hmt_lj[1] : histos->hmt_aj[1];
  TH1F* hrec  =
    (channel.empty()) ? histos->hmt[2] :
    (channel=="lj") ? histos->hmt_lj[2] : histos->hmt_aj[2];
  TH1F* hcorr =
    (channel.empty()) ? histos->hmt[3] :
    (channel=="lj") ? histos->hmt_lj[3] : histos->hmt_aj[3];
  TH1F* hl5 =
    (channel.empty()) ? histos->hmt[4] :
    (channel=="lj") ? histos->hmt_lj[4] : histos->hmt_aj[4];
  
  hgen->Draw();
  hrec->Draw("SAME");
  hcorr->Draw("SAME");
  hl5->Draw("SAME");
  
  hgen->SetTitle("top mass: GEN-REC-CORR-L5");
  hgen ->SetLineStyle(4);
  hrec ->SetLineStyle(3);
  hcorr->SetLineStyle(2);
  hl5  ->SetLineStyle(1);
  hgen ->SetFillStyle(1001);
  hgen ->SetFillColor(kGreen);
  hcorr ->SetFillStyle(1001);
  hcorr ->SetFillColor(kOrange+1);
  hl5  ->SetFillStyle(1001);
  hl5  ->SetFillColor(kRed);
  stringstream genmean; genmean <<"Mean: "<<setprecision(4)<<hgen->GetMean();
  stringstream genrms;  genrms  <<"RMS : "<<setprecision(3)<<hgen->GetRMS();
  stringstream recmean; recmean <<"Mean: "<<setprecision(4)<<hrec->GetMean();
  stringstream recrms;  recrms  <<"RMS : "<<setprecision(3)<<hrec->GetRMS();
  stringstream corrmean;corrmean<<"Mean: "<<setprecision(4)<<hcorr->GetMean();
  stringstream corrrms; corrrms <<"RMS : "<<setprecision(3)<<hcorr->GetRMS();
  stringstream l5mean;  l5mean  <<"Mean: "<<setprecision(4)<<hl5->GetMean();
  stringstream l5rms;   l5rms   <<"RMS : "<<setprecision(3)<<hl5->GetRMS();

  TText gentxt; gentxt.SetTextColor(kGreen);
  TText rectxt; rectxt.SetTextColor(kBlack);
  TText corrtxt;corrtxt.SetTextColor(kOrange+1);
  TText l5txt;  l5txt.SetTextColor(kRed);

  gentxt .DrawTextNDC(0.35,0.80,"GEN");
  gentxt .DrawTextNDC(0.35,0.75,genmean.str().c_str());
  gentxt .DrawTextNDC(0.35,0.70,genrms.str().c_str());
  rectxt .DrawTextNDC(0.17,0.50,"REC");
  rectxt .DrawTextNDC(0.17,0.45,recmean.str().c_str());
  rectxt .DrawTextNDC(0.17,0.40,recrms.str().c_str());
  corrtxt.DrawTextNDC(0.73,0.35,"CORR");
  corrtxt.DrawTextNDC(0.73,0.30,corrmean.str().c_str());
  corrtxt.DrawTextNDC(0.73,0.25,corrrms.str().c_str());
  l5txt  .DrawTextNDC(0.66,0.53,"L5(CORR+FLV)");
  l5txt  .DrawTextNDC(0.66,0.48,l5mean.str().c_str());
  l5txt  .DrawTextNDC(0.66,0.43,l5rms.str().c_str());
  
  TLatex algtxt;
  algtxt.SetNDC();
  algtxt.DrawLatex(0.66,0.885,alg_label(histos->algname).c_str());

  hgen->SetMaximum(1.1*hgen->GetMaximum());
  TLine* mtopline = new TLine(175,hgen->GetMinimum(),175,hgen->GetMaximum());
  mtopline->SetLineWidth(2);
  mtopline->Draw("SAME");
  
  ctopmass->Modified();
}


//______________________________________________________________________________
void plot_Wmass(histograms* histos,const string& channel)
{
  stringstream sscname;
  string ch = (channel.empty()) ? channel : "_"+channel;
  sscname<<"Wmass"<<ch<<"_"<<histos->algname
         <<"_pt"<<histos->ptmin<<"_dr0"<<histos->drmax*10;
  TCanvas* cWmass = new TCanvas(sscname.str().c_str(),
                                sscname.str().c_str(),
                                700,700);
  cWmass->SetLeftMargin(0.15);
  cWmass->SetRightMargin(0.05);
  
  TH1F* hgen  =
    (channel.empty()) ? histos->hmw[1] :
    (channel=="lj") ? histos->hmw_lj[1] : histos->hmw_aj[1];
  TH1F* hrec  =
    (channel.empty()) ? histos->hmw[2] :
    (channel=="lj") ? histos->hmw_lj[2] : histos->hmw_aj[2];
  TH1F* hcorr =
    (channel.empty()) ? histos->hmw[3] :
    (channel=="lj") ? histos->hmw_lj[3] : histos->hmw_aj[3];
  TH1F* hl5 =
    (channel.empty()) ? histos->hmw[4] :
    (channel=="lj") ? histos->hmw_lj[4] : histos->hmw_aj[4];
  
  hgen ->Draw();
  hrec ->Draw("SAME");
  hcorr->Draw("SAME");
  hl5  ->Draw("SAME");
  
  hgen->SetTitle("W mass: GEN-REC-CORR-L5");
  hgen ->SetLineStyle(4);
  hrec ->SetLineStyle(3);
  hcorr->SetLineStyle(2);
  hl5  ->SetLineStyle(1);
  hgen ->SetFillStyle(1001);
  hgen ->SetFillColor(kGreen);
  hcorr ->SetFillStyle(1001);
  hcorr ->SetFillColor(kAzure+1);
  hl5  ->SetFillStyle(1001);
  hl5  ->SetFillColor(kBlue);
  stringstream genmean; genmean <<"Mean: "<<setprecision(3)<<hgen->GetMean();
  stringstream genrms;  genrms  <<"RMS : "<<setprecision(3)<<hgen->GetRMS();
  stringstream recmean; recmean <<"Mean: "<<setprecision(3)<<hrec->GetMean();
  stringstream recrms;  recrms  <<"RMS : "<<setprecision(3)<<hrec->GetRMS();
  stringstream corrmean;corrmean<<"Mean: "<<setprecision(3)<<hcorr->GetMean();
  stringstream corrrms; corrrms <<"RMS : "<<setprecision(3)<<hcorr->GetRMS();
  stringstream l5mean;  l5mean  <<"Mean: "<<setprecision(3)<<hl5->GetMean();
  stringstream l5rms;   l5rms   <<"RMS : "<<setprecision(3)<<hl5->GetRMS();

  TText gentxt;  gentxt.SetTextColor(kGreen);
  TText rectxt;  rectxt.SetTextColor(kBlack);
  TText corrtxt; corrtxt.SetTextColor(kAzure+1);
  TText l5txt;   l5txt.SetTextColor(kBlue);
 
  gentxt .DrawTextNDC(0.30,0.80,"GEN");
  gentxt .DrawTextNDC(0.30,0.75,genmean.str().c_str());
  gentxt .DrawTextNDC(0.30,0.70,genrms.str().c_str());
  rectxt .DrawTextNDC(0.17,0.50,"REC");
  rectxt .DrawTextNDC(0.17,0.45,recmean.str().c_str());
  rectxt .DrawTextNDC(0.17,0.40,recrms.str().c_str());
  corrtxt.DrawTextNDC(0.72,0.35,"CORR");
  corrtxt.DrawTextNDC(0.72,0.30,corrmean.str().c_str());
  corrtxt.DrawTextNDC(0.72,0.25,corrrms.str().c_str());
  l5txt  .DrawTextNDC(0.61,0.53,"L5(CORR+FLV)");
  l5txt  .DrawTextNDC(0.61,0.48,l5mean.str().c_str());
  l5txt  .DrawTextNDC(0.61,0.43,l5rms.str().c_str());
  
  TLatex algtxt;
  algtxt.SetNDC();
  algtxt.DrawLatex(0.66,0.885,alg_label(histos->algname).c_str());

  hgen->SetMaximum(1.1*hgen->GetMaximum());
  TLine* mWline = new TLine(80.42,hgen->GetMinimum(),80.42,hgen->GetMaximum());
  mWline->SetLineWidth(2);
  mWline->Draw("SAME");
  
  cWmass->Modified();
}


//______________________________________________________________________________
void plot_jetmult(histograms* histos,const string& channel)
{
  stringstream sscname;
  string ch = (channel.empty()) ? channel : "_"+channel;
  sscname<<"jetmult"<<ch<<"_"<<histos->algname
         <<"_pt"<<histos->ptmin<<"_dr0"<<histos->drmax*10;
  TCanvas* cjetmult = new TCanvas(sscname.str().c_str(),
                                  sscname.str().c_str(),
                                  700,700);
  
  TH1F* h(0);
  if (channel.empty()) h=histos->hjetmult;
  else if (channel=="lj") h=histos->hjetmult_lj;
  else if (channel=="aj") h=histos->hjetmult_aj;
  h->Draw("BAR");
  h->SetTitle("jet multiplicity");
  
  stringstream ssmean; ssmean<<"Mean: "<<setprecision(2)<<h->GetMean();
  TText meantxt; meantxt.DrawTextNDC(0.6,0.6,ssmean.str().c_str());
  
  TLatex algtxt;
  algtxt.SetNDC();
  algtxt.DrawLatex(0.62,0.885,alg_label(histos->algname).c_str());
  if (!channel.empty()) {
    string  ch  = (channel=="lj") ? "lepton+jets" : "alljets";
    Color_t clr = (channel=="lj") ? kYellow-9 : kGreen-9;
    TPaveText* ptch = new TPaveText(0.6,0.7,0.85,0.8,"brNDC");
    ptch->SetFillColor(clr);
    ptch->AddText(ch.c_str());
    ptch->Draw();
  }
  cjetmult->Modified();
}


//______________________________________________________________________________
void table_seleff(vector<histograms*>vec_histos,const string& channel)
{
  if (channel.empty()||channel=="lj") {
    float         nlj = vec_histos[0]->nlepjets;
    vector<float> eff_njt;
    vector<float> eff_bqq;
    for (vector<histograms*>::iterator it=vec_histos.begin();
         it!=vec_histos.end();++it) {
      eff_njt.push_back(100.*(*it)->nlj_njt4/nlj);
      eff_bqq.push_back(100.*(*it)->ntop_lj/(float)((*it)->nlj_njt4));
    }
    cout<<"lepton+jets, eff(njt>=4):"<<endl;
    for (vector<float>::iterator it=eff_njt.begin();it!=eff_njt.end();++it)
      cout<<" & "<<setprecision(3)<<(*it)<<flush;
    cout<<" \\\\\n"<<endl;
    cout<<"lepton+jets, eff(t->bqq'):"<<endl;
    for (vector<float>::iterator it=eff_bqq.begin();it!=eff_bqq.end();++it)
      cout<<" & "<<setprecision(3)<<(*it)<<flush;
    cout<<" \\\\\n"<<endl;
  }
  if (channel.empty()||channel=="aj") {
    float         naj = vec_histos[0]->nalljets;
    vector<float> eff_njt;
    vector<float> eff_bqq;
    vector<float> eff_bbqqqq;
    for (vector<histograms*>::iterator it=vec_histos.begin();
         it!=vec_histos.end();++it) {
      eff_njt.push_back(100.*(*it)->naj_njt6/naj);
      eff_bqq.push_back(100.*(*it)->ntop_aj/2./(float)((*it)->naj_njt6));
      eff_bbqqqq.push_back(100.*(*it)->nttbar_aj/(float)((*it)->naj_njt6));
    }
    cout<<"alljets, eff(njt>=6):"<<endl;
    for (vector<float>::iterator it=eff_njt.begin();it!=eff_njt.end();++it)
      cout<<" & "<<setprecision(3)<<(*it)<<flush;
    cout<<" \\\\\n"<<endl;
    cout<<"alljets, eff(t->bqq'):"<<endl;
    for (vector<float>::iterator it=eff_bqq.begin();it!=eff_bqq.end();++it)
      cout<<" & "<<setprecision(3)<<(*it)<<flush;
    cout<<" \\\\\n"<<endl;
    cout<<"alljets, eff(ttbar->bbqq'qq'):"<<endl;
    for (vector<float>::iterator it=eff_bbqqqq.begin();it!=eff_bbqqqq.end();++it)
      cout<<" & "<<setprecision(3)<<(*it)<<flush;
    cout<<" \\\\\n"<<endl;
  }

}


//______________________________________________________________________________
void plot_seleff(vector<histograms*> vec_histos,const string& channel)
{
  stringstream sscname_lj;
  sscname_lj<<"seleff_lj_pt"<<vec_histos[0]->ptmin<<"_dr0"<<vec_histos[0]->drmax*10;
  stringstream sscname_aj;
  sscname_aj<<"seleff_aj_pt"<<vec_histos[0]->ptmin<<"_dr0"<<vec_histos[0]->drmax*10;
  
  if (channel.empty()||channel=="lj") {
    TCanvas* cseleff_lj = new TCanvas(sscname_lj.str().c_str(),
                                      sscname_lj.str().c_str(),
                                      1000,350);
    cseleff_lj->SetLogy();
    cseleff_lj->SetGridy();
    cseleff_lj->SetLeftMargin(0.12);
    cseleff_lj->SetRightMargin(0.05);
    cseleff_lj->SetTopMargin(0.23);
    
    TLegend* leg_lj = new TLegend(0.50,0.785,0.95,0.99);
    leg_lj->SetFillColor(10);
    
    int nbins = vec_histos.size()*2+1;
    TH1F*    hnsel_lj=new TH1F("nsel_lj","",nbins,0,nbins+1);
    TH1F*    hntop_lj=new TH1F("ntop_lj","",nbins,0,nbins+1);
    THStack* stack_lj=new THStack("seleff_lj","#epsilon(t#rightarrow bq#bar{q}')");
    stack_lj->Add(hntop_lj);
    stack_lj->Add(hnsel_lj);

    hnsel_lj->SetFillStyle(1001);
    hntop_lj->SetFillStyle(1001);
    hnsel_lj->SetFillColor(kBlue+2);
    hntop_lj->SetFillColor(kYellow);

    for (unsigned int i=0;i<vec_histos.size();i++) {
      hnsel_lj->SetBinContent(i*2+2,vec_histos[i]->nlj_njt4);
      hntop_lj->SetBinContent(i*2+2,vec_histos[i]->ntop_lj);
    }
    stack_lj->Draw();
    
    float ptmin = vec_histos[0]->ptmin;
    stringstream ssnsel;ssnsel<<"N(jet)#geq4 (jet pT#geq"<<ptmin<<")";
    stringstream ssntop;ssntop<<"N(t#rightarrow bq#bar{q}') "
                              <<"matched (jet pT#geq"<<ptmin<<")";
    leg_lj->AddEntry(hnsel_lj,ssnsel.str().c_str(),"F");
    leg_lj->AddEntry(hntop_lj,ssntop.str().c_str(),"F");
    leg_lj->Draw();

    for (unsigned int i=0;i<vec_histos.size();i++)
      stack_lj->GetHistogram()->GetXaxis()
        ->SetBinLabel(i*2+2,alg_label(vec_histos[i]->algname).c_str());
    stack_lj->GetHistogram()->LabelsOption("u");    
    stack_lj->GetHistogram()->GetXaxis()->SetLabelOffset(0.01);
    stack_lj->GetHistogram()->GetXaxis()->SetLabelSize(0.065);

    TPaveText* pt_lj = new TPaveText(0.12,0.82,0.3,0.90,"brNDC");
    pt_lj->SetFillColor(kYellow-9);
    pt_lj->AddText("lepton+jets");
    pt_lj->Draw();
    
    cseleff_lj->Modified();
  }
  if (channel.empty()||channel=="aj") {
    TCanvas* cseleff_aj = new TCanvas(sscname_aj.str().c_str(),
                                      sscname_aj.str().c_str(),
                                      1000,350);
    cseleff_aj->SetLogy();
    cseleff_aj->SetGridy();
    cseleff_aj->SetLeftMargin(0.12);
    cseleff_aj->SetRightMargin(0.05);
    cseleff_aj->SetTopMargin(0.23);
    
    TLegend* leg_aj = new TLegend(0.50,0.785,0.95,0.99);
    leg_aj->SetFillColor(10);
    
    int nbins = vec_histos.size()*2+1;
    TH1F*    hnsel_aj  =new TH1F("nsel_aj",  "",nbins,0,nbins+1);
    TH1F*    hntop_aj  =new TH1F("ntop_aj",  "",nbins,0,nbins+1);
    TH1F*    hnttbar_aj=new TH1F("nttbar_aj","",nbins,0,nbins+1);
    THStack* stack_aj  =new THStack("seleff_aj","#epsilon(t#rightarrow bq#bar{q}')");
    stack_aj->Add(hnttbar_aj);
    stack_aj->Add(hntop_aj);
    stack_aj->Add(hnsel_aj);

    hnsel_aj  ->SetFillStyle(1001);
    hntop_aj  ->SetFillStyle(1001);
    hnttbar_aj->SetFillStyle(1001);
    hnsel_aj  ->SetFillColor(kBlue+2);
    hntop_aj  ->SetFillColor(kYellow);
    hnttbar_aj->SetFillColor(kRed);

    for (unsigned int i=0;i<vec_histos.size();i++) {
      hnsel_aj  ->SetBinContent(i*2+2,vec_histos[i]->nlj_njt4);
      hntop_aj  ->SetBinContent(i*2+2,vec_histos[i]->ntop_aj);
      hnttbar_aj->SetBinContent(i*2+2,vec_histos[i]->nttbar_aj);
    }
    stack_aj->Draw();
    
    float ptmin = vec_histos[0]->ptmin;
    stringstream ssnsel;  ssnsel<<"N(jet)#geq6 (jet pT#geq"<<ptmin<<")";
    stringstream ssntop;  ssntop<<"N(t#rightarrow bq#bar{q}') "
                                <<"matched (jet pT#geq"<<ptmin<<")";
    stringstream ssnttbar;ssnttbar<<"N(t#bar{t}#rightarrow "
                                  <<"b#bar{b}q#bar{q}'q#bar{q}') "
                                  <<"matched (jet pT#geq"<<ptmin<<")";
    leg_aj->AddEntry(hnsel_aj,  ssnsel  .str().c_str(),"F");
    leg_aj->AddEntry(hntop_aj,  ssntop  .str().c_str(),"F");
    leg_aj->AddEntry(hnttbar_aj,ssnttbar.str().c_str(),"F");
    leg_aj->Draw();

    for (unsigned int i=0;i<vec_histos.size();i++)
      stack_aj->GetHistogram()->GetXaxis()
        ->SetBinLabel(i*2+2,alg_label(vec_histos[i]->algname).c_str());
    stack_aj->GetHistogram()->LabelsOption("u");    
    stack_aj->GetHistogram()->GetXaxis()->SetLabelOffset(0.01);
    stack_aj->GetHistogram()->GetXaxis()->SetLabelSize(0.065);

    TPaveText* pt_aj = new TPaveText(0.12,0.82,0.3,0.90,"brNDC");
    pt_aj->SetFillColor(kGreen-9);
    pt_aj->AddText("alljets");
    pt_aj->Draw();
    
    cseleff_aj->Modified();
  }
}


//______________________________________________________________________________
void plot_seleff_vspt(vector<histograms*> vec_histos,const string& channel)
{
  stringstream sscname_lj;
  sscname_lj<<"seleff_vspt_lj_dr0"<<vec_histos[0]->drmax*10;
  stringstream sscname_aj;
  sscname_aj<<"seleff_vspt_aj_dr0"<<vec_histos[0]->drmax*10;

  std::vector<histograms*>::iterator it;
  
  if (channel.empty()||channel=="lj") {
    TCanvas* cseleff_lj = new TCanvas(sscname_lj.str().c_str(),
                                      sscname_lj.str().c_str(),
                                      700,700);
    TLegend* leg_lj = new TLegend(0.55,0.85,0.85,0.85-vec_histos.size()*0.045);
    leg_lj->SetFillColor(10);
    string draw_option="PEX0";
    for (it=vec_histos.begin();it!=vec_histos.end();++it) {
      TH1F* h = (*it)->hseleff_lj;
      h->SetTitle("Jet Selection Effificency");
      h->SetLineColor(alg_color((*it)->algname));
      h->SetMarkerColor(alg_color((*it)->algname));
      h->SetMarkerStyle(alg_style((*it)->algname));
      h->Draw(draw_option.c_str());
      if (draw_option=="PEX0") draw_option+="SAME";
      leg_lj->AddEntry(h,alg_label((*it)->algname).c_str());
    }
    leg_lj->Draw();

    TPaveText* pt_lj = new TPaveText(0.67,0.9,0.9,0.98,"brNDC");
    pt_lj->SetFillColor(kYellow-9);
    pt_lj->AddText("lepton+jets");
    pt_lj->Draw();

    cseleff_lj->Modified();
  }
  if (channel.empty()||channel=="aj") {
    TCanvas* cseleff_aj = new TCanvas(sscname_aj.str().c_str(),
                                      sscname_aj.str().c_str(),
                                      700,700);
    TLegend* leg_aj = new TLegend(0.55,0.85,0.85,0.85-vec_histos.size()*0.045);
    leg_aj->SetFillColor(10);
    string draw_option="PEX0";
    for (it=vec_histos.begin();it!=vec_histos.end();++it) {
      TH1F* h = (*it)->hseleff_aj;
      h->SetTitle("Jet Selection Efficiency");
      h->SetLineColor(alg_color((*it)->algname));
      h->SetMarkerColor(alg_color((*it)->algname));
      h->SetMarkerStyle(alg_style((*it)->algname));
      h->Draw(draw_option.c_str());
      if (draw_option=="PEX0") draw_option+="SAME";
      leg_aj->AddEntry(h,alg_label((*it)->algname).c_str());
    }
    leg_aj->Draw();

    TPaveText* pt_aj = new TPaveText(0.7,0.9,0.9,0.98,"brNDC");
    pt_aj->SetFillColor(kGreen-9);
    pt_aj->AddText("alljets");
    pt_aj->Draw();

    cseleff_aj->Modified();
  }
}


//______________________________________________________________________________
void plot_Wmass_mean(std::vector<histograms*> vec_histos,
                     const std::string& channel)
{
  stringstream sscname;
  string ch = (channel.empty()) ? channel : "_"+channel;
  sscname<<"WmassMean"<<ch<<"_pt"<<vec_histos[0]->ptmin
         <<"_dr0"<<vec_histos[0]->drmax*10;
  TCanvas*    cWmean = new TCanvas(sscname.str().c_str(),
                                   sscname.str().c_str(),
                                   1000,350);
  cWmean->SetTopMargin(0.09);
  cWmean->SetLeftMargin(0.1);
  cWmean->SetRightMargin(0.15);
  
  TMultiGraph*    mg = new TMultiGraph();
  TLegend*       leg = new TLegend(0.87,0.17,0.98,0.91);
  leg->SetFillColor(10);
  
  TGraphErrors* ggen = Wmass_mean(vec_histos,1,channel);
  TGraphErrors* grec = Wmass_mean(vec_histos,2,channel);
  TGraphErrors* gcorr= Wmass_mean(vec_histos,3,channel);
  TGraphErrors* gl5  = Wmass_mean(vec_histos,4,channel);

  mg->Add(ggen);
  mg->Add(grec);
  mg->Add(gcorr);
  mg->Add(gl5);
  
  ggen ->SetLineWidth(2);
  grec ->SetLineWidth(2);
  gcorr->SetLineWidth(2);
  gl5  ->SetLineWidth(2);

  ggen->SetLineColor(kGreen+1);
  ggen->SetMarkerColor(kGreen+1);
  grec->SetLineColor(kBlue+2);
  grec->SetMarkerColor(kBlue+2);
  gcorr->SetLineColor(kMagenta);
  gcorr->SetMarkerColor(kMagenta);
  gl5->SetLineColor(kRed);
  gl5->SetMarkerColor(kRed);
  
  ggen->SetMarkerStyle(kFullStar);
  ggen->SetMarkerSize(1.5);
  grec->SetMarkerStyle(kFullSquare);
  gcorr->SetMarkerStyle(kOpenCircle);
  gl5->SetMarkerStyle(kFullCircle);
  
  leg->AddEntry(gcorr,"CORR","lp");
  leg->AddEntry(gl5,  "L5",  "lp");
  leg->AddEntry(ggen, "GEN", "lp");
  leg->AddEntry(grec, "REC", "lp");

  mg->Draw("AP");
  mg->GetHistogram()->SetTitle("W mass (mean)");
  mg->GetHistogram()->SetYTitle("m_{W} [GeV]");
  
  for (unsigned int i=0;i<vec_histos.size();i++) {
    int ibin=0; while ((i+1)>mg->GetHistogram()->GetBinLowEdge(ibin)) { ibin++;}
    mg->GetHistogram()->GetXaxis()
      ->SetBinLabel(ibin,alg_label(vec_histos[i]->algname).c_str());
  }
  mg->GetHistogram()->GetXaxis()->LabelsOption("u");
  mg->GetHistogram()->GetXaxis()->SetLabelOffset(0.01);
  mg->GetHistogram()->GetXaxis()->SetLabelSize(0.065);
  mg->GetHistogram()->GetYaxis()->SetTitleOffset(0.75);
  leg->Draw();

  TLatex mWtxt;
  mWtxt.DrawLatex(0.5,81.0,"m_{W}=80.42 GeV");

  TLatex pttxt;
  stringstream sspt; sspt<<"jet p_{T}#geq "<<vec_histos[0]->ptmin<<" GeV";
  pttxt.SetNDC();
  pttxt.DrawLatex(0.7,0.92,sspt.str().c_str());

  TAxis* xaxis = mg->GetHistogram()->GetXaxis();
  TLine* mWline=new TLine(xaxis->GetXmin(),80.42,xaxis->GetXmax(),80.42);
  mWline->SetLineWidth(2);
  mWline->Draw("SAME");

  cWmean->Modified();
}


//______________________________________________________________________________
void plot_topmass_mean(std::vector<histograms*> vec_histos,
                       const std::string& channel)
{
  stringstream sscname;
  string ch = (channel.empty()) ? channel : "_"+channel;
  sscname<<"topmassMean"<<ch<<"_pt"<<vec_histos[0]->ptmin
         <<"_dr0"<<vec_histos[0]->drmax*10;
  TCanvas* ctopmean = new TCanvas(sscname.str().c_str(),
                                  sscname.str().c_str(),
                                  1000,350);
  
  ctopmean->SetTopMargin(0.09);
  ctopmean->SetLeftMargin(0.1);
  ctopmean->SetRightMargin(0.15);
  
  TMultiGraph*    mg = new TMultiGraph();
  TLegend*       leg = new TLegend(0.87,0.17,0.98,0.91);
  leg->SetFillColor(10);

  TGraphErrors* ggen = topmass_mean(vec_histos,1,channel);
  TGraphErrors* grec = topmass_mean(vec_histos,2,channel);
  TGraphErrors* gcorr= topmass_mean(vec_histos,3,channel);
  TGraphErrors* gl5  = topmass_mean(vec_histos,4,channel);

  mg->Add(ggen);
  mg->Add(grec);
  mg->Add(gcorr);
  mg->Add(gl5);

  ggen ->SetLineWidth(2);
  grec ->SetLineWidth(2);
  gcorr->SetLineWidth(2);
  gl5  ->SetLineWidth(2);
  
  ggen->SetLineColor(kGreen+1);
  ggen->SetMarkerColor(kGreen+1);
  grec->SetLineColor(kBlue+2);
  grec->SetMarkerColor(kBlue+2);
  gcorr->SetLineColor(kMagenta);
  gcorr->SetMarkerColor(kMagenta);
  gl5->SetLineColor(kRed);
  gl5->SetMarkerColor(kRed);
  
  ggen->SetMarkerStyle(kFullStar);
  ggen->SetMarkerSize(1.5);
  grec->SetMarkerStyle(kFullSquare);
  gcorr->SetMarkerStyle(kOpenCircle);
  gl5->SetMarkerStyle(kFullCircle);
  
  leg->AddEntry(gcorr,"CORR","lp");
  leg->AddEntry(gl5,  "L5",  "lp");
  leg->AddEntry(ggen, "GEN", "lp");
  leg->AddEntry(grec, "REC", "lp");

  mg->Draw("AP");
  mg->GetHistogram()->SetTitle("top mass (mean)");
  mg->GetHistogram()->SetYTitle("m_{top} [GeV]");

  for (unsigned int i=0;i<vec_histos.size();i++) {
    int ibin=0; while ((i+1)>mg->GetHistogram()->GetBinLowEdge(ibin)) { ibin++;}
    mg->GetHistogram()->GetXaxis()
      ->SetBinLabel(ibin,alg_label(vec_histos[i]->algname).c_str());
  }
  mg->GetHistogram()->GetXaxis()->LabelsOption("u");
  mg->GetHistogram()->GetXaxis()->SetLabelOffset(0.01);
  mg->GetHistogram()->GetXaxis()->SetLabelSize(0.065);
  mg->GetHistogram()->GetYaxis()->SetTitleOffset(0.75);
  leg->Draw();
  
  TLatex mtoptxt;
  mtoptxt.DrawLatex(0.5,178.6,"m_{top}=175 GeV");
  
  TLatex pttxt;
  stringstream sspt; sspt<<"jet p_{T}#geq "<<vec_histos[0]->ptmin<<" GeV";
  pttxt.SetNDC();
  pttxt.DrawLatex(0.7,0.92,sspt.str().c_str());

  TAxis* xaxis = mg->GetHistogram()->GetXaxis();
  TLine* mtopline = new TLine(xaxis->GetXmin(),175.,xaxis->GetXmax(),175.);
  mtopline->SetLineWidth(2);
  mtopline->Draw("SAME");
  
  ctopmean->Modified();
}


//______________________________________________________________________________
void plot_Wmass_width(std::vector<histograms*> vec_histos,
                      const std::string& channel)
{
  stringstream sscname;
  string ch = (channel.empty()) ? channel : "_"+channel;
  sscname<<"WmassWidth"<<ch<<"_pt"<<vec_histos[0]->ptmin
         <<"_dr0"<<vec_histos[0]->drmax*10;
  TCanvas* cWwidth = new TCanvas(sscname.str().c_str(),
                                 sscname.str().c_str(),
                                 1000,350);
  
  cWwidth->SetTopMargin(0.09);
  cWwidth->SetLeftMargin(0.1);
  cWwidth->SetRightMargin(0.15);
  
  TMultiGraph*    mg = new TMultiGraph();
  TLegend*       leg = new TLegend(0.87,0.17,0.98,0.91);
  leg->SetFillColor(10);

  TGraphErrors* ggen = Wmass_width(vec_histos,1,channel);
  TGraphErrors* grec = Wmass_width(vec_histos,2,channel);
  TGraphErrors* gcorr= Wmass_width(vec_histos,3,channel);
  TGraphErrors* gl5  = Wmass_width(vec_histos,4,channel);

  mg->Add(ggen);
  mg->Add(grec);
  mg->Add(gcorr);
  mg->Add(gl5);

ggen->SetLineColor(kGreen+1);
ggen->SetMarkerColor(kGreen+1);
grec->SetLineColor(kBlue+2);
  grec->SetMarkerColor(kBlue+2);
  gcorr->SetLineColor(kMagenta);
  gcorr->SetMarkerColor(kMagenta);
  gl5->SetLineColor(kRed);
  gl5->SetMarkerColor(kRed);
  
  ggen->SetMarkerStyle(kFullStar);
  ggen->SetMarkerSize(1.5);
  grec->SetMarkerStyle(kFullSquare);
  gcorr->SetMarkerStyle(kOpenCircle);
  gl5->SetMarkerStyle(kFullCircle);
  
  leg->AddEntry(gcorr,"CORR","lp");
  leg->AddEntry(gl5,  "L5",  "lp");
  leg->AddEntry(ggen, "GEN", "lp");
  leg->AddEntry(grec, "REC", "lp");

  mg->Draw("AP");
  mg->GetHistogram()->SetTitle("W mass width (RMS)");
  mg->GetHistogram()->SetYTitle("RMS(m_{W}) [GeV]");

  for (unsigned int i=0;i<vec_histos.size();i++) {
    int ibin=0; while ((i+1)>mg->GetHistogram()->GetBinLowEdge(ibin)) { ibin++;}
    mg->GetHistogram()->GetXaxis()
      ->SetBinLabel(ibin,alg_label(vec_histos[i]->algname).c_str());
  }
  mg->GetHistogram()->GetXaxis()->LabelsOption("u");
  mg->GetHistogram()->GetXaxis()->SetLabelOffset(0.01);
  mg->GetHistogram()->GetXaxis()->SetLabelSize(0.065);
  mg->GetHistogram()->GetYaxis()->SetTitleOffset(0.75);
  leg->Draw();

  TLatex pttxt;
  stringstream sspt; sspt<<"jet p_{T}#geq "<<vec_histos[0]->ptmin<<" GeV";
  pttxt.SetNDC();
  pttxt.DrawLatex(0.7,0.92,sspt.str().c_str());

  cWwidth->Modified();
}


//______________________________________________________________________________
void plot_topmass_width(std::vector<histograms*> vec_histos,
                        const std::string& channel)
{
  stringstream sscname;
  string ch = (channel.empty()) ? channel : "_"+channel;
  sscname<<"topmassWidth"<<ch<<"_pt"<<vec_histos[0]->ptmin
         <<"_dr0"<<vec_histos[0]->drmax*10;
  TCanvas* ctopwidth = new TCanvas(sscname.str().c_str(),
                                   sscname.str().c_str(),
                                   1000,350);

  ctopwidth->SetTopMargin(0.09);
  ctopwidth->SetLeftMargin(0.1);
  ctopwidth->SetRightMargin(0.15);

  TMultiGraph*    mg = new TMultiGraph();
  TLegend*       leg = new TLegend(0.87,0.17,0.98,0.91);
  leg->SetFillColor(10);

  TGraphErrors* ggen = topmass_width(vec_histos,1,channel);
  TGraphErrors* grec = topmass_width(vec_histos,2,channel);
  TGraphErrors* gcorr= topmass_width(vec_histos,3,channel);
  TGraphErrors* gl5  = topmass_width(vec_histos,4,channel);

  mg->Add(ggen);
  mg->Add(grec);
  mg->Add(gcorr);
  mg->Add(gl5);
  
  ggen->SetLineColor(kGreen+1);
  ggen->SetMarkerColor(kGreen+1);
  grec->SetLineColor(kBlue+2);
  grec->SetMarkerColor(kBlue+2);
  gcorr->SetLineColor(kMagenta);
  gcorr->SetMarkerColor(kMagenta);
  gl5->SetLineColor(kRed);
  gl5->SetMarkerColor(kRed);
  
  ggen->SetMarkerStyle(kFullStar);
  ggen->SetMarkerSize(1.5);
  grec->SetMarkerStyle(kFullSquare);
  gcorr->SetMarkerStyle(kOpenCircle);
  gl5->SetMarkerStyle(kFullCircle);
  
  leg->AddEntry(gcorr,"CORR","lp");
  leg->AddEntry(gl5,  "L5",  "lp");
  leg->AddEntry(ggen, "GEN", "lp");
  leg->AddEntry(grec, "REC", "lp");

  mg->Draw("AP");
  mg->GetHistogram()->SetTitle("top mass width (RMS)");
  mg->GetHistogram()->SetYTitle("RMS(m_{top}) [GeV]");
  
  for (unsigned int i=0;i<vec_histos.size();i++) {
    int ibin=0; while ((i+1)>mg->GetHistogram()->GetBinLowEdge(ibin)) { ibin++;}
    mg->GetHistogram()->GetXaxis()
      ->SetBinLabel(ibin,alg_label(vec_histos[i]->algname).c_str());
  }
  mg->GetHistogram()->GetXaxis()->LabelsOption("u");
  mg->GetHistogram()->GetXaxis()->SetLabelOffset(0.01);
  mg->GetHistogram()->GetXaxis()->SetLabelSize(0.065);
  mg->GetHistogram()->GetYaxis()->SetTitleOffset(0.75);
  leg->Draw();

  TLatex pttxt;
  stringstream sspt; sspt<<"jet p_{T}#geq "<<vec_histos[0]->ptmin<<" GeV";
  pttxt.SetNDC();
  pttxt.DrawLatex(0.7,0.92,sspt.str().c_str());

  ctopwidth->Modified();
}


//______________________________________________________________________________
TGraphErrors* Wmass_mean(std::vector<histograms*> vec_histos,int level,
                          const std::string& channel)
{
  int     n  = (int)vec_histos.size();
  double *x  = new double[n];
  double *y  = new double[n];
  double *ex = new double[n];
  double *ey = new double[n];
  for (int i=0;i<n;i++) {
    x[i]      = i+1;
    ex[i]     = 0.0;
    if (channel.empty()) {
      y[i] =vec_histos[i]->hmw[level]->GetMean();
      ey[i]=vec_histos[i]->hmw[level]->GetMeanError();
    }
    else if (channel=="lj") {
      y[i] =vec_histos[i]->hmw_lj[level]->GetMean();
      ey[i]=vec_histos[i]->hmw_lj[level]->GetMeanError();
    }
    else if (channel=="aj") {
      y[i] =vec_histos[i]->hmw_aj[level]->GetMean();
      ey[i]=vec_histos[i]->hmw_aj[level]->GetMeanError();
    }
    else assert(0);
  }
  
  return new TGraphErrors(n,x,y,ex,ey);
}


//______________________________________________________________________________
TGraphErrors* topmass_mean(std::vector<histograms*> vec_histos,int level,
                           const std::string& channel)
{
  int     n  = (int)vec_histos.size();
  double *x  = new double[n];
  double *y  = new double[n];
  double *ex = new double[n];
  double *ey = new double[n];
  for (int i=0;i<n;i++) {
    x[i]      = i+1;
    ex[i]     = 0.0;
    if (channel.empty()) {
      y[i] =vec_histos[i]->hmt[level]->GetMean();
      ey[i]=vec_histos[i]->hmt[level]->GetMeanError();
    }
    else if (channel=="lj") {
      y[i] =vec_histos[i]->hmt_lj[level]->GetMean();
      ey[i]=vec_histos[i]->hmt_lj[level]->GetMeanError();
    }
    else if (channel=="aj") {
      y[i] =vec_histos[i]->hmt_aj[level]->GetMean();
      ey[i]=vec_histos[i]->hmt_aj[level]->GetMeanError();
    }
    else assert(0);
  }
  
  return new TGraphErrors(n,x,y,ex,ey);
}


//______________________________________________________________________________
TGraphErrors* Wmass_width(std::vector<histograms*> vec_histos,int level,
                          const std::string& channel)
{
  int     n  = (int)vec_histos.size();
  double *x  = new double[n];
  double *y  = new double[n];
  double *ex = new double[n];
  double *ey = new double[n];
  for (int i=0;i<n;i++) {
    x[i]      = i+1;
    ex[i]     = 0.0;
    if (channel.empty()) {
      y[i] =vec_histos[i]->hmw[level]->GetRMS();
      ey[i]=vec_histos[i]->hmw[level]->GetRMSError();
    }
    else if (channel=="lj") {
      y[i] =vec_histos[i]->hmw_lj[level]->GetRMS();
      ey[i]=vec_histos[i]->hmw_lj[level]->GetRMSError();
    }
    else if (channel=="aj") {
      y[i] =vec_histos[i]->hmw_aj[level]->GetRMS();
      ey[i]=vec_histos[i]->hmw_aj[level]->GetRMSError();
    }
    else assert(0);
  }
  
  return new TGraphErrors(n,x,y,ex,ey);
}


//______________________________________________________________________________
TGraphErrors* topmass_width(std::vector<histograms*> vec_histos,int level,
                            const std::string& channel)
{
  int     n  = (int)vec_histos.size();
  double *x  = new double[n];
  double *y  = new double[n];
  double *ex = new double[n];
  double *ey = new double[n];
  for (int i=0;i<n;i++) {
    x[i]      = i+1;
    ex[i]     = 0.0;
    if (channel.empty()) {
      y[i] =vec_histos[i]->hmt[level]->GetRMS();
      ey[i]=vec_histos[i]->hmt[level]->GetRMSError();
    }
    else if (channel=="lj") {
      y[i] =vec_histos[i]->hmt_lj[level]->GetRMS();
      ey[i]=vec_histos[i]->hmt_lj[level]->GetRMSError();
    }
    else if (channel=="aj") {
      y[i] =vec_histos[i]->hmt_aj[level]->GetRMS();
      ey[i]=vec_histos[i]->hmt_aj[level]->GetRMSError();
    }
    else assert(0);
  }
  
  return new TGraphErrors(n,x,y,ex,ey);
}


//______________________________________________________________________________
bool parse_filename(const std::string& filename,float& ptmin,float& drmax)
{
  unsigned int pos;
  string tmp=filename;
  while ((pos=tmp.find("/"))!=string::npos) tmp=tmp.substr(pos+1);
  pos = tmp.find("topmass_pt");
  if (pos==string::npos) {
    cout<<"ERROR parsing "<<filename<<endl;
    return false;
  }
  tmp = tmp.substr(pos+10);
  
  pos = tmp.find("_dr");
  if (pos==string::npos) {
    cout<<"ERROR parsing "<<filename<<endl;
    return false;
  }
  stringstream sspt; sspt<<tmp.substr(0,pos);tmp=tmp.substr(pos+3);
  sspt>>ptmin;
  cout<<"ptmin="<<ptmin<<endl;
  
  pos = tmp.find(".root");
  if (pos==string::npos) {
    cout<<"ERROR parsing "<<filename<<endl;
    return false;
  }
  stringstream ssdr; ssdr<<tmp.substr(0,pos);
  ssdr>>drmax;
  cout<<"drmax="<<drmax<<endl;
  
  return true;
}


//______________________________________________________________________________
string alg_label(const string& alg_name)
{
  string alg=alg_name.substr(0,2);
  string par=alg_name.substr(2,alg_name.find('.'));
  if      (alg=="kt") return "Fast k_{T}, D=0."+par;
  else if (alg=="sc") return "SISCone, R=0."+par;
  else if (alg=="mc") return "MidCone, R=0."+par;
  else if (alg=="ic") return "ItCone,  R=0."+par;
  else return "UNKNOWN";
  
}


//______________________________________________________________________________
Color_t alg_color(const string& alg_name)
{
  if (alg_name=="kt3") return kBlue-9;
  if (alg_name=="kt4") return kBlue;
  if (alg_name=="kt6") return kBlue+3;
  if (alg_name=="sc4") return kMagenta-9;
  if (alg_name=="sc5") return kMagenta;
  if (alg_name=="sc7") return kMagenta+4;
  if (alg_name=="mc4") return kCyan-9;
  if (alg_name=="mc5") return kCyan;
  if (alg_name=="mc7") return kCyan+4;
  if (alg_name=="ic4") return kGreen-9;
  if (alg_name=="ic5") return kGreen;
  if (alg_name=="ic7") return kGreen+4;
  return kBlack;
}


//______________________________________________________________________________
Style_t alg_style(const string& alg_name)
{
  if (alg_name=="kt3") return kFullTriangleDown;
  if (alg_name=="kt4") return kFullCircle;
  if (alg_name=="kt6") return kFullTriangleUp;
  if (alg_name=="sc4") return kFullTriangleDown;
  if (alg_name=="sc5") return kFullCircle;
  if (alg_name=="sc7") return kFullTriangleUp;
  if (alg_name=="mc4") return kFullTriangleDown;
  if (alg_name=="mc5") return kFullCircle;
  if (alg_name=="mc7") return kFullTriangleUp;
  if (alg_name=="ic4") return kFullTriangleDown;
  if (alg_name=="ic5") return kFullCircle;
  if (alg_name=="ic7") return kFullTriangleUp;
  return kOpenCircle;
}
Revision:	1.6
Committed:	Sun Dec 9 18:12:17 2007 UTC (17 years, 4 months ago) by schiefer
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.5:	+107 -30 lines
Log Message:	* empty log message *