/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/TopQuarkAnalysis/TopTools/src/LRHelpFunctions.cc

Go to the documentation of this file.

#include "TopQuarkAnalysis/TopTools/interface/LRHelpFunctions.h"
#include "TopQuarkAnalysis/TopTools/test/tdrstyle.C"

// constructors
LRHelpFunctions::LRHelpFunctions() {
  constructPurity = false;
  setTDRStyle();
  gStyle->SetCanvasDefW(900);
}


LRHelpFunctions::LRHelpFunctions(std::vector<int> obsNr, int nrBins, std::vector<double> obsMin, std::vector<double> obsMax, std::vector<const char*> functions) { 
  obsNumbers = obsNr;
  constructPurity = false;
  setTDRStyle();
  gStyle->SetCanvasDefW(900);
  for(size_t o=0; o<obsNr.size(); o++){
    // create Signal, background and s/(s+b) histogram
    TString htS  = "Obs"; htS  += obsNr[o]; htS += "_S"; 
    TString htB  = "Obs"; htB  += obsNr[o]; htB += "_B"; 
    TString htSB = "Obs"; htSB += obsNr[o]; htSB += "_SoverSplusB"; 
    hObsS.push_back( new TH1F(htS,htS,nrBins,obsMin[o],obsMax[o]) );
    hObsB.push_back( new TH1F(htB,htB,nrBins,obsMin[o],obsMax[o]) );
    hObsSoverSplusB.push_back( new TH1F(htSB,htSB,nrBins,obsMin[o],obsMax[o]) );
    
    // create the correlation 2D plots among the observables (for signal)
    for (size_t o2 = o+1; o2 < obsNr.size(); o2++) {
      TString hcorr  = "Corr_Obs"; hcorr  += obsNr[o]; hcorr += "_Obs";  hcorr  += obsNr[o2]; 
      hObsCorr.push_back( new TH2F(hcorr,hcorr,nrBins,obsMin[o],obsMax[o],nrBins,obsMin[o2],obsMax[o2]) );
    }    
    
    // create fit functions
    TString ftSB = "F_Obs"; ftSB += obsNr[o]; ftSB += "_SoverSplusB"; 
    fObsSoverSplusB.push_back( new TF1(ftSB,functions[o],hObsS[o]->GetXaxis()->GetXmin(),hObsS[o]->GetXaxis()->GetXmax()) );
  }
}

void LRHelpFunctions::recreateFitFct(std::vector<int> obsNr, std::vector<const char*> functions) {
  if (!fObsSoverSplusB.empty()) fObsSoverSplusB.clear();
  for(size_t o=0; o<obsNr.size(); o++){
    // create fit functions
    TString ftSB = "F_Obs"; ftSB += obsNr[o]; ftSB += "_SoverSplusB";
    fObsSoverSplusB.push_back( new TF1(ftSB,functions[o],hObsS[o]->GetXaxis()->GetXmin(),hObsS[o]->GetXaxis()->GetXmax()) );
  }
}

LRHelpFunctions::LRHelpFunctions(int nrLRbins, double LRmin, double LRmax, const char* LRfunction) { 
  constructPurity = true;
  setTDRStyle();
  gStyle->SetCanvasDefW(900);

  // create LR histograms
  hLRtotS = new TH1F("hLRtotS","hLRtotS",nrLRbins,LRmin,LRmax);
  hLRtotS->GetXaxis()->SetTitle("Combined LR");
  hLRtotB = new TH1F("hLRtotB","hLRtotB",nrLRbins,LRmin,LRmax);
  hLRtotB->GetXaxis()->SetTitle("Combined LR");
  hLRtotSoverSplusB = new TH1F("hLRtotSoverSplusB","hLRtotSoverSplusB",nrLRbins,LRmin,LRmax);
    
  // create LR fit function
  fLRtotSoverSplusB = new TF1("fLRtotSoverSplusB",LRfunction,LRmin,LRmax);
}


// destructor
LRHelpFunctions::~LRHelpFunctions() {}


// member function to initialize the LR hists and fits
void LRHelpFunctions::initLRHistsAndFits(int nrLRbins, double LRmin, double LRmax, const char* LRfunction){ 
  constructPurity = true;
  // create LR histograms
  hLRtotS = new TH1F("hLRtotS","hLRtotS",nrLRbins,LRmin,LRmax);
  hLRtotB = new TH1F("hLRtotB","hLRtotB",nrLRbins,LRmin,LRmax);
  hLRtotSoverSplusB = new TH1F("hLRtotSoverSplusB","hLRtotSoverSplusB",nrLRbins,LRmin,LRmax);    
  // create LR fit function
  fLRtotSoverSplusB = new TF1("fLRtotSoverSplusB",LRfunction,LRmin,LRmax);
}



// member function to set initial values to the observable fit function
void LRHelpFunctions::setObsFitParameters(int obs,std::vector<double> fitPars){
  for(size_t fit=0; fit<fObsSoverSplusB.size(); fit++){
    TString fn = "_Obs"; fn += obs;
    if(((TString)fObsSoverSplusB[fit]->GetName()).Contains(fn)){
      for(size_t p=0; p<fitPars.size(); p++){
        fObsSoverSplusB[fit]->SetParameter(p,fitPars[p]);
      }
    }
  }
}





// member function to add observable values to the signal histograms
void LRHelpFunctions::fillToSignalHists(std::vector<double> obsVals, double weight){
  int hIndex = 0;
  for(size_t o=0; o<obsVals.size(); o++) {
    hObsS[o]->Fill(obsVals[o], weight);
    for(size_t o2=o+1; o2<obsVals.size(); o2++) {
      hObsCorr[hIndex] -> Fill(obsVals[o],obsVals[o2], weight);
      ++hIndex;
    }
  }
}

// member function to add observable values to the signal histograms
void LRHelpFunctions::fillToSignalHists(int obsNbr, double obsVal, double weight){
  TString obs = "Obs"; obs += obsNbr; obs += "_";
  for(size_t f = 0; f<hObsS.size(); f++){
    if(((TString)(hObsS[f]->GetName())).Contains(obs)) {
      hObsS[f]->Fill(obsVal, weight);
      return;
    }
  }
}

// member function to add observable values to the signal histograms
void LRHelpFunctions::fillToSignalCorrelation(int obsNbr1, double obsVal1, int obsNbr2,
        double obsVal2, double weight){
  TString hcorr  = "Corr_Obs"; hcorr  += std::min(obsNbr1,obsNbr2) ; hcorr += "_Obs";  hcorr  += std::max(obsNbr1, obsNbr2);
  for(size_t i=0; i<hObsCorr.size(); i++) {
    if(((TString)(hObsCorr[i]->GetName())) == hcorr) {
      if (obsNbr1 < obsNbr2) {
        hObsCorr[i] -> Fill(obsVal1,obsVal2, weight);
      } else {
        hObsCorr[i] -> Fill(obsVal2,obsVal1, weight);
      }
      return;
    }
  }
}


// member function to add observable values to the background histograms
void LRHelpFunctions::fillToBackgroundHists(std::vector<double> obsVals, double weight){
  for(size_t o=0; o<obsVals.size(); o++) hObsB[o]->Fill(obsVals[o], weight);
}

// member function to add observable values to the signal histograms
void LRHelpFunctions::fillToBackgroundHists(int obsNbr, double obsVal, double weight){
  TString obs = "Obs"; obs += obsNbr; obs += "_";
  for(size_t f = 0; f<hObsB.size(); f++){
    if(((TString)(hObsB[f]->GetName())).Contains(obs)) {
      hObsB[f]->Fill(obsVal, weight);
      return;
    }
  }
}

// member function to normalize the S and B spectra
void LRHelpFunctions::normalizeSandBhists(){
  for(size_t o=0; o<hObsS.size(); o++){
    // count entries in each histo. Do it this way instead of GetEntries method,
    // since the latter does not account for the weights.
    double nrSignEntries = 0, nrBackEntries = 0;
    for (int i=0; i<=hObsS[o]->GetNbinsX()+1; i++){
      nrSignEntries += hObsS[o]->GetBinContent(i);
      nrBackEntries += hObsB[o]->GetBinContent(i);
     }
    for(int b=0; b <= hObsS[o]->GetNbinsX()+1; b++){
      hObsS[o]->SetBinContent(b,hObsS[o]->GetBinContent(b)/(nrSignEntries));
      hObsB[o]->SetBinContent(b,hObsB[o]->GetBinContent(b)/(nrBackEntries));
      hObsS[o]->SetBinError(b,hObsS[o]->GetBinError(b)/(nrSignEntries));
      hObsB[o]->SetBinError(b,hObsB[o]->GetBinError(b)/(nrBackEntries));
    }
    //std::cout<<"Integral for obs"<<o<<" S: "<<hObsS[o]->Integral(0,10000)<<" & obs"<<o<<" B: "<<hObsB[o]->Integral(0,10000)<<std::endl;
  }
}

 
    
// member function to produce and fit the S/(S+B) histograms
void LRHelpFunctions::makeAndFitSoverSplusBHists(){
  for(size_t o=0; o<hObsS.size(); o++){
    for(int b=0; b <= hObsS[o]->GetNbinsX()+1; b++){
      if ((hObsS[o]->GetBinContent(b)+ hObsB[o]->GetBinContent(b)) > 0) {
        hObsSoverSplusB[o]->SetBinContent(b, hObsS[o]->GetBinContent(b) / (hObsS[o]->GetBinContent(b) + hObsB[o]->GetBinContent(b)));
        double error = sqrt( pow(hObsS[o]->GetBinError(b) * hObsB[o]->GetBinContent(b)/pow(hObsS[o]->GetBinContent(b) + hObsB[o]->GetBinContent(b),2),2)
                           + pow(hObsB[o]->GetBinError(b) * hObsS[o]->GetBinContent(b)/pow(hObsS[o]->GetBinContent(b) + hObsB[o]->GetBinContent(b),2),2) );
        hObsSoverSplusB[o]->SetBinError(b,error);
      }
    }
    hObsSoverSplusB[o]->Fit(fObsSoverSplusB[o]->GetName(),"RQ");
  }
}


// member function to read the observable hists & fits from a root-file
void LRHelpFunctions::readObsHistsAndFits(TString fileName, std::vector<int> observables, bool readLRplots){
  hObsS.clear();
  hObsB.clear();
  hObsSoverSplusB.clear();
  fObsSoverSplusB.clear();
  TFile *fitFile = new TFile(fileName, "READ");
  if(observables[0] == -1){  
    std::cout<<" ... will read hists and fit for all available observables in file "<<fileName<<std::endl;
    TList *list = fitFile->GetListOfKeys();
    TIter next(list);
    TKey *el;
    while ((el = (TKey*)next())) {
      TString keyName = el->GetName();
      if(keyName.Contains("F_") && keyName.Contains("_SoverSplusB")){
        fObsSoverSplusB.push_back( new TF1(*((TF1*) el -> ReadObj())));
      }
      else if(keyName.Contains("_SoverSplusB")){
        hObsSoverSplusB.push_back( new TH1F(*((TH1F*) el -> ReadObj())));
      }
      else if(keyName.Contains("_S")){
        hObsS.push_back( new TH1F(*((TH1F*) el -> ReadObj())));
      }
      else if(keyName.Contains("_B")){
        hObsB.push_back( new TH1F(*((TH1F*) el -> ReadObj())));
      }
      else if(keyName.Contains("Corr")){
        hObsCorr.push_back( new TH2F(*((TH2F*) el -> ReadObj())));
      }
    }
  }
  else{
    obsNumbers = observables;
    for(unsigned int obs = 0; obs < observables.size(); obs++){
      std::cout<<"  ... will read hists and fit for obs "<<observables[obs]<<" from file "<<fileName<<std::endl;
      TString hStitle  = "Obs";   hStitle  += observables[obs];   hStitle += "_S";
      hObsS.push_back( new TH1F(*((TH1F*)fitFile->GetKey(hStitle)->ReadObj())));
      TString hBtitle  = "Obs";   hBtitle  += observables[obs];   hBtitle += "_B";
      hObsB.push_back( new TH1F(*((TH1F*)fitFile->GetKey(hBtitle)->ReadObj())));
      TString hSBtitle = "Obs";   hSBtitle += observables[obs];   hSBtitle += "_SoverSplusB";
      TString fSBtitle = "F_";  fSBtitle += hSBtitle;
      hObsSoverSplusB.push_back( new TH1F(*((TH1F*)fitFile->GetKey(hSBtitle)->ReadObj())));
      fObsSoverSplusB.push_back( new TF1(*((TF1*)fitFile->GetKey(fSBtitle)->ReadObj())));
      for(unsigned int obs2 = obs+1; obs2 < observables.size(); obs2++){
        TString hCorrtitle  = "Corr_Obs"; hCorrtitle  += observables[obs];   
                hCorrtitle += "_Obs";     hCorrtitle  += observables[obs2]; 
        hObsCorr.push_back( new TH2F(*((TH2F*)fitFile->GetKey(hCorrtitle)->ReadObj())));
      }
    }
  }
  
  if(readLRplots){
    constructPurity = true;
    std::cout<<"  ... will LR s and B histos from file "<<fileName<<std::endl;
    hLRtotS = new TH1F(*((TH1F*)fitFile->GetKey("hLRtotS")->ReadObj()));
    hLRtotB = new TH1F(*((TH1F*)fitFile->GetKey("hLRtotB")->ReadObj()));
    hLRtotSoverSplusB = new TH1F(*((TH1F*)fitFile->GetKey("hLRtotSoverSplusB")->ReadObj()));
    fLRtotSoverSplusB = new TF1(*((TF1*)fitFile->GetKey("fLRtotSoverSplusB") -> ReadObj()));
  }
}




// member function to store all observable plots and fits to a ROOT-file
void  LRHelpFunctions::storeToROOTfile(TString fname){
  TFile fOut(fname,"RECREATE");
  fOut.cd();
  for(size_t o=0; o<hObsS.size(); o++){
    hObsS[o]            -> Write();
    hObsB[o]            -> Write();
    hObsSoverSplusB[o]  -> Write();
    fObsSoverSplusB[o]  -> Write();
  }
  int hIndex = 0;
  for(size_t o=0; o<hObsS.size(); o++) {
    for(size_t o2=o+1; o2<hObsS.size(); o2++) {
      hObsCorr[hIndex] -> Write();
      ++ hIndex;
    }
  }
  if(constructPurity){
    hLRtotS             -> Write();
    hLRtotB             -> Write();
    hLRtotSoverSplusB   -> Write();
    fLRtotSoverSplusB   -> Write();
    hEffvsPur           -> Write();
    hLRValvsPur         -> Write();
    hLRValvsEff         -> Write();
  }
  fOut.cd();
  fOut.Write();
  fOut.Close();
}





// member function to make some simple control plots and store them in a ps-file
void  LRHelpFunctions::storeControlPlots(TString fname){  
  TCanvas c("dummy","",1);
  c.Print(fname + "[","landscape");
  for(size_t o=0; o<hObsS.size(); o++) {
     TCanvas c2("c2","",1);
     c2.Divide(2,1);
     c2.cd(1);
     hObsS[o] -> SetLineColor(2);
     if(hObsB[o]->GetMaximum() > hObsS[o]->GetMaximum()){
       hObsB[o] -> Draw("hist");
       hObsS[o] -> Draw("histsame");
     }
     else
     {
       hObsS[o] -> Draw("hist");
       hObsB[o] -> Draw("histsame");
     }
     c2.cd(2);
     hObsSoverSplusB[o] -> Draw();
     fObsSoverSplusB[o] -> Draw("same");
     c2.Print(fname,"Landscape");
  }
  
  int hIndex = 0;
  for(size_t o=0; o<hObsS.size(); o++) {
    for(size_t o2=o+1; o2<hObsS.size(); o2++) {
      TCanvas cc("cc","",1);
      hObsCorr[hIndex] -> Draw("box");
      TPaveText pt(0.5,0.87,0.98,0.93,"blNDC");
      pt.SetFillStyle(1);
      pt.SetFillColor(0);
      pt.SetBorderSize(0);
      TString tcorr = "Corr. of "; tcorr += (int)(100.*hObsCorr[hIndex] -> GetCorrelationFactor()); tcorr += " %";
      //TText *text = pt.AddText(tcorr);
      pt.Draw("same");
      ++hIndex;
      cc.Print(fname,"Landscape");
    }
  }
  
  if(constructPurity){
    TCanvas c3("c3","",1);
    c3.Divide(2,1);
    c3.cd(1);
    hLRtotB -> Draw();
    hLRtotS -> SetLineColor(2);
    hLRtotS -> Draw("same");
    c3.cd(2);
    hLRtotSoverSplusB -> Draw();
    c3.Print(fname,"Landscape");
  
    TCanvas c4("c4","",1);
    hEffvsPur -> Draw("AL*");
    c4.Print(fname,"Landscape");

    TCanvas c5("c5","",1);
    hLRValvsPur -> Draw("AL*");
    c5.Print(fname,"Landscape");

    TCanvas c6("c6","",1);
    hLRValvsEff -> Draw("AL*");
    c6.Print(fname,"Landscape");
  } 
  c.Print(fname + "]","landscape");
}
   



// member function to calculate the LR value, using the S/N definition
double  LRHelpFunctions::calcLRval(std::vector<double> vals){
  double logLR = 0.;
  for(size_t o=0; o<fObsSoverSplusB.size(); o++){
    double SoverSplusN = fObsSoverSplusB[o]->Eval(vals[o]);
    double SoverN = 1./((1./SoverSplusN) - 1.);
    logLR += log(SoverN);
  }  
  return logLR;
}


// member function to calculate the LR value, using the definition that was
// used in the P-TDR: S/(S+N)
 
double LRHelpFunctions::calcPtdrLRval(std::vector<double> vals, bool useCorrelation) {
  double logLR = 1.;
  for(size_t o=0; o<fObsSoverSplusB.size(); o++){
    double SoverSplusN = fObsSoverSplusB[o]->Eval(vals[o]);
    if (SoverSplusN<0.0001) SoverSplusN=0.0001;
    if (useCorrelation){
      double corr = 0;
      for(size_t j=0; j<fObsSoverSplusB.size(); j++){
        if (o==j) ++corr;
        else {
          TString hcorr  = "Corr_Obs"; hcorr  += std::min(obsNumbers[o],obsNumbers[j]) ; hcorr += "_Obs";  hcorr  += std::max(obsNumbers[o],obsNumbers[j]);
          for(size_t i=0; i<hObsCorr.size(); i++) {
            if(((TString)(hObsCorr[i]->GetName())) == hcorr)
              corr += fabs(hObsCorr[i]->GetCorrelationFactor());
          }
        }
      }
     logLR *= pow(SoverSplusN/fObsSoverSplusB[o]->GetMaximum(), 1./corr);
    }
    else logLR *= SoverSplusN/fObsSoverSplusB[o]->GetMaximum();
  }
  //std::cout <<logLR<<std::endl;
  return logLR;
}


// member function to fill a signal contribution to the LR histogram
void  LRHelpFunctions::fillLRSignalHist(double val, double weight){
  //  std::cout<< "@@@===> LRHelpf Signal LRval = "<< val<<std::endl;
  hLRtotS -> Fill(val, weight);
}



// member function to fill a background contribution to the LR histogram
void  LRHelpFunctions::fillLRBackgroundHist(double val, double weight){
  // std::cout<< "@@@===> LRHelpf Backgroud LRval = "<< val<<std::endl;
  hLRtotB -> Fill(val, weight);
}
  
  
   
// member function to make and fit the purity vs. LRval histogram, and the purity vs. efficiency graph
void  LRHelpFunctions::makeAndFitPurityHists(){

  for(int b=0; b<=hLRtotS->GetNbinsX(); b++) {
    float Sint = hLRtotS->Integral(b, hLRtotS->GetNbinsX()+1);
    float Bint = hLRtotB->Integral(b, hLRtotB->GetNbinsX()+1);
    if (Sint + Bint > 0) {
      hLRtotSoverSplusB->SetBinContent(b,1. * Sint / (Sint + Bint));
      hLRtotSoverSplusB->SetBinError(b,sqrt((pow(Sint*sqrt(Bint),2)+pow(Bint*sqrt(Sint),2)))/pow((Sint+Bint),2));
    }
  }

  hLRtotS->GetXaxis()->SetTitle("Combined LR ratio");
  hLRtotB->GetXaxis()->SetTitle("Combined LR ratio");
  hLRtotSoverSplusB->GetXaxis()->SetTitle("Cut on Combined LR");
  hLRtotSoverSplusB->GetYaxis()->SetTitle("Purity");

  hLRtotSoverSplusB -> Fit(fLRtotSoverSplusB->GetName(),"RQ");  
  double totSignal = hLRtotS->Integral(0,hLRtotS->GetNbinsX()+1);
  double Eff[200], Pur[200], LRVal[200];
  if (hLRtotS->GetNbinsX()>200) {
    std::cout << "Number of bins of LR histograms can not execeed 200!"<<std::endl;
    return;
  }
  for(int cut=0; (cut<=hLRtotS->GetNbinsX())&&(cut<200) ; cut ++){
        double LRcutVal = hLRtotS->GetBinLowEdge(cut);
        Eff[cut]   = hLRtotS->Integral(cut,hLRtotS->GetNbinsX()+1)/totSignal;
        Pur[cut]   = fLRtotSoverSplusB->Eval(LRcutVal);
        LRVal[cut] = LRcutVal;
  }
  hEffvsPur = new TGraph(hLRtotS->GetNbinsX(),Eff,Pur);
  hEffvsPur -> SetName("hEffvsPur");
  hEffvsPur -> SetTitle("");
  hEffvsPur -> GetXaxis() -> SetTitle((TString)("Efficiency of cut on log combined LR"));
  hEffvsPur -> GetYaxis() -> SetTitle((TString)("Purity"));
  hEffvsPur -> GetYaxis() -> SetRangeUser(0,1.1);

  hLRValvsPur = new TGraph(hLRtotS->GetNbinsX(),LRVal,Pur);
  hLRValvsPur -> SetName("hLRValvsPur");
  hLRValvsPur -> SetTitle("");
  hLRValvsPur -> GetXaxis() -> SetTitle((TString)("Cut on the log combined LR value"));
  hLRValvsPur -> GetYaxis() -> SetTitle((TString)("Purity"));
  hLRValvsPur -> GetYaxis() -> SetRangeUser(0,1.1);

  hLRValvsEff = new TGraph(hLRtotS->GetNbinsX(),LRVal,Eff);
  hLRValvsEff -> SetName("hLRValvsEff");
  hLRValvsEff -> SetTitle("");
  hLRValvsEff -> GetXaxis() -> SetTitle((TString)("Cut on the log combined LR value"));
  hLRValvsEff -> GetYaxis() -> SetTitle((TString)("Efficiency of cut on log combined LR"));
  hLRValvsEff -> GetYaxis() -> SetRangeUser(0,1.1);

}   




// member function to calculate the probability for signal given a logLR value
double  LRHelpFunctions::calcProb(double logLR){
  return fLRtotSoverSplusB -> Eval(logLR);
}


// member to check if a certain S/S+B fit function is read from the root-file
bool    LRHelpFunctions::obsFitIncluded(int o){
  bool included = false;
  TString obs = "_Obs"; obs += o; obs += "_";
  for(size_t f = 0; f<fObsSoverSplusB.size(); f++){    
    if(((TString)(fObsSoverSplusB[f]->GetName())).Contains(obs)) included = true;
  }
  return included;
}

void LRHelpFunctions::setXlabels(const std::vector<std::string> & xLabels)
{
  if (hObsS.size() != xLabels.size()) {
    std::cout << "LRHelpFunctions::setXlabels: Number of labels ("
         << xLabels.size()<< ") does not match number of obervables("
         << hObsS.size()<<").\n";
    return;
  }
  for(size_t i = 0; i<hObsS.size(); ++i){
    hObsS[i]->GetXaxis()->SetTitle(TString(xLabels[i]));
    hObsB[i]->GetXaxis()->SetTitle(TString(xLabels[i]));
    hObsSoverSplusB[i]->GetXaxis()->SetTitle(TString(xLabels[i]));
    hObsSoverSplusB[i]->GetYaxis()->SetTitle(TString("S/(S+B)"));
  }
}

void LRHelpFunctions::setYlabels(const std::vector<std::string> & yLabels)
{
  if (hObsS.size() != yLabels.size()) {
    std::cout << "LRHelpFunctions::setYlabels: Number of labels ("
         << yLabels.size()<< ") does not match number of obervables("
         << hObsS.size()<<").\n";
    return;
  }
  for(size_t i = 0; i<hObsS.size(); ++i){
    hObsS[i]->GetYaxis()->SetTitle(TString(yLabels[i]));
    hObsB[i]->GetYaxis()->SetTitle(TString(yLabels[i]));
  }
}

void  LRHelpFunctions::singlePlot(TString fname, int obsNbr, TString extension) {
  if (!obsFitIncluded(obsNbr)) return;

  TStyle *tdrStyle = gROOT->GetStyle("tdrStyle");
  tdrStyle->SetOptFit(0);
  tdrStyle->SetOptStat(0);
  tdrStyle->SetOptTitle(0);
//   tdrStyle->SetPadTopMargin(0.01);
//   tdrStyle->SetPadBottomMargin(0.01);
//   tdrStyle->SetPadLeftMargin(0.01);
//   tdrStyle->SetPadRightMargin(0.01);

  TCanvas c2("c2","",600,300);
  c2.SetTopMargin(0.01);
  c2.SetBottomMargin(0.01);
  c2.SetLeftMargin(0.01);
  c2.SetRightMargin(0.01);
  std::cout <<fname<<std::endl;
  c2.Divide(2,1);

  TString obs = "Obs"; obs += obsNbr; obs += "_";
  for(size_t o = 0; o<hObsB.size(); ++o){
    if(((TString)(hObsB[o]->GetName())).Contains(obs)) {
      c2.cd(1);
      hObsS[o] -> SetLineColor(2);
      if(hObsB[o]->GetMaximum() > hObsS[o]->GetMaximum()){
        hObsB[o] -> Draw("hist");
        hObsS[o] -> Draw("histsame");
      }
      else
      {
        hObsS[o] -> Draw("hist");
        hObsB[o] -> Draw("histsame");
      }
      c2.cd(2);

      hObsSoverSplusB[o] -> Draw();
      fObsSoverSplusB[o] -> Draw("same");
    }
  }
  std::cout << fname+"."+extension<<std::endl;
  c2.Print(fname+"."+extension);
}

void  LRHelpFunctions::purityPlot(TString fname, TString extension)
{
  TStyle *tdrStyle = gROOT->GetStyle("tdrStyle");
  tdrStyle->SetOptFit(0);
  tdrStyle->SetOptStat(0);
  tdrStyle->SetOptTitle(0);
//   tdrStyle->SetPadTopMargin(0.01);
//   tdrStyle->SetPadBottomMargin(0.01);
//   tdrStyle->SetPadLeftMargin(0.01);
//   tdrStyle->SetPadRightMargin(0.01);

  TCanvas c2("c2","",600,300);
  c2.SetTopMargin(0.01);
  c2.SetBottomMargin(0.01);
  c2.SetLeftMargin(0.01);
  c2.SetRightMargin(0.01);
  std::cout <<fname<<std::endl;
  c2.Divide(2,1);

  c2.cd(1);
  hLRValvsPur->Draw("AP");
  c2.cd(2);
  hEffvsPur->Draw("AP");
  c2.Print(fname+"Purity."+extension);

  hLRtotS->GetXaxis()->SetNdivisions(505);
  hLRtotB->GetXaxis()->SetNdivisions(505);
  TCanvas c3("c2","",300,300);
//   c3.SetTopMargin(0.01);
//   c3.SetBottomMargin(0.01);
//   c3.SetLeftMargin(0.01);
//   c3.SetRightMargin(0.01);
  hLRtotS->GetXaxis()->SetTitle("Combined LR");
  hLRtotB->GetXaxis()->SetTitle("Combined LR");

    hLRtotB -> Draw();
    hLRtotS -> SetLineColor(2);
    hLRtotS -> Draw("same");
  c3.Print(fname+"Dist."+extension);


    TCanvas c5("c3","",900,600);
    c5.Divide(2,2);
    c5.cd(1);
    hLRtotB -> Draw();
    hLRtotS -> SetLineColor(2);
    hLRtotS -> Draw("same");
    c5.cd(3);
    hLRValvsEff -> Draw("AP");
    c5.cd(2);
    hEffvsPur -> Draw("AP");
    c5.cd(4);
    hLRtotSoverSplusB -> Draw();
    c5.Print(fname+"all."+extension);
}