Jeng/scripts/fitRelIso_Data.C

#include <vector>
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TH2.h>
#include <THStack.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TMultiGraph.h>
#include <TGraphErrors.h>
#include <TLegend.h>
#include <TPaveStats.h>
#include <TLatex.h>
#include "Math/GSLIntegrator.h"
#include "Math/WrappedTF1.h"
#include "map.h"
#include <sstream>
#include <fstream>
#include "TopStyle/tdrstyle.C"
#include "TopStyle/CMSTopStyle.cc"

using namespace std;

double sf_mc = 1.;

// Dynamic(true) or static(false) range
bool dynamic      = true;
TString isoname   = "New";
TString s0        = "fuser";//fuser,Landau,gauss
TString fitOpt    = "QLLMER0";
bool debug        = true;
double yMax       = 0.;
bool DrawOverFlow = true;
bool drawtail     = true;
double wxmin      = 1.;
double wxmax      = 1.;
double plot_xMax  = 1.4; // max of x axis for output plots; 1.0 for Old reliso
double plot_xMin  = 0.3; // min of x axis for for Old reliso
int scaleType     = 0; // 0: scale MC w.r.t. total number of events; 1: scale QCD number of event to fitted result
bool drawStats    = false;

map<TString,TCanvas*> cvs; // map of usual histogram

double round(double num, int places)
{
    double temp, mult;
    mult = pow(10.0, places);
    temp = floor(num * mult + 0.5);
    temp = temp / mult;
    return temp;
}

double stdDev(int data[], int n, double mean)
{
  double dev = 0;
  double sum2 = 0;

  for ( int i = 0; i < n; i++ ) {
    sum2 += pow((data[i] - mean),2);
  }
  dev = sqrt(sum2/(n - 1));
  return dev;
}

void cmsPrel(double intLumi){
  TLatex latex;
  latex.SetNDC();
  latex.SetTextSize(0.05);

  latex.SetTextAlign(31); // align right
  latex.DrawLatex(0.98,0.945,"#sqrt{s} = 7 TeV");
  if (intLumi > 0.) {
    latex.SetTextAlign(31); // align right
    latex.DrawLatex(0.98,0.68,Form("#int #font[12]{L}dt = %.1fnb^{-1}",intLumi));
  }
  latex.SetTextAlign(11); // align left
  latex.DrawLatex(0.02,0.945,"#font[22]{CMS preliminary 2010}");
  return;
}

void fitRelIso_Data()
{
  //setTDRStyle();
  CMSTopStyle style;
  gROOT->SetStyle("CMS");
  style.setupICHEPv1();

  //New RelIso
  if (isoname == "New") {
    //s0 = "landau";
    wxmax = 4.;//4
    wxmin = 2.;
    plot_xMin = 0.;
    // Select fitting function:
  }
  //Old RelIso
  if (isoname == "Old") {
    s0 = "gaus"; // for old
    //s0 = "fuser";
    wxmin = 2.;
    wxmax = 1.;
    DrawOverFlow = false;
    plot_xMax = 1.;
  }

  cout << " =========================== " << endl;
  cout << " = Start fitting           = " << endl;
  cout << " =========================== " << endl;

  // ////////////////////////
  // Input/output fileName
  // ////////////////////////
  TString infile1,infile2,infile3,infile4,outFile;
  if (isoname == "New") {//skimmed361p3_20100624: Loose ; skimmed361p3_20100623: tight
    infile1 = "skimmed361p3_20100720/Data_250_Incl_NewRelIso_wErrors_TandL_4";
    infile2 = "skimmed361p3_MC/MC_Mu15_Incl_NewRelIso_wErrors_TandL_4"; //InclusiveMu15
    //infile2 = "skimmed361p3_MC/MC_ppMuX_Incl_NewRelIso_wErrors_TandL_4";
    infile3 = "skimmed361p3_MC/MC_Wjets_Incl_NewRelIso_wErrors_TandL_4";
    infile4 = "skimmed361p3_MC/MC_Zjets_Incl_NewRelIso_wErrors_TandL_4";

    //outFile = "Plots/20100720/QCD_data_LandauFit_4_Incl_ppMuX_Scale0_250_125_"; //ppMuX
    outFile = "Plots/20100720/QCD_data_LandauFit_4_Incl_Scale0_250_125_"; //InclusiveMu15
  }
  if (isoname == "Old") {
    infile1 = "skimmed361p3_20100720/Data_250_Incl1_OldRelIso_wErrors_TandL_4";
    infile2 = "skimmed361p3_MC/MC_Mu15_Incl1_OldRelIso_wErrors_TandL_4"; //InclusiveMu15
    //infile2 = "skimmed361p3_MC/MC_ppMuX_Incl_OldRelIso_wErrors_TandL_4";
    infile3 = "skimmed361p3_MC/MC_Wjets_Incl1_OldRelIso_wErrors_TandL_4";
    infile4 = "skimmed361p3_MC/MC_Zjets_Incl1_OldRelIso_wErrors_TandL_4";

    outFile = "Plots/20100720/QCD_data_GaussFit_4_Incl1_Scale0_250_";
  }
  // Output file

  ofstream outTxtFile;
  outTxtFile.open(TString(outFile+"Results.txt"));
  
  double binwds[5] = {0.05,0.02,0.01,0.005,0.002};
  double binwidth[4];
  if (isoname == "New") {
    for (int i=0; i<4;++i)
      binwidth[i] = binwds[i];
  }
  else {
    for (int i=0; i<4;++i)
      binwidth[i] = binwds[i+1];
  }

  // ////////////////////////
  // Looping of files of different binning
  // ////////////////////////
  int res_nbkg[4] = {0,0,0,0};
  int res_nsig[4] = {0,0,0,0};
  int nfiles = sizeof(binwidth)/sizeof(double);

  for(int ibin = 0; ibin < nfiles ; ++ibin) {
    // ////////////////////////
    // define fit region:
    // ////////////////////////
    double ax0 = 0.0;
    double axf = 0.0;
    double bx0 = 0.0;
    double bxf = 0.0;
    double loosecut = 0.;

    cout<<">>>>> Use "<<isoname<<" RelIso"<<endl;

    if ( isoname.Contains("Old") ) {
      // old
      loosecut = 0.9; // 10/11
      ax0 = 20./21.; // 0.95
      axf = 1.0 + binwidth[ibin];
      bx0 = 0.35; // Must be smaller than x at peak //0.35
      bxf = 0.9; // Should not overlap signal region //0.9
    }
    else if ( isoname.Contains("New") ) {
      // new
      ax0 = 0.0;
      axf = 0.05; // 1./19. = old 0.95
      loosecut = 0.125; // 1./9. = old 0.9
      bx0 = 0.1;//0.2
      bxf = 1.;//1.2
    }

    //if (isoname == "Old" && ibin == 0) continue;
    TString plotFile = outFile;
    if (drawStats) gStyle->SetOptFit(100);
    cout << " =========================== " << endl;
    cout << " = Bin width = " << binwidth[ibin] << endl;
    cout << " =========================== " << endl;

    map<TString,TFile*> fs;
    map<TString,TH1*> hs; // map of usual histogram
    map<TString,TH1*> hso; // map of histo with overflow bin drawn
    map<TString,std::pair<TString,double> > m_fs;
    map<TString,THStack*> hst; // map of stacked histogram

    ostringstream suffix1,suffix2;
    suffix1 << binwidth[ibin];
    suffix2 << ibin;

    // Input files
    // Data
    m_fs["Data"] = pair<TString,double>(TString(infile1+"_"+suffix1.str()+".root"),1.);
    // MC
    m_fs["QCD"] = pair<TString,double>(TString(infile2+"_"+suffix1.str()+".root"),0.0062348); //InclusiveMu15
    //m_fs["QCD"] = pair<TString,double>(TString(infile2+"_"+suffix1.str()+".root"),2.9647348);
    m_fs["WJets"] = pair<TString,double>(TString(infile3+"_"+suffix1.str()+".root"),0.0007934);
    m_fs["ZJets"] = pair<TString,double>(TString(infile4+"_"+suffix1.str()+".root"),0.0007024);

    // User fit function
    //     TF1 *fuser = new TF1("fuser","1./(1.+[0]*TMath::Landau(x,[1],[2],0))",ax0,axf);
    //     TF1 *fu0 = new TF1("fu0","landau",0.,ax0);
    //     TF1 *fu1 = new TF1("fu1","gaus",0.,ax0);
    //     TF1 *fu2 = new TF1("fu2","landau",ax0,axf);
    //     TF1 *fu3 = new TF1("fu3","landau(3)+gaus(2)",0.,ax0);
    //     TF1 *fuser = new TF1("fuser","[0]*TMath::Landau(x,[1],[2],1)",ax0,axf);
    TF1 *fuser = new TF1("fuser","landau",ax0,axf);
    fuser->SetParameters(10.,0.3,1000.);//const ==> Remember to change according to L_{int}
    fuser->SetParLimits(2,0.15,1.);//MPV
    fuser->SetParLimits(3,0.,1.);//sigma
    //fuser->SetParLimits(3,0.,0.5);

    // File type
    TString fileName = m_fs["Data"].first;
    //cout << fileName << endl;

    Int_t jbin;
    Double_t na[4],na1[4];
    int j = 0;
    for (std::map<TString,pair<TString,double> >::iterator itFile = m_fs.begin();
         itFile != m_fs.end(); ++itFile, ++j) {
      cout << (*itFile).first << "; " << ((*itFile).second).first << endl;
      fs[(*itFile).first] = TFile::Open(((*itFile).second).first,"READ");
      hs[(*itFile).first] = (TH1D*) fs[(*itFile).first]->Get("histos/h_RelIso_all");
      TTree *tree = (TTree*) fs[(*itFile).first]->Get("myTree");
      tree->SetBranchAddress("na",&na[j]);
      if (fileName.Contains("TandL"))
        tree->SetBranchAddress("na1",&na1[j]);
      if ((*itFile).first == "Data") {
        tree->SetBranchAddress("jbin",&jbin);
        if (isoname == "New") sf_mc = hs[(*itFile).first]->Integral(1,round(0.8/binwidth[ibin],0));
        if (isoname == "Old") sf_mc = hs[(*itFile).first]->Integral(round(0.45/binwidth[ibin],0),hs[(*itFile).first]->GetMaximumBin());
        //      sf_mc = hs[(*itFile).first]->GetEntries();
        if (scaleType == 0) cout << ">>>>> Total events of data        = " << sf_mc << " <<<<<<<<<<<<<<<<" << endl;
      }
      else
        hs[(*itFile).first]->Scale(((*itFile).second).second);
      tree->GetEntry(0);
    }

    // Calculate scale factor for MC; normalized to total number of events from data
    hs["MCBKG"] = (TH1D*) hs["QCD"]->Clone();
    hs["MCBKG"]->Add((TH1D*) hs["WJets"]->Clone());
    hs["MCBKG"]->Add((TH1D*) hs["ZJets"]->Clone());
    double nTotMC = 0;
    if (isoname == "New") nTotMC = hs["MCBKG"]->Integral(1,round(0.8/binwidth[ibin],0));
    if (isoname == "Old") nTotMC = hs["MCBKG"]->Integral(round(0.45/binwidth[ibin],0),hs["MCBKG"]->GetMaximumBin());
    if (scaleType == 0) cout << ">>>>> Total weighted MC BKG       = " << nTotMC  << " <<<<<<<<<<<<<<<<" << endl;
    //double nTotMC = hs["QCD"]->GetEntries()*m_fs["QCD"].second + hs["WJets"]->GetEntries()*m_fs["WJets"].second;
    //cout << ">>>>> Total weighted MC BKG       = " << nTotMC << " <<<<<<<<<<<<<<<<" << endl;
    sf_mc /= nTotMC;
    if (scaleType == 0) cout << ">>>>> Scale factor for MC (sf_mc) = " << sf_mc << " <<<<<<<<<<<<<<<<" << endl;

    int niter = 1;
    double pass[3] = {0.,0.,0.};
    double bound[2] = {0.,0.};
    double ns_tight[3] = {0.,0.,0.};
    double ns_loose[3] = {0.,0.,0.};
    double ns2 = 0.;
    double chi2[2] = {0.,0.};
    int ndof[2] = {0,0};

    // Fit QCD in background region
    hs["data_all_clone"] = (TH1D*) hs["Data"]->Clone();
    //hs["data_all_clone"]->SetLineWidth(2);
    hs["data_all_clone"]->SetMinimum(0);
    hs["mc_qcd_clone"] = (TH1D*) hs["QCD"]->Clone();
    hs["mc_wj_clone"] = (TH1D*) hs["WJets"]->Clone();
    hs["mc_zj_clone"] = (TH1D*) hs["ZJets"]->Clone();

    // Get num of bin at peak and bin width
    Int_t nbxMax = hs["data_all_clone"]->GetXaxis()->FindBin(2.0);
    hs["data_all_clone"]->GetXaxis()->SetRange(round(0.05/binwidth[ibin],0),hs["data_all_clone"]->GetXaxis()->FindBin(bxf));
    Int_t nbMax = hs["data_all_clone"]->GetMaximumBin();
    cout << ">>>>> Bin number at peak = " << nbMax << endl;
    hs["data_all_clone"]->GetXaxis()->SetRange(1,nbxMax);
    TAxis *axis0 = hs["data_all_clone"]->GetXaxis();
    double bwidth = axis0->GetBinWidth(nbMax);
    cout << ">>>>> Histo bin width = " << bwidth << endl;

    TString title;
    if (isoname == "New") title = "Relative Isolation";
    if (isoname == "Old") title = "Relative Isolation (old)";
//     if ( jbin == -1 )
//       title += " (N_{jets} #geq 0)";
//     else if ( jbin == 5 ) {
//       title += "(#geq 1 jets (up to 3 jets))";
//     }
//     else if ( jbin > 10) {
//       ostringstream tmp_;
//       tmp_ << (jbin-10);
//       title += " (N_{jets} = ";
//       title += tmp_.str();
//       title += ")";
//     }
//     else if ( jbin <= 4 ) {
//       title += " (N_{jets} #geq ";
//       title += jbin;
//       title += ")";
//     }
//     else {
//       cout<<"jbin = " << jbin << " not supported!"<<endl;
//       return;
//     }

    cvs[TString("cv"+suffix2.str())] = new TCanvas(TString("cv"+suffix2.str()),TString("cv"+suffix2.str()),10,10,700,500);
    cvs[TString("cv"+suffix2.str())]->cd();

    // Data
    hs["data_all_clone"]->GetXaxis()->SetTitle(title);
    hs["data_all_clone"]->GetXaxis()->SetLabelFont(42);
    hs["data_all_clone"]->GetXaxis()->SetLabelSize(0.04);
    hs["data_all_clone"]->GetXaxis()->SetTitleFont(42);
    hs["data_all_clone"]->GetXaxis()->SetTitleSize(0.05);
    hs["data_all_clone"]->GetXaxis()->SetTitleOffset(1.2);

    hs["data_all_clone"]->GetYaxis()->SetTitle(TString("Events / "+suffix1.str()));
    //hs["data_all_clone"]->GetYaxis()->SetTitle(TString("Events"));
    hs["data_all_clone"]->GetYaxis()->SetLabelFont(42);
    hs["data_all_clone"]->GetYaxis()->SetLabelSize(0.04);
    hs["data_all_clone"]->GetYaxis()->SetTitleFont(42);
    hs["data_all_clone"]->GetYaxis()->SetTitleSize(0.05);
    hs["data_all_clone"]->GetYaxis()->SetTitleOffset(1.4);
    hs["data_all_clone"]->GetXaxis()->SetRangeUser(plot_xMin,plot_xMax);
    hs["data_all_clone"]->SetMarkerStyle(20);
    hs["data_all_clone"]->SetMarkerSize(0.8);
    hs["data_tmp_clone"] = (TH1D*) hs["data_all_clone"]->Clone(); // for zoomed in plot

    // Find good y range first
    if (isoname == "New") {
      if ((hs["data_all_clone"]->GetBinContent(1) + hs["data_all_clone"]->GetBinError(1))*1.05 > yMax)
        yMax = (hs["data_all_clone"]->GetBinContent(1) + hs["data_all_clone"]->GetBinError(1))*1.05;
      if ((hs["mc_qcd_clone"]->GetBinContent(1)+hs["mc_wj_clone"]->GetBinContent(1))*1.05 > yMax) 
        yMax = (hs["mc_qcd_clone"]->GetBinContent(1)+hs["mc_wj_clone"]->GetBinContent(1))*1.05;
    }
    if (isoname == "Old") {
      if ((hs["data_all_clone"]->GetBinContent(nbMax) + hs["data_all_clone"]->GetBinError(nbMax))*1.05 > yMax)
        yMax = (hs["data_all_clone"]->GetBinContent(nbMax) + hs["data_all_clone"]->GetBinError(nbMax))*1.05;
      if ((hs["mc_qcd_clone"]->GetBinContent(nbMax)+hs["mc_wj_clone"]->GetBinContent(nbMax))*1.05 > yMax) 
        yMax = (hs["mc_qcd_clone"]->GetBinContent(nbMax)+hs["mc_wj_clone"]->GetBinContent(nbMax))*1.05;
    }

    cout << " ===============> YMax = " << yMax << endl;
    //hs["data_all_clone"]->GetYaxis()->SetRangeUser(0.,yMax);
    hs["data_all_clone"]->Draw("e");

    cout << "\n =========================== " << endl;
    cout << " = Iteration step: "<< niter << endl;
    cout << " =========================== " << endl;
    cout << ">>>>> Start fitting control region: " << bx0 << " to " << bxf << endl;

    // Very first fit
    cout << ">>>>> Fitting with " << s0 << " function!" << endl;
    fuser->SetParameters(10.,0.3,100.);//const
    fuser->SetParLimits(2,0.15,1.);//MPV
    fuser->SetParLimits(3,0.,1.);//sigma
    hs["data_all_clone"]->Fit(s0,fitOpt,"",bx0,bxf);

    // Very first fit function
    TF1 *f0 = (TF1*) hs["data_all_clone"]->GetFunction(s0)->Clone();
    Int_t npar = f0->GetNpar();
    chi2[0] = f0->GetChisquare();
    ndof[0] = f0->GetNDF();
    pass[0] = chi2[0]/ndof[0]; // Very first fit norm chi2. Fixed range.
    cout << ">>> Norm chi2 = " << pass[0] << endl;
    if (!dynamic)
      pass[1] = pass[0];
    else
      pass[1] = 999.; // require at least 2 iterations

    cout << " >> Peak at " << f0->GetMaximumX(bx0,bxf) << endl;

    double * par0 = new double [npar];
    double * parErr0 = new double [npar];

    cout << " >> Parameters of initial fit: " << endl;
    for (Int_t i=0;i<npar;i++) {
      printf(" >> %s = %g +- %g\n",
             f0->GetParName(i),
             f0->GetParameter(i),
             f0->GetParError(i)
             );
      par0[i] = f0->GetParameter(i);
      parErr0[i] = f0->GetParError(i);
    }
    if(dynamic) {
      bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
      bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
      cout << "[Before] bx0 = " << bx0 << endl;
      cout << "[Before] bxf = " << bxf << endl;
    }
    if (isoname == "New") {
      bx0 = ( bx0 <= loosecut ? loosecut : bx0 );
      bxf = ( bxf >= hs["data_all_clone"]->GetXaxis()->GetXmax() ?
              hs["data_all_clone"]->GetXaxis()->GetXmax() :
              ( bxf > bx0 ? bxf : 
                (bx0 + wxmax*f0->GetParameter(2) >= 2. ? 2. : bx0 + wxmax*f0->GetParameter(2) ) ) );
    }
    if (isoname == "Old") {
      bxf = ( bxf > loosecut ? loosecut : bxf );
      bx0 = ( bx0 < 0. ? 0. : bx0 );
    }
//     if (isoname == "New") {
//       bx0 = ( bx0 <= loosecut ? loosecut : 
//            ( bx0 < hs["data_all_clone"]->GetXaxis()->GetXmax() ? bx0 : 
//              ( hs["data_all_clone"]->GetXaxis()->GetXmax() + loosecut )/2. ) );
//       bxf = ( bxf >= hs["data_all_clone"]->GetXaxis()->GetXmax() ?
//            hs["data_all_clone"]->GetXaxis()->GetXmax() :
//            ( bxf > bx0 ? bxf :
//              ( bx0 + hs["data_all_clone"]->GetXaxis()->GetXmax() )/2. ) );
//     }
//     if (isoname == "Old") {
//       bxf = ( bxf > loosecut ? loosecut :
//            ( bxf > hs["data_all_clone"]->GetXaxis()->GetXmin() ? bxf :
//              ( hs["data_all_clone"]->GetXaxis()->GetXmin() + loosecut)/2. ) );
//       bx0 = ( bx0 < hs["data_all_clone"]->GetXaxis()->GetXmin() ? hs["data_all_clone"]->GetXaxis()->GetXmin() :
//            ( bx0 > bxf ? ( bxf + hs["data_all_clone"]->GetXaxis()->GetXmin() )/2. : bx0) );
//     }
    bound[0] = bx0;
    bound[1] = bxf;
    if (dynamic) {
      cout << ">>> Range for next iteration (bx0, bxf) = (" << bx0;
      cout << ", " << bxf << ")" << endl;
    }

    double delta = pass[1]-pass[2];

    while (dynamic && niter <= 20 && abs(delta) > 0.000001 && ndof[0] >= 4) {
      if (niter > 1)
        pass[1] = pass[2];
      if (delta != 0){
        niter++;
        cout << "\n =========================== " << endl;
        cout << " = Iteration step: "<< niter << endl;
        cout << " =========================== " << endl;
        cout << ">>>>> Previous step norm. chi2 = " << pass[1];
        if (niter == 2)
          cout << "; pass[0] = " << pass[0];

        cout << "\n>>> Start fitting new range (bx0, bxf) = (" << bx0
             << ", " << bxf << ")" << endl;

        fuser->SetParameters(10.,0.3,100.);//const
        fuser->SetParLimits(2,0.15,1.);//MPV
        fuser->SetParLimits(3,0.,1.);//sigma

        hs["data_all_clone"]->Fit(s0,fitOpt,"",bx0,bxf);
        TF1 *fi = (TF1*) hs["data_all_clone"]->GetFunction(s0)->Clone();
        cout << " >> fi->GetChisquare() = " << fi->GetChisquare() << "; fi->GetNDF() = " << fi->GetNDF() << endl;
        if ( fi->GetNDF() > 0 )
          pass[2] = fi->GetChisquare()/fi->GetNDF();
        else {
          chi2[1] = 999.;
          ndof[1] = 1.;
          pass[2] = 999.;
        }
        cout << " >> Current step norm. chi2  = " << pass[2] << endl;
        delta = pass[1]-pass[2];
        cout << " >> Delta = " << delta << endl;

        if (delta > 0 && bx0 >= loosecut && bxf <= hs["data_all_clone"]->GetXaxis()->GetXmax() && bx0 < bxf) {
          chi2[1] = fi->GetChisquare();
          ndof[1] = fi->GetNDF();
          bound[0] = bx0;
          bound[1] = bxf;
          printf(" >> Function has %i parameters. Chisquare = %g\n",
                 npar,
                 fi->GetChisquare());
          for (Int_t i=0;i<npar;i++) {
            printf(" >> %s = %g +- %g\n",
                   fi->GetParName(i),
                   fi->GetParameter(i),
                   fi->GetParError(i)
                   );
            par0[i] = fi->GetParameter(i);
            parErr0[i] = fi->GetParError(i);
            //       cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
          }
          bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2);
          bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2);
          cout << "[Before] bx0 = " << bx0 << endl;
          cout << "[Before] bxf = " << bxf << endl;

          if (isoname == "New") {
            bx0 = ( bx0 < loosecut ? loosecut : bx0 );
            bxf = ( bxf > 2. ? 2. : ( bxf > bx0 ? bxf : 1.1*bx0 ) );
          }
          if (isoname == "Old") {
            bxf = ( bxf > loosecut ? loosecut : bxf );
            bxf = ( bx0 < 0. ? 0. : bx0 );
          }
//        if (isoname == "New") {
//          bx0 = ( bx0 <= loosecut ? loosecut : 
//                  ( bx0 < hs["data_all_clone"]->GetXaxis()->GetXmax() ? bx0 : 
//                    ( hs["data_all_clone"]->GetXaxis()->GetXmax() + loosecut )/2. ) );
//          bxf = ( bxf >= hs["data_all_clone"]->GetXaxis()->GetXmax() ?
//                  hs["data_all_clone"]->GetXaxis()->GetXmax() :
//                  ( bxf > bx0 ? bxf :
//                    ( bx0 + hs["data_all_clone"]->GetXaxis()->GetXmax() )/2. ) );
//        }
//        if (isoname == "Old") {
//          bxf = ( bxf > loosecut ? loosecut :
//                  ( bxf > hs["data_all_clone"]->GetXaxis()->GetXmin() ? bxf :
//                    ( hs["data_all_clone"]->GetXaxis()->GetXmin() + loosecut)/2. ) );
//          bx0 = ( bx0 < hs["data_all_clone"]->GetXaxis()->GetXmin() ? hs["data_all_clone"]->GetXaxis()->GetXmin() :
//                  ( bx0 > bxf ? ( bxf + hs["data_all_clone"]->GetXaxis()->GetXmin() )/2. :bx0 ) );
//        }
          cout << ">>> Range for next iteration (bx0, bxf) = (" << bx0;
          cout << ", " << bxf << ")" << endl;
        }
        else if (delta < 0) {
          cout << ">>> Use previous fit parameters for extrapolation fit" << endl;
          delta = 0;
        }
        else { // delta == 0
          cout << ">> Fit converges." << endl;
          cout << ">> Best fitting region = (" << bound[0];
          cout << ", " << bound[1] << ")" << endl;
        }
        if (niter == 2)
          pass[0]=pass[2]; // First dynamic fit norm. chi2
      } // delta != 0
    } // End while

    // Output equivalent range
    if (dynamic) {
      if (bound[0]/binwidth[ibin] - floor(bound[0]/binwidth[ibin]) > 0.5)
        bound[0] = round(bound[0]/binwidth[ibin],0) * binwidth[ibin];
      else
        bound[0] = (floor(bound[0]/binwidth[ibin])+0.5) * binwidth[ibin];
      if (bound[1]/binwidth[ibin] - floor(bound[1]/binwidth[ibin]) >= 0.5)
        bound[1] = round(bound[1]/binwidth[ibin],0) * binwidth[ibin];
      else
        bound[1] = round(bound[1]/binwidth[ibin],0) * binwidth[ibin];
    }

    if(debug){
      cout << ">>> First dynamic fit norm. chi2   = " << pass[0] << endl;
      cout << ">>> Min dynamic norm. chi2         = " << pass[1] << endl;
      cout << ">>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
      cout << ">>> Best fitting region (bx0, bxf) = (" << bound[0] << ", ";
      cout << bound[1] << ")"<< endl;
    }

    // Get number of events within fitted region
    TAxis *axis = hs["data_all_clone"]->GetXaxis();
    int bmin = axis->FindBin(bound[0]);
    int bmax = axis->FindBin(bound[1]);
    double nfac1 = hs["data_all_clone"]->Integral(bmin,bmax);
    cout << ">>> Initial Nfac = " << nfac1 << endl;
    nfac1 -= (hs["data_all_clone"]->GetBinContent(bmin))*(bound[0]-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
    nfac1 -= (hs["data_all_clone"]->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bound[1])/axis->GetBinWidth(bmax);
    cout << ">>> Final Nfac = " << nfac1 << endl;

    TF1 *f1 = (TF1*) f0->Clone();
    f1->SetRange(bound[0],bound[1]);

    cout << ">>> Parameters used for " << s0 << " function: " << endl;
    for (int i=0;i<npar;i++) {
      cout << " >> Par[" << i << "] = " << par0[i] << "; Error = " << parErr0[i] << endl;
      f1->SetParameter(i,par0[i]);
    }
    f1->SetLineColor(2);

    TF1 *f2 = (TF1*) f1->Clone();
    f2->SetRange(ax0,axf);
    f2->SetLineColor(4);
    //f2->SetLineWidth(3.5);

    // ////////////////////////
    // Extract number of qcd events
    // ////////////////////////  
    // Method 1: Ratio
    int np = 100;
    double *x=new double[np];
    double *w=new double[np];
    f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);

    double nfac0 = 0.;
    if (fileName.Contains("0053") || 
        fileName.Contains("TandL") ||
        fileName.Contains("Tight")) {
      ns_tight[0] = f2->IntegralFast(np,x,w,ax0,axf);
      nfac0 = f2->Integral(bound[0],bound[1]);
      cout << " >> Area of signal region = " << ns_tight[0] << endl;
      cout << " >> Area of control region = " << nfac0 << endl;
      ns_tight[0]/=nfac0;
      ns_tight[0]*=nfac1;
    }
    else if (fileName.Contains("011") || 
             fileName.Contains("Loose")) {
      ns_tight[0] = f2->IntegralFast(np,x,w,ax0,loosecut);
      nfac0 = f2->Integral(bound[0],bound[1]);
      ns_tight[0]/=nfac0;
      ns_tight[0]*=nfac1;
    }

    ns_loose[0] = f2->IntegralFast(np,x,w,ax0,loosecut);
    ns_loose[0]/=(f2->Integral(bound[0],bound[1]));
    ns_loose[0]*=nfac1;

    cout << ">>> Tight Ns by usual integral    = " << ns_tight[0] << endl;
    if (fileName.Contains("TandL") || fileName.Contains("011"))
      cout << ">>> Loose Ns by usual integral    = " << ns_loose[0] << endl;

    delete [] x;
    delete [] w;

    // Another way to do integration
    TF1 *g = (TF1*)f1->Clone();
    ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000); 
    ROOT::Math::WrappedTF1 wf(*g);
    ig.SetFunction(wf);
    if (fileName.Contains("0053") || 
        fileName.Contains("TandL") || 
        fileName.Contains("Tight"))
      ns2 = ig.Integral(ax0,axf);
    else if (fileName.Contains("011") ||
             fileName.Contains("Loose"))
      ns2 = ig.Integral(ax0,loosecut);

    ns2/=(g->Integral(bound[0],bound[1]));
    ns2*=nfac1;
    cout<<">>> Ns by MathMore integral       = "<<ns2<<endl;

    // Method 2: area of trapezoid
    ns_tight[1] = (TMath::Landau(ax0,par0[1],par0[2],kFALSE)*par0[0] + TMath::Landau(axf,par0[1],par0[2],kFALSE)*par0[0])
      * fabs(axf-ax0)/(2*bwidth);
    ns_loose[1] = (TMath::Landau(ax0,par0[1],par0[2],kFALSE)*par0[0] + TMath::Landau(loosecut,par0[1],par0[2],kFALSE)*par0[0])
      * fabs(loosecut-ax0)/(2*bwidth);
    // Method 3: area of Landau function (default)
    ns_tight[2] = f2->Integral(ax0,axf)/bwidth;
    ns_loose[2] = f2->Integral(ax0,loosecut)/bwidth;

    // ////////////////////////
    // Histos and extrapolation
    // ////////////////////////

    // Draw stacked MC
    double sf_mc1 = 0.;
    if (isoname == "New") sf_mc1 = ns_tight[2]/hs["QCD"]->Integral(1,round(0.05/binwidth[ibin],0));
    if (isoname == "Old") sf_mc1 = ns_tight[2]/hs["QCD"]->Integral(round(0.95/binwidth[ibin],0)+1,hs["QCD"]->GetMaximumBin());
    if (scaleType == 0) {
      hs["mc_qcd_clone"]->Scale(sf_mc);
      hs["mc_wj_clone"]->Scale(sf_mc);
    }
    if (scaleType == 1) {
      sf_mc = sf_mc1;
      hs["mc_qcd_clone"]->Scale(sf_mc1);
      hs["mc_wj_clone"]->Scale(sf_mc1);
      cout << ">>>>> Scale factor for MC (sf_mc) = " << sf_mc1 << " <<<<<<<<<<<<<<<<" << endl;
    }

    if (DrawOverFlow) { // store histo with overflow bin drawn
      int nbins_ = hs["data_all_clone"]->GetXaxis()->FindBin(plot_xMax);
      for (std::map<TString,TH1*>::iterator itHS = hs.begin();
           itHS != hs.end(); ++itHS) {
        hso[(*itHS).first] = (TH1D*) ((*itHS).second)->Clone();
        hso[(*itHS).first]->SetBins(nbins_,0.,plot_xMax+bwidth);
        double v_ovf = 0.;
        for (int i=0; i < ((*itHS).second)->GetNbinsX()+1 ; i++) {
          if (i < nbins_) {
            hso[(*itHS).first]->SetBinContent(i+1,hs[(*itHS).first]->GetBinContent(i+1));
            hso[(*itHS).first]->SetBinError(i+1,hs[(*itHS).first]->GetBinError(i+1));
          }
          else {
            v_ovf += hs[(*itHS).first]->GetBinContent(i+1);
          }
        }
        hso[(*itHS).first]->SetBinContent(nbins_,v_ovf);
        hso[(*itHS).first]->SetBinError(nbins_,TMath::Sqrt(v_ovf));
      }
    }

    hst[TString("Data"+suffix2.str())] = new THStack(TString("Data"+suffix2.str()),TString("Data"+suffix2.str()+"stacked"));
    if (DrawOverFlow) {
      hso["mc_qcd_clone"]->SetStats(0);
      //hso["mc_qcd_clone"]->SetLineWidth(2);
      //hso["mc_qcd_clone"]->SetLineColor(style.QCDColor);
      hso["mc_qcd_clone"]->SetFillColor(style.QCDColor);
      hso["mc_qcd_clone"]->SetFillStyle(style.QCDFill);

      hso["mc_zj_clone"]->SetStats(0);
      hso["mc_zj_clone"]->SetFillStyle(style.DYZJetsFill);
      hso["mc_zj_clone"]->SetFillColor(style.DYZJetsColor);

      hso["mc_wj_clone"]->SetStats(0);
      //hso["mc_wj_clone"]->SetLineWidth(2);
      //hso["mc_wj_clone"]->SetLineColor(style.WJetsColor);
      hso["mc_wj_clone"]->SetFillStyle(style.WJetsFill);
      hso["mc_wj_clone"]->SetFillColor(style.WJetsColor);
      hst[TString("Data"+suffix2.str())]->Add(hso["mc_qcd_clone"]);
      hst[TString("Data"+suffix2.str())]->Add(hso["mc_zj_clone"]);
      hst[TString("Data"+suffix2.str())]->Add(hso["mc_wj_clone"]);
    }
    else {
      hs["mc_qcd_clone"]->SetStats(0);
      //hs["mc_qcd_clone"]->SetLineWidth(2);
      //hs["mc_qcd_clone"]->SetLineColor(style.QCDColor);
      hs["mc_qcd_clone"]->SetFillStyle(style.QCDFill);
      hs["mc_qcd_clone"]->SetFillColor(style.QCDColor);

      hs["mc_zj_clone"]->SetStats(0);
      hs["mc_zj_clone"]->SetFillStyle(style.DYZJetsFill);
      hs["mc_zj_clone"]->SetFillColor(style.DYZJetsColor);

      hs["mc_wj_clone"]->SetStats(0);
      //hs["mc_wj_clone"]->SetLineWidth(2);
      //hs["mc_wj_clone"]->SetLineColor(style.WJetsColor);
      hs["mc_wj_clone"]->SetFillStyle(style.WJetsFill);
      hs["mc_wj_clone"]->SetFillColor(style.WJetsColor);

      hst[TString("Data"+suffix2.str())]->Add(hs["mc_qcd_clone"]);
      hst[TString("Data"+suffix2.str())]->Add(hs["mc_zj_clone"]);
      hst[TString("Data"+suffix2.str())]->Add(hs["mc_wj_clone"]);
    }

    // Draw data
    if (DrawOverFlow) {
      hso["data_all_clone"]->GetYaxis()->SetRangeUser(0.,yMax);
      hso["data_all_clone"]->Draw("esame");
    }

    hst[TString("Data"+suffix2.str())]->Draw("sames hist");

    // Draw data again
    gStyle->SetErrorX(0.5);

    if (DrawOverFlow) {
      hso["data_all_clone"]->Draw("esame");
    }
    else {
      hs["data_all_clone"]->DrawClone("esame"); // Draw on top of others
    }

    // Draw fit function

    // fit function in control region
    //f1->SetLineWidth(3.5);
    f1->Draw("same");

    // fit fuction in signal region
    f2->Draw("same");

    // Connect fit function between regions and tail
    TF1 *fin = (TF1*) f1->Clone();
    if (isoname == "New") {
      fin->SetRange(axf,bound[0]);
      fin->SetLineStyle(2);
      fin->SetLineColor(9);
      fin->SetLineWidth(3.5);
      fin->Draw("same");

      if (drawtail){
        TF1 *fine = (TF1*) f1->Clone();
        fine->SetRange(bound[1],2.0);
        fine->SetLineStyle(2);
        fine->SetLineColor(9);
        fine->SetLineWidth(3.5);
        fine->Draw("same");
      }
    }
    if (isoname == "Old") {
      fin->SetRange(bound[1],ax0);
      fin->SetLineStyle(2);
      fin->SetLineColor(9);
      fin->SetLineWidth(3.5);
      fin->Draw("same");

      if (drawtail){
        TF1 *fine = (TF1*) f1->Clone();
        fine->SetRange(0.,bound[0]);
        fine->SetLineStyle(2);
        fine->SetLineColor(9);
        fine->SetLineWidth(3.5);
        fine->Draw("same");
      }
    }

    TString mytxt;
    if (isoname == "New") mytxt = "Landau Fit";
    if (isoname == "Old") mytxt = "Gaussian Fit";

    if (drawStats) {
      // Draw Stats box
      TPaveStats * tps1;
      if(!dynamic){
        cvs[TString("cv"+suffix2.str())]->Update();
        tps1 = (TPaveStats*) cvs[TString("cv"+suffix2.str())]->GetPrimitive("stats");
      }
      else {
        hs["data_all_clone"]->Fit(s0,fitOpt,"",bound[0],bound[1]);
        cvs[TString("cv"+suffix2.str())]->Update();
        tps1 = (TPaveStats*) cvs[TString("cv"+suffix2.str())]->GetPrimitive("stats");
      }

      tps1->SetName(TString(s0+" fit"+suffix2.str()));

      TList *list = tps1->GetListOfLines();
      //   TText *tconst = tps1->GetLineWith("Constant");
      //   list->Remove(tconst);
      TLatex *myt = new TLatex(0,0,mytxt);
      //TLatex *myt = new TLatex(0,0,"#font[22]{Landau Fit}"); 
      list->AddFirst(myt);
      hs["data_all_clone"]->SetStats(0);
      gStyle->SetOptFit(0);
      cvs[TString("cv"+suffix2.str())]->Modified();

      tps1->SetTextFont(42);
      tps1->SetTextSize(0.03);
      tps1->SetFillColor(0);
      tps1->SetFillStyle(0);
      tps1->SetBorderSize(0);
      tps1->SetX1NDC(0.63);
      tps1->SetX2NDC(0.88);
      tps1->SetY1NDC(0.74);
      tps1->SetY2NDC(0.87);
      cvs[TString("cv"+suffix2.str())]->cd();
      gPad->Update();
    }

    TLatex *text1, *text2;
    text1 = new TLatex(3.570061, 23.08044, "CMS Preliminary");
    text1->SetNDC();
    text1->SetTextAlign(13);
    if (isoname == "New")
      text1->SetX(0.24);
    else
      text1->SetX(0.54);
    text1->SetY(0.92);
    //text1->SetLineWidth(2);
    text1->SetTextFont(42);
    text1->SetTextSize(0.05);
    text1->SetTextSizePixels(24);// dflt=28
    text1->Draw();

    text2 = new TLatex(13.570061, 23.08044, "0.25 pb^{-1} at #sqrt{s} = 7 TeV");
    text2->SetNDC();
    text2->SetTextAlign(13);
    if (isoname == "New")
      text2->SetX(0.24);
    else 
      text2->SetX(0.54);
    text2->SetY(0.85);
    //text2->SetLineWidth(2);
    text2->SetTextFont(42);
    text2->SetTextSize(0.05);
    text2->SetTextSizePixels(24);// dflt=28
    text2->Draw();

    //cmsPrel(35);
    // Label histo
    //TLegend * leg1 = new TLegend(0.58,0.45,0.85,0.8);
    TLegend * leg1;
    if (isoname == "New")
      leg1 = new TLegend(0.66,0.52,0.95,0.92);
    else
      leg1 = new TLegend(0.28,0.57,0.55,0.92);
    //leg1->SetHeader("#font[22]{CMS preliminary 2010   #sqrt{s} = 7 TeV }");
    leg1->SetTextSize(0.03);
    leg1->SetFillColor(0);
    leg1->SetFillStyle(0);
    if (DrawOverFlow) {
      //leg1->AddEntry(hso["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
      leg1->AddEntry(hso["data_all_clone"],"Data","lp");
//       leg1->AddEntry(f1,TString(mytxt+" in control region"),"l");
//       leg1->AddEntry(f2,"Extrapolation to signal region","l");
      leg1->AddEntry(f1,mytxt,"l");
      leg1->AddEntry(f2,"Extrapolation","l");
      leg1->AddEntry(hso["mc_wj_clone"],style.WJetsText,"f");
      leg1->AddEntry(hso["mc_zj_clone"],style.DYZJetsText,"f");
      leg1->AddEntry(hso["mc_qcd_clone"],style.QCDText,"f");
    }
    else {
      //leg1->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
      leg1->AddEntry(hs["data_all_clone"],"Data","lp");
      leg1->AddEntry(f1,mytxt,"l");
      leg1->AddEntry(f2,"Extrapolation","l");
      leg1->AddEntry(hs["mc_wj_clone"],style.WJetsText,"f");
      leg1->AddEntry(hs["mc_zj_clone"],style.DYZJetsText,"f");
      leg1->AddEntry(hs["mc_qcd_clone"],style.QCDText,"f");
    }
    leg1->Draw();

    plotFile = plotFile + suffix1.str() + ".pdf";
//     if (ibin == 1 || ibin == 2)
    cvs[TString("cv"+suffix2.str())]->Print(plotFile);

    // Print results
    cout << "\n";
    cout << " =========================== " << endl;
    cout << " = Results                 = " << endl;
    cout << " =========================== " << endl;

    if ( jbin == -1 )
      cout << ">>>>> Inclusive jet bin " << ":" << endl;
    else 
      cout << ">>>>> Jet bin " << jbin << ":" << endl;
    cout << ">>>>> Bin number at peak = " << nbMax << endl;
    cout << ">>>>> Bin width = " << bwidth << endl;
    if (dynamic) {
      cout << ">>>>> First dynamic fit norm. chi2   = " << pass[0] << endl;
      cout << ">>>>> Min dynamic fit norm. chi2     = " << pass[1] << endl;
      cout << ">>>>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
      cout << ">>>>> Best fitting region (bx0, bxf) = (" << bound[0] << ", " << bound[1] << ")"<< endl;
    }
    cout << ">>>>> Calculated Ns (ratio)     = " << ns_tight[0] << endl;
    cout << ">>>>> Calculated Ns (trapezoid) = " << ns_tight[1] << endl;
    cout << ">>>>> Calculated Ns (integral)  = " << ns_tight[2] << endl;
    cout << ">>>>> Total events expected Na  = " << na[0] << endl;
    cout << ">>>>> Num of Signal events      = " << na[0] - round(ns_tight[2],0) << endl;

    res_nbkg[ibin] = round(ns_tight[2],0);
    res_nsig[ibin] = na[0]-round(ns_tight[2],0);

    //   if (fileName.Contains("TandL")) {
    //     cout << "====================================" << endl;
    //     cout << ">>>>> Calculated Loose Ns (ratio)     = " << ns_loose[0] << endl;
    //     cout << ">>>>> Calculated Loose Ns (trapezoid) = " << ns_loose[1] << endl;
    //     cout << ">>>>> Calculated Loose Ns (integral)  = " << ns_loose[2] << endl;
    //     cout << ">>>>> Total expected events Loose Na  = " << na1[0] << endl;
    //     cout << ">>>>> Num of Signal events            = " << na1[0] - round(ns_loose[2],0) << endl;
    //   }
    //   cout << "Landau(0.05) = " << TMath::Landau(0.05,par0[1],par0[2],kFALSE)*par0[0] << endl;
    //   cout << "Landau(0.1) = " << TMath::Landau(0.1,par0[1],par0[2],kFALSE)*par0[0] << endl;
    cout << " ===========================\n\n " << endl;

    outTxtFile << "===========================" << endl;
    outTxtFile << ">>>>> Bin width                  = " << bwidth << endl;
    //outTxtFile << ">>>>> Calculated Ns (ratio)      = " << ns_tight[0] << endl;
    //outTxtFile << ">>>>> Calculated Ns (trapezoid)  = " << ns_tight[1] << endl;
    if (dynamic) {
      outTxtFile << ">>>>> Min dynamic norm. chi2     = " << pass[1] << endl;
      outTxtFile << ">>>>> chi2/ndof                  = " << chi2[1] << "/" << ndof[1] << " = "
                 << chi2[1]/ndof[1] << endl;
    }
    else {
      outTxtFile << ">>>>> Min norm. chi2             = " << pass[1] << endl;
      outTxtFile << ">>>>> chi2/ndof                  = " << chi2[0] << "/" << ndof[0] << " = " << chi2[0]/ndof[0] << endl;
    }
    outTxtFile << ">>>>> Best fitting region        = (" << bound[0];
    outTxtFile << " - " << bound[1] << ")" << endl;
    outTxtFile << ">>>>> Total number of events     = " << na[0] << endl;
    outTxtFile << ">>>>> Calculated bkg (integral)  = " << round(ns_tight[2],0) << endl;
    outTxtFile << ">>>>> Num of Signal events       = " << na[0] - round(ns_tight[2],0) << endl;
    outTxtFile << ">>>>> Total MC QCD expected(N)   = " << round(na[1]*((m_fs["QCD"]).second)*sf_mc,1) << " +/- " <<
      round(TMath::Sqrt(na[1])*((m_fs["QCD"]).second)*sf_mc,1) << endl;
    outTxtFile << ">>>>> Total MC WJets expected(N) = " << round(na[2]*((m_fs["WJets"]).second)*sf_mc,1) << " +/- " <<
      round(TMath::Sqrt(na[2])*((m_fs["WJets"]).second)*sf_mc,1) << endl;
    outTxtFile << ">>>>> Total MC ZJets expected(N) = " << round(na[3]*((m_fs["ZJets"]).second)*sf_mc,1) << " +/- " <<
      round(TMath::Sqrt(na[3])*((m_fs["ZJets"]).second)*sf_mc,1) << endl;
    outTxtFile << ">>>>> Total MC QCD expected(R)   = " << round(na[1]*((m_fs["QCD"]).second),1) << " +/- " <<
      round(TMath::Sqrt(na[1])*((m_fs["QCD"]).second),1)<< endl;
    outTxtFile << ">>>>> Total MC WJets expected(R) = " << round(na[2]*((m_fs["WJets"]).second),1) << " +/- " <<
      round(TMath::Sqrt(na[2])*((m_fs["WJets"]).second),1) << endl;
    outTxtFile << ">>>>> Total MC ZJets expected(R) = " << round(na[3]*((m_fs["ZJets"]).second),1) << " +/- " <<
      round(TMath::Sqrt(na[3])*((m_fs["ZJets"]).second),2) << endl;
    outTxtFile << "===========================\n\n " << endl;


    // /////////////
    // Zoom-in plot
    // /////////////
    cvs[TString("cv1"+suffix2.str())] = new TCanvas(TString("cv1"+suffix2.str()),TString("cv1"+suffix2.str()),10,10,700,500);
    cvs[TString("cv1"+suffix2.str())]->cd();

    if (isoname == "New") {
      hs["data_tmp_clone"]->GetXaxis()->SetRangeUser(0.,0.06);
      hs["data_tmp_clone"]->GetYaxis()->SetRangeUser(0.3,10000.*binwidth[ibin]/binwidth[3]);
    }
    if (isoname == "Old") {
      hs["data_tmp_clone"]->GetXaxis()->SetRangeUser(0.9,axf);
      hs["data_tmp_clone"]->GetYaxis()->SetRangeUser(0.3,10000.*binwidth[ibin]/binwidth[3]);
    }
    hs["data_tmp_clone"]->Draw("e");
    hst[TString("Data"+suffix2.str())]->Draw("sames hist");
    hs["data_tmp_clone"]->Draw("esame");

    text1->Draw();
    text2->Draw();

    TLegend * leg2;
    if (isoname == "New")
      leg2 = new TLegend(0.66,0.52,0.95,0.92);
    else
      leg2 = new TLegend(0.28,0.57,0.55,0.92);
    //leg2->SetHeader("#font[22]{CMS preliminary 2010   #sqrt{s} = 7 TeV }");
    leg2->SetTextSize(0.03);
    leg2->SetFillColor(0);
    leg2->SetFillStyle(0);
    //leg2->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
    if (DrawOverFlow) {
      //leg2->AddEntry(hso["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
      leg2->AddEntry(hso["data_all_clone"],"Data","lp");
      //leg2->AddEntry(f2,TString("Extrapolation of "+mytxt),"l");
      leg2->AddEntry(f2,"Extrapolation","l");
      leg2->AddEntry(hso["mc_wj_clone"],style.WJetsText,"f");
      leg2->AddEntry(hso["mc_zj_clone"],style.DYZJetsText,"f");
      leg2->AddEntry(hso["mc_qcd_clone"],style.QCDText,"f");
    }
    else {
      //leg2->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
      leg2->AddEntry(hs["data_all_clone"],"Data","lp");
      //leg2->AddEntry(f2,TString("Extrapolation of "+mytxt),"l");
      leg2->AddEntry(f2,"Extrapolation","l");
      leg2->AddEntry(hs["mc_wj_clone"],style.WJetsText,"f");
      leg2->AddEntry(hs["mc_zj_clone"],style.DYZJetsText,"f");
      leg2->AddEntry(hs["mc_qcd_clone"],style.QCDText,"f");
    }
    leg2->Draw();

    //   TF1 *fb2 = new TF1("fa3","[0]*TMath::Landau(x,[1],[2],0)",ax0,axf);
    //   fb2->SetParameters(par0[0],par0[1],par0[2]);
    //   fb2->SetLineColor(5);
    //   fb2->Draw("same");
    f2->Draw("same");
    fin->Draw("same");
    //hs["Data"]->GetXaxis()->SetRangeUser(0.,0.06);
    //hs["Data"]->GetYaxis()->SetRangeUser(0.3,300.);
    hs["Data"]->SetLineWidth(2);
    cvs[TString("cv1"+suffix2.str())]->SetLogy();
    hs["Data"]->DrawClone("esame"); // Draw on top of others

    plotFile = plotFile(0,plotFile.Length()-4) + "_zoomed.pdf";
    //if (ibin == 2 || ibin == 3)
    cvs[TString("cv1"+suffix2.str())]->Print(plotFile);

    hs.clear();
    hso.clear();
    fs.clear();
    m_fs.clear();
  }
  // Calculate errors:

  float sum[2] = {0.,0.};
  float average[2] = {0.,0.};
  float errors = 0.;

  if ( nfiles > 1) {
    for (int ii=0; ii<nfiles; ii++) {
      sum[0] += res_nbkg[ii];
      sum[1] += res_nsig[ii];
    }
    average[0] = round(sum[0] / nfiles,0);
    average[1] = round(sum[1] / nfiles,0);
    errors = stdDev(res_nbkg,nfiles,average[0]);

    cout << "Mean bkg estimated = " << average[0] << endl;
    cout << "Mean sig estimated = " << average[1] << endl;
    cout << "Standard deviation = " << errors << endl;
    cout << "Uncertainties      = " << errors/average[0]*100. << "% " << endl;

    outTxtFile << "Estimated bkg      = ";
    for (int jj=0; jj<nfiles; ++jj)
      outTxtFile << res_nbkg[jj] << "  ";
    outTxtFile << endl;
    outTxtFile << "Mean bkg estimated = " << average[0] << endl;
    outTxtFile << "Mean sig estimated = " << average[1] << endl;
    outTxtFile << "Standard deviation = " << errors << endl;
    outTxtFile << "Uncertainties      = " << errors/average[0]*100. << "% " << endl;
  }
  outTxtFile.close();
}
Revision:	1.3
Committed:	Wed Jul 21 00:55:58 2010 UTC (14 years, 9 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.2:	+552 -305 lines
Log Message:	update