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>

using namespace std;

double sf_mc = 1.;

// Dynamic(true) or static(false) range
bool dynamic      = true;
TString isoname   = "New";
TString s0        = "fuser";
TString fitOpt    = "QLM0";
bool debug        = true;
double yMax       = 0.;
bool DrawOverFlow = true;
bool drawtail     = false;
double wxmin      = 1.;
double wxmax      = 1.;
double plot_xMax  = 1.4; // max of x axis for output plots; 1.0 for Old reliso

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()
{
  //New RelIso
  if (isoname == "New") {
    s0 = "landau";
    wxmin = 10.;
    wxmax = wxmin/3.;
    // Select fitting function:
  }
  //Old RelIso
  if (isoname == "Old") {
    s0 = "gaus"; // for old
    wxmin = 1.;
    wxmax = 1./3.;
    DrawOverFlow = false;
    plot_xMax = 1.;
  }

  gROOT->SetStyle("CMS");

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

  // Input fileName
  TString infile1 = "skimmed361p3_20100623/Data_Incl_NewRelIso_wErrors_TandL_3";
  TString infile2 = "skimmed361p3_20100623/MC_ppMuX_Incl_NewRelIso_wErrors_TandL_3";
  TString infile3 = "skimmed361p3_20100623/MC_Wjets_Incl_NewRelIso_wErrors_TandL_3";
//   TString infile1 = "skimmed361p3_20100623/Data_Incl_OldRelIso_wErrors_TandL_3";
//   TString infile2 = "skimmed361p3_20100623/MC_ppMuX_Incl_OldRelIso_wErrors_TandL_3";
//   TString infile3 = "skimmed361p3_20100623/MC_Wjets_Incl_OldRelIso_wErrors_TandL_3";

  // Output file
  TString outFile = "Plots/20100623/QCD_data_LandauFit_3_";
//   TString outFile = "Plots/20100623/QCD_data_GaussFit_3_";
  ofstream outTxtFile;
  outTxtFile.open(TString(outFile+"Results.txt"));

  double binwidth[5] = {0.05,0.02,0.01,0.005,0.002};
  int res_nbkg[5] = {0,0,0,0,0};
  int res_nsig[5] = {0,0,0,0,0};
  int nfiles = sizeof(binwidth)/sizeof(double);
  for(int ibin = 0; ibin < nfiles ; ++ibin) {
    TString plotFile = outFile;
    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.154166);
    m_fs["WJets"] = pair<TString,double>(TString(infile3+"_"+suffix1.str()+".root"),0.0000369);

    // User fit function
    TF1 *fuser = new TF1("fuser","[0]*TMath::Landau(x,[1],[2],0)",0.,2.);
    fuser->SetParameters(1.,0.5,1.);
    fuser->SetParLimits(2,0.125,1.);
    fuser->SetParLimits(3,0.,0.5);

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

    Int_t jbin;
    Double_t na[3],na1[3];
    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);
        sf_mc = hs[(*itFile).first]->Integral();
        //      sf_mc = hs[(*itFile).first]->GetEntries();
        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());
    double nTotMC = hs["MCBKG"]->Integral();
    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;
    cout << ">>>>> Scale factor for MC (sf_mc) = " << sf_mc << " <<<<<<<<<<<<<<<<" << endl;

    hs["QCD"]->Scale(sf_mc);
    hs["WJets"]->Scale(sf_mc);

    // 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 = 10./11.; // 0.9
      ax0 = 20./21.; // 0.95
      axf = 1.0 + binwidth[ibin];
      bx0 = 0.3; // Must be smaller than x at peak
      bxf = 0.9; // Should not overlap signal region
    }
    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.2;
      bxf = 1.2;
    }

    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.;

    // Fit QCD in background region
    cout<<">>>>> Fitting with "<<s0<<" function!"<<endl;
    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();

    // 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(2,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 = "RelIso";
    if ( jbin == -1 )
      title += " (inclusive)";
    else if ( jbin == 4 )
      title += " (#geq 4 jets)";
    else if ( jbin == 1 ) {
      title += " (";
      title += jbin;
      title += " jet)";
    }
    else if ( jbin < 4 ) {
      title += " (";
      title += jbin;
      title += " jets)";
      hs["data_all_clone"]->SetTitle(title);
    }
    else if ( jbin == 5 ) {
      hs["data_all_clone"]->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
    }
    else {
      cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
      return;
    }

    cvs[TString("cv"+suffix2.str())] = new TCanvas(TString("cv"+suffix2.str()),TString("cv"+suffix2.str()),600,600);
    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.03);
    hs["data_all_clone"]->GetXaxis()->SetTitleFont(42);
    hs["data_all_clone"]->GetXaxis()->SetTitleSize(0.035);
    hs["data_all_clone"]->GetXaxis()->SetTitleOffset(1.2);

    hs["data_all_clone"]->GetYaxis()->SetTitle("Number of Events");
    hs["data_all_clone"]->GetYaxis()->SetLabelFont(42);
    hs["data_all_clone"]->GetYaxis()->SetLabelSize(0.03);
    hs["data_all_clone"]->GetYaxis()->SetTitleFont(42);
    hs["data_all_clone"]->GetYaxis()->SetTitleSize(0.035);
    hs["data_all_clone"]->GetYaxis()->SetTitleOffset(1.4);
    hs["data_all_clone"]->GetXaxis()->SetRangeUser(0.,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

    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(40);
      hso["mc_qcd_clone"]->SetFillStyle(3018);
      hso["mc_qcd_clone"]->SetFillColor(38);

      hso["mc_wj_clone"]->SetStats(0);
      hso["mc_wj_clone"]->SetLineWidth(2);
      hso["mc_wj_clone"]->SetLineColor(8);
      //hso["mc_wj_clone"]->SetFillStyle(3001);
      hso["mc_wj_clone"]->SetFillColor(3);
      hst[TString("Data"+suffix2.str())]->Add(hso["mc_qcd_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(40);
      hs["mc_qcd_clone"]->SetFillStyle(3018);
      hs["mc_qcd_clone"]->SetFillColor(38);

      hs["mc_wj_clone"]->SetStats(0);
      hs["mc_wj_clone"]->SetLineWidth(2);
      hs["mc_wj_clone"]->SetLineColor(8);
      //hs["mc_wj_clone"]->SetFillStyle(3001);
      hs["mc_wj_clone"]->SetFillColor(3);

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

    // 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
    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();
    Double_t chi2 = f0->GetChisquare();
    Int_t ndof = f0->GetNDF();
    pass[0] = chi2/ndof; // 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);
    }
    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 );
    }
    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 >= 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;

        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 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) {
          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);
          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 );
          }
          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
    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;

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

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

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

    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);

    // Draw on top of everything
    gStyle->SetErrorX(0.5);
    if (DrawOverFlow) {
      hso["data_all_clone"]->Draw("esame");
    }
    else {
      hs["data_all_clone"]->Draw("esame"); // Draw on top of others
    }
    //f1->SetLineWidth(3.5);
    f1->Draw("same");

    // Connect fitted and extrapolated region
    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],axf);
      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");
      }
    }

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

    // ////////////////////////
    // 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
    ns_tight[2] = f2->Integral(ax0,axf)/bwidth;
    ns_loose[2] = f2->Integral(ax0,loosecut)/bwidth;

    // 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("Landau fit"+suffix2.str()));


    TList *list = tps1->GetListOfLines();
    //   TText *tconst = tps1->GetLineWith("Constant");
    //   list->Remove(tconst);
    TLatex *myt = new TLatex(0,0,"Landau Fit"); 
    //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();

    //cmsPrel(13);
    // Label histo
    TLegend * leg1 = new TLegend(0.18,0.47,0.45,0.87);
    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 = 13 nb^{-1}","p");
      leg1->AddEntry(f1,"Fit of all events in control region","l");
      leg1->AddEntry(f2,"Extrapolation to signal region","l");
      leg1->AddEntry(hso["mc_qcd_clone"],"QCD","f");
      leg1->AddEntry(hso["mc_wj_clone"],"W+Jets","f");
    }
    else {
      leg1->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 13 nb^{-1}","p");
      leg1->AddEntry(f1,"Fit of all events in control region","l");
      leg1->AddEntry(f2,"Extrapolation to signal region","l");
      leg1->AddEntry(hs["mc_qcd_clone"],"QCD","f");
      leg1->AddEntry(hs["mc_wj_clone"],"W+Jets","f");
    }
    leg1->Draw();

    plotFile = plotFile + suffix1.str() + ".pdf";
    //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;
    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     = " << round(na[1]*((m_fs["QCD"]).second)*sf_mc,0) << endl;
    outTxtFile << ">>>>> Total MC WJets expected   = " << round(na[2]*((m_fs["WJets"]).second)*sf_mc,0) << endl;
    outTxtFile << "===========================\n\n " << endl;


    cvs[TString("cv1"+suffix2.str())] = new TCanvas(TString("cv1"+suffix2.str()),TString("cv1"+suffix2.str()),600,600);
    cvs[TString("cv1"+suffix2.str())]->cd();
    if (isoname == "New")
      hs["data_tmp_clone"]->GetXaxis()->SetRangeUser(0.,0.1);
    if (isoname == "Old")
      hs["data_tmp_clone"]->GetXaxis()->SetRangeUser(0.9,axf);
    hs["data_tmp_clone"]->GetYaxis()->SetRangeUser(0.,yMax);
    hs["data_tmp_clone"]->Draw("e");
    hst[TString("Data"+suffix2.str())]->Draw("sames hist");
    hs["data_tmp_clone"]->Draw("esame");

    TLegend * leg2 = new TLegend(0.38,0.6,0.65,0.85);
    leg2->SetHeader("#font[22]{CMS preliminary 2010   #sqrt{s} = 7 TeV }");
    leg2->SetTextSize(0.03);
    leg2->SetFillColor(0);
    leg2->SetFillStyle(0);
    if (DrawOverFlow) {
      leg2->AddEntry(hso["data_all_clone"],"Data   #int#font[12]{L}dt = 13 nb^{-1}","p");
      leg2->AddEntry(hso["mc_qcd_clone"],"QCD","f");
      leg2->AddEntry(hso["mc_wj_clone"],"W+Jets","f");
    }
    else { 
      leg2->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 13 nb^{-1}","p");
      leg2->AddEntry(hs["mc_qcd_clone"],"QCD","f");
      leg2->AddEntry(hs["mc_wj_clone"],"W+Jets","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");

    plotFile = plotFile(0,plotFile.Length()-4) + "_zoomed.pdf";
    //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.;

  int nmax = sizeof(res_nbkg)/sizeof(int);
  for (int ii=0; ii<nmax; ii++) {
    sum[0] += res_nbkg[ii];
    sum[1] += res_nsig[ii];
  }
  average[0] = round(sum[0] / nmax,0);
  average[1] = round(sum[1] / nmax,0);
  errors = stdDev(res_nbkg,nmax,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<nmax;++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.1
Committed:	Thu Jun 24 02:16:25 2010 UTC (14 years, 10 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
Log Message:	Add file