Plotting/Modules/ShapeLimit.C

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

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

//#include "TTbar_stuff.C"
using namespace std;

using namespace PlottingSetup;

namespace SUSYScanSpace {
  vector<string> loaded_files;
  int SUSYscantype;
  float SavedMGlu;
  float SavedMLSP;
  string SavedMLSPname;
  string SavedMGluname;
}

const char* strip_flip_away(string flipped_name) {
  int pos = flipped_name.find("flipped_");
  if(pos>=0&&pos<1000) flipped_name=flipped_name.substr(8,1000);
  return flipped_name.c_str();
}

void EliminateNegativeEntries(TH1F *histo) {
  for(int i=1;i<=histo->GetNbinsX();i++) {
    if(histo->GetBinContent(i)<0) histo->SetBinContent(i,0);
  }
}

vector<float> get_expected_points(float observed,int npoints, string RunDirectory, bool doPlot=false) {
  TF1 *gaus = new TF1("gaus","gaus(0)",0,10*observed);
  gaus->SetParameters(1,observed,2*observed);
  gaus->SetParameter(0,1/(gaus->Integral(0,10*observed)));
  float currentpoint=0.01;
  float stepsize=observed/100;
  vector<float> points;
  while(currentpoint<25*observed && points.size()<npoints) {
    
    if(gaus->Integral(0,currentpoint)>((points.size()+1.0)/(npoints+2.0))) {
      if(Verbosity>0) dout << "Added points for calculation at " << currentpoint << " (integral is " << gaus->Integral(0,currentpoint) << ")" << endl;
      points.push_back(currentpoint);
      if(points.size()==npoints) break;
    }
    currentpoint+=stepsize;
  }
  if(doPlot) {
    gaus->SetName("Expected limit computation points");
    gaus->SetTitle("Expected limit computation points");
    TCanvas *can = new TCanvas("can","can");
    gaus->Draw();
    TLine *lines[npoints];
    for(int i=0;i<npoints;i++) {
      lines[i] = new TLine(points[i],0,points[i],gaus->GetMaximum());
      lines[i]->SetLineColor(kBlue);
      lines[i]->Draw();
    }
    can->SaveAs((RunDirectory+"/DistributionOfComputedPoints.png").c_str());
    delete can;
    for(int i=0;i<npoints;i++) delete lines[i];
  }
  delete gaus;
  
  return points;
}

void prepare_MC_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
  TH1F *dataob = (TH1F*)f->Get("data_obs");
  TH1F *signal = (TH1F*)f->Get("signal");
  assert(dataob);
  assert(signal);

  write_warning(__FUNCTION__,"The following two defs for th1's are fake");TH1F *Tbackground; TH1F *Zbackground;
  write_warning(__FUNCTION__,"Not correctly implemented yet for MC");
  ofstream datacard;
  write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
  datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
  dout << "Writing this to a file.\n";
  datacard << "imax 1\n"; // number of channels
  datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
  datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
  datacard << "---------------\n";
  datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
  datacard << "---------------\n";
  datacard << "bin 1\n";
  datacard << "observation "<<dataob->Integral()<<"\n";
  datacard << "------------------------------\n";
  datacard << "bin             1          1          1\n";
  datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
  datacard << "process         0          1          2\n";
  datacard << "rate            "<<signal->Integral()<<"         "<<Tbackground->Integral()<<"         "<<Zbackground->Integral()<<"\n";
  datacard << "--------------------------------\n";
  datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
  datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
  datacard << "Stat   shape    -          -          1    statistical uncertainty (zjets)\n";
  datacard << "Sys    shape    -          1          -    systematic uncertainty on ttbar\n";
  datacard << "Sys    shape    -          -          1    systematic uncertainty on zjets\n";
  datacard << "Stat   shape    1          -          -    statistical uncertainty (signal)\n";
  datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
  datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
  if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
  datacard.close();
}

void prepare_datadriven_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
  TH1F *dataob = (TH1F*)f->Get("data_obs");
  TH1F *signal = (TH1F*)f->Get("signal");
  TH1F *Tbackground = (TH1F*)f->Get("TTbarBackground");
  TH1F *Zbackground  = (TH1F*)f->Get("ZJetsBackground");
  
  assert(dataob);
  assert(signal);
  assert(Tbackground);
  assert(Zbackground);
 

   ofstream datacard;
   write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
   datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
   datacard << "Writing this to a file.\n";
   datacard << "imax " << dataob->GetNbinsX() << "\n"; // number of channels
   datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
   datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
   datacard << "---------------\n";
   datacard << "bin\t";
   for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << " \t";
   datacard << "\n";
   datacard << "observation\t";
   for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinContent(i) << " \t";
   datacard<<"\n";
   datacard << "------------------------------\n";
   datacard << "bin\t";
   for(int i=1;i<=dataob->GetNbinsX();i++) {
     datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" << 
                 " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" << 
                 " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t";
   }
   datacard << "\n";
   datacard << "process\t";
   for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t signal \t TTbarBackground \t ZJetsBackground";
   datacard << "\n";
   datacard << "process\t";
   for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t 0 \t 1 \t 2";
   datacard << "\n";
   
   
   datacard << "rate\t";
   for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t " << signal->GetBinContent(i) << " \t " << Tbackground->GetBinContent(i) << " \t " << Zbackground->GetBinContent(i) << "\t";
   datacard<<"\n";
   
   datacard << "--------------------------------\n";
   
   
   datacard << "lumi     lnN    \t";
   for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+ PlottingSetup::lumiuncert << "\t - \t -";
   datacard << "       luminosity uncertainty\n"; // only affects MC -> signal!
   
   // Statistical uncertainty
   datacard << "Stat     lnN    \t";
   for(int i=1;i<=dataob->GetNbinsX();i++) {
     //Signal
     float central = signal->GetBinContent(i);
     float up      = ((TH1F*)f->Get("signal_StatDown"))->GetBinContent(i);
     float down    = ((TH1F*)f->Get("signal_StatUp"))->GetBinContent(i);
     float suncert=0;
     if(central>0) {
       if(abs(up-central)>abs(down-central)) suncert=abs(up-central)/central;
       else suncert=abs(central-down)/central;
     }

     //TTbar
     central = Tbackground->GetBinContent(i);
     up      = ((TH1F*)f->Get("TTbarBackground_StatUp"))->GetBinContent(i);
     down    = ((TH1F*)f->Get("TTbarBackground_StatDown"))->GetBinContent(i);
     float tuncert=0;
     if(central>0) {
       if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
       else tuncert=abs(central-down)/central;
     }
     //ZJets
     central = Zbackground->GetBinContent(i);
     up      = ((TH1F*)f->Get("ZJetsBackground_StatUp"))->GetBinContent(i);
     down    = ((TH1F*)f->Get("ZJetsBackground_StatDown"))->GetBinContent(i);
     float zuncert=0;
     if(central>0) {
       if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
       else zuncert=abs(central-down)/central;
     }
     datacard << " " << 1+suncert << " \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
   }
   
   datacard << "       statistical uncertainty\n";
   
   // Statistical uncertainty
   datacard << "Sys     lnN    \t";
   for(int i=1;i<=dataob->GetNbinsX();i++) {
     float central = Tbackground->GetBinContent(i);
     float up      = ((TH1F*)f->Get("TTbarBackground_SysUp"))->GetBinContent(i);
     float down    = ((TH1F*)f->Get("TTbarBackground_SysDown"))->GetBinContent(i);
     float tuncert=0;
     if(central>0) {
       if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
       else tuncert=abs(central-down)/central;
     }
     central = Zbackground->GetBinContent(i);
     up      = ((TH1F*)f->Get("ZJetsBackground_SysUp"))->GetBinContent(i);
     down    = ((TH1F*)f->Get("ZJetsBackground_SysDown"))->GetBinContent(i);
     float zuncert=0;
     if(central>0) {
       if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
       else zuncert=abs(central-down)/central;
     }
     datacard << " - \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
   }
   
   datacard << "       systematic uncertainty\n";
   
   
   float JESup = ((TH1F*)f->Get("signal_JESUp"))->Integral();
   float JESdn = ((TH1F*)f->Get("signal_JESDown"))->Integral();
   float central = signal->Integral();
   float uJES=0;
   if(abs(JESup-central)>abs(JESdn-central)) uJES=abs(JESup-central)/central;
   else uJES=abs(central-JESdn)/central;
   datacard << "JES    lnN";
   for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJES << "\t - \t  - \t";
   datacard << "uncertainty on Jet Energy Scale\n";
   
   datacard << "JSU      lnN    ";
   for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJSU << "\t - \t  - \t";
   datacard << "JZB Scale Uncertainty\n";
   
   if(uPDF>0) {
     datacard << "PDF      lnN    ";
     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uPDF << "\t - \t - \t";
     datacard << "uncertainty from PDFs\n";
   }
   
   datacard.close();
}
  
void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
 
 if(FullMCAnalysis) {
   //dealing with MC analysis - need to use shapes.
   prepare_MC_datacard(RunDirectory,f,uJSU,uPDF);
 } else {
   //doing mutibin analysis
   prepare_datadriven_datacard(RunDirectory,f,uJSU,uPDF);
 }
   
}

float QuickDrawNevents=0;

TH1F* QuickDraw(TTree *events, string hname, string variable, vector<float> binning, string xlabel, string ylabel, TCut cut, string addcut, bool dataormc, float luminosity, map <  pair<float, float>, map<string, float>  > &xsec) {
  TH1F *histo = new TH1F(hname.c_str(),hname.c_str(),binning.size()-1,&binning[0]);
  histo->Sumw2();
  stringstream drawcommand;
  drawcommand << variable << ">>" << hname;
  if(SUSYScanSpace::SUSYscantype!=mSUGRA) {
    events->Draw(drawcommand.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight);
  } else {
      float totalxs=0;
    for(int i=0;i<12;i++) {
      stringstream drawcommand2;
      drawcommand2 << variable << ">>+" << hname;
      stringstream mSUGRAweight;
      float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i);
      totalxs+=xs;
      mSUGRAweight << "(0.5*(2*(process==" << i << "))*" << xs << ")";
      events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str());
    }
    histo->Scale(1.0/totalxs); // to get back to 1 /pb normalization
  }
  events->Draw("eventNum",addcut.c_str(),"goff");
  float nevents = events->GetSelectedRows();
  QuickDrawNevents=nevents;
  histo->Scale(luminosity/nevents); // this will give a histogram which is normalized to 1 pb so that any limit will be with respect to 1 pb
  
  return histo;
}


void SQRT(TH1F *h) {
  for (int i=1;i<=h->GetNbinsX();i++) {
    h->SetBinContent(i,TMath::Sqrt(h->GetBinContent(i)));
  }
}

TH1F* Multiply(TH1F *h1, TH1F *h2) {
  TH1F *h = (TH1F*)h1->Clone();
  for(int i=1;i<=h1->GetNbinsX();i++) {
    h->SetBinContent(i,h1->GetBinContent(i) * h2->GetBinContent(i));
  }
  return h;
}


string ReduceNumericHistoName(string origname) {
  int pos = origname.find("__h_");
  if(pos == -1) return origname;
  return origname.substr(0,pos);
}
  

void generate_mc_shapes(bool dotree, TTree *signalevents, string mcjzb,string datajzb,vector<float> binning, TCut weight, TCut JetCut, string addcut, string identifier, map <  pair<float, float>, map<string, float>  > &xsec) {
    //weight can be varied to take care of efficiency errors (weightEffUp, weightEffDown)
  vector<float> mbinning;
  mbinning.push_back(-8000);
  for(int i=binning.size()-1;i>=0;i--) mbinning.push_back(-binning[i]);
  mbinning.push_back(0);
  for(int i=0;i<binning.size();i++) mbinning.push_back(binning[i]);
  mbinning.push_back(8000);

  
  TCanvas *shapecan = new TCanvas("shapecan","shapecan");
  TH1F* Signal;
  
  stringstream sInverseCutWeight;
  sInverseCutWeight << "(1.0/" << (const char*)cutWeight<<")";  
  TCut InverseCutWeight(sInverseCutWeight.str().c_str());
  if(weight==cutWeight) InverseCutWeight = TCut("1.0");  
  if(!dotree) {
      //dealing with ALL MC samples (when we save the standard file)
      THStack McPredicted = allsamples.DrawStack("McPred",mcjzb,mbinning, "JZB [GeV]", "events", ((cutmass&&cutOSSF&&JetCut)*(weight*InverseCutWeight)),mc,luminosity);
      int nhists = McPredicted.GetHists()->GetEntries();
      for (int i = 0; i < nhists; i++) { 
          TH1* hist = static_cast<TH1*>(McPredicted.GetHists()->At(i));
//          cout << hist->GetName() << " : " << hist->Integral() << endl;
//          cout << "     " << ReduceNumericHistoName(hist->GetName()) << endl;
          if(identifier=="StatUp") {
            float appliedweight = hist->Integral() / hist->GetEntries();
            for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)+appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
          }
          if(identifier=="StatDown") {
            float appliedweight = hist->Integral() / hist->GetEntries();
            for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)-appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
          }
          hist->SetName(("mc_"+ReduceNumericHistoName(hist->GetName())+"_"+identifier).c_str());
          hist->Write();
      }
  } else {
      string histoname="mc_signal";
      if(identifier!="") histoname=histoname+"_"+identifier;
      Signal = QuickDraw(signalevents,histoname,mcjzb,binning,"JZB", "events",((cutmass&&cutOSSF&&JetCut&&basiccut))*(weight*InverseCutWeight),addcut,mc,luminosity,xsec);
      if(identifier=="StatUp") {
        float appliedweight = Signal->Integral() / Signal->GetEntries();
        for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)+appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
      }
      if(identifier=="StatDown") {
        float appliedweight = Signal->Integral() / Signal->GetEntries();
        for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)-appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
      }
      Signal->Write();
  }
  
    /*
   * 
   * Jet energy scale (easy, use different JetCut)
   * JZB Scale Uncertainty (will be a number - signal only)
   * Efficiency (will be weightEffUp, weightEffDown)
   * PDFs (signal only) -> Will be a number yet again, computed in the traditional way
   */
  delete shapecan;

}
  
  
void generate_shapes_for_systematic(bool signalonly, bool dataonly,TFile *limfile, TTree *signalevents, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan, string addcut, map <  pair<float, float>, map<string, float>  > &xsec) {
  
  binning.push_back(8000);
  limfile->cd();
  dout << "Creating shape templates ";
  if(identifier!="") dout << "for "<<identifier;
  dout << " ... " << endl;
  
  TCut limitnJetcut;
  
  if(JES==noJES) limitnJetcut=cutnJets;
  else {
    if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
    if(JES==JESup) limitnJetcut=cutnJetsJESup;
  }
  
  float zjetsestimateuncert=zjetsestimateuncertONPEAK;
  float emuncert=emuncertONPEAK;
  float emsidebanduncert=emsidebanduncertONPEAK;
  float eemmsidebanduncert=eemmsidebanduncertONPEAK;
  
  if(!PlottingSetup::RestrictToMassPeak) {
    zjetsestimateuncert=zjetsestimateuncertOFFPEAK;
    emuncert=emuncertOFFPEAK;
    emsidebanduncert=emsidebanduncertOFFPEAK;
    eemmsidebanduncert=eemmsidebanduncertOFFPEAK;
  } 
    
  if(signalonly) {
    dout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: '" << identifier << "')" << endl;
    TH1F *ZOSSFP = QuickDraw(signalevents,"ZOSSFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
    TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
    TH1F *ZOSOFP;
    TH1F *ZOSOFN;

    if(!PlottingSetup::FullMCAnalysis) {
      ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
      ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
    }
    
    TH1F *SBOSSFP;
    TH1F *SBOSOFP;
    TH1F *SBOSSFN;
    TH1F *SBOSOFN;
    
    if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
      SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
      SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
      SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
      SBOSOFN = QuickDraw(signalevents,"SBOSOFN","-"+mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
    }
    
    string signalname="signal";
    if(identifier!="") {
      signalname="signal_"+identifier;
    }
    
    TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
    TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
    
    TH1F *flippedLobs = new TH1F("flippedLobs","flippedLobs",binning.size()-1,&binning[0]);
    TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
    
    Lobs->Add(ZOSSFP);
    if(!PlottingSetup::FullMCAnalysis) Lpred->Add(ZOSSFN);
    
    dout << "SITUATION FOR SIGNAL: " << endl;
    
    
    if(PlottingSetup::FullMCAnalysis) {
      dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << endl;
    } else {
      dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
      dout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
      if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
        dout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
        dout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
      }
    }

    
    flippedLobs->Add(ZOSSFN);
    if(!PlottingSetup::FullMCAnalysis) flippedLpred->Add(ZOSSFP);
    
    if(!PlottingSetup::FullMCAnalysis) {
      if(PlottingSetup::RestrictToMassPeak) {
        Lpred->Add(ZOSOFP,1.0/3);
        Lpred->Add(ZOSOFN,-1.0/3);
        Lpred->Add(SBOSSFP,1.0/3);
        Lpred->Add(SBOSSFN,-1.0/3);
        Lpred->Add(SBOSOFP,1.0/3);
        Lpred->Add(SBOSOFN,-1.0/3);
        
        //flipped prediction
        flippedLpred->Add(ZOSOFP,-1.0/3);
        flippedLpred->Add(ZOSOFN,1.0/3);
        flippedLpred->Add(SBOSSFP,-1.0/3);
        flippedLpred->Add(SBOSSFN,1.0/3);
        flippedLpred->Add(SBOSOFP,-1.0/3);
        flippedLpred->Add(SBOSOFN,1.0/3);
      } else {
        Lpred->Add(ZOSOFP,1.0);
        Lpred->Add(ZOSOFN,-1.0);
        
        //flipped prediction
        flippedLpred->Add(ZOSOFP,-1.0);
        flippedLpred->Add(ZOSOFN,1.0);
      }
    }
    
    TH1F *signal = (TH1F*)Lobs->Clone("signal");
    if(!PlottingSetup::FullMCAnalysis) signal->Add(Lpred,-1);
    signal->SetName(signalname.c_str());
    signal->SetTitle(signalname.c_str());
    signal->Write();
    
    TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
    if(!PlottingSetup::FullMCAnalysis) flippedsignal->Add(flippedLpred,-1);
    flippedsignal->SetName(("flipped_"+signalname).c_str());
    flippedsignal->Write();
    
    if(identifier=="") {
      TH1F *signalStatUp = (TH1F*)signal->Clone();
      signalStatUp->SetName(((string)signal->GetName()+"_StatUp").c_str());
      signalStatUp->SetTitle(((string)signal->GetTitle()+"_StatUp").c_str());
      
      TH1F *signalStatDn = (TH1F*)signal->Clone();
      signalStatDn->SetName(((string)signal->GetName()+"_StatDown").c_str());
      signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
      
      for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
        float staterr;
        if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
        else staterr = TMath::Sqrt(Lobs->GetBinContent(i));
        
        if(!PlottingSetup::FullMCAnalysis) {
          dout << "Stat err in bin " << i << " : " << staterr << endl;
          dout << "    prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
          dout << "    we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
        }
        if(signal->GetBinContent(i)-staterr>0) signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
        else signalStatDn->SetBinContent(i,0);
        signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr);
        signal->SetBinError(i,staterr);
      }
      
      signalStatDn->Write();
      signalStatUp->Write();
      
      delete signalStatDn;
      delete signalStatUp;

      TH1F *flippedsignalStatUp = (TH1F*)flippedsignal->Clone();
      flippedsignalStatUp->SetName(((string)flippedsignal->GetName()+"_StatUp").c_str());
      flippedsignalStatUp->SetTitle(((string)flippedsignal->GetTitle()+"_StatUp").c_str());
      
      TH1F *flippedsignalStatDn = (TH1F*)flippedsignal->Clone();
      flippedsignalStatDn->SetName(((string)flippedsignal->GetName()+"_StatDown").c_str());
      flippedsignalStatDn->SetTitle(((string)flippedsignal->GetTitle()+"_StatDown").c_str());
      
      for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
        float staterr;
        if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
        else staterr = TMath::Sqrt(flippedLobs->GetBinContent(i));
        if(flippedsignal->GetBinContent(i)-staterr>0) flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
        else flippedsignalStatDn->SetBinContent(i,0);
        flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+staterr);
        flippedsignal->SetBinError(i,staterr);
      }
      
      flippedsignalStatDn->Write();
      flippedsignalStatUp->Write();
      
      delete flippedsignalStatDn;
      delete flippedsignalStatUp;
    }
    
    delete ZOSSFP;
    delete ZOSOFP;
    delete ZOSSFN;
    delete ZOSOFN;
    
    if(PlottingSetup::RestrictToMassPeak) {
      delete SBOSSFP;
      delete SBOSOFP;
      delete SBOSSFN;
      delete SBOSOFN;
    }
    
    delete Lobs;
    delete Lpred;
    delete flippedLobs;
    delete flippedLpred;
  }
  
    
  if(dataonly) {
    dout << "Processing data with datajzb: " << datajzb << endl;
    TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
    TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
    TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
    TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
    
    TH1F *SBOSSFP;
    TH1F *SBOSOFP;
    TH1F *SBOSSFN;
    TH1F *SBOSOFN;
    
    if(PlottingSetup::RestrictToMassPeak) {
      SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
      SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
      SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
      SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
    }
    
    string obsname="data_obs";
    string predname="background";
    string Zpredname="ZJetsBackground";
    string Tpredname="TTbarBackground";
    if(identifier!="") {
      obsname=("data_"+identifier);
      predname=("background_"+identifier);
      Zpredname=("ZJetsBackground_"+identifier);
      Tpredname=("TTbarBackground_"+identifier);
    }
    
    TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
    obs->SetName(obsname.c_str());
    obs->Write();
    
    TH1F *flippedobs = (TH1F*)ZOSSFN->Clone("flipped_observation");
    flippedobs->SetName(("flipped_"+obsname).c_str());
    flippedobs->Write();
    
    TH1F *Tpred = (TH1F*)ZOSSFN->Clone("TTbarprediction");
    TH1F *Zpred = (TH1F*)ZOSSFN->Clone("ZJetsprediction");
    
    TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
    
    TH1F *flippedTpred = (TH1F*)ZOSSFP->Clone("flipped_TTbarprediction");
    TH1F *flippedZpred = (TH1F*)ZOSSFP->Clone("flipped_ZJetsprediction");
    
    TH1F *flippedpred = (TH1F*)ZOSSFP->Clone("flipped_prediction");
    if(PlottingSetup::RestrictToMassPeak) {
      pred->Add(ZOSOFP,1.0/3);
      pred->Add(ZOSOFN,-1.0/3);
      pred->Add(SBOSSFP,1.0/3);
      pred->Add(SBOSSFN,-1.0/3);
      pred->Add(SBOSOFP,1.0/3);
      pred->Add(SBOSOFN,-1.0/3);
      
      Tpred->Add(ZOSSFN,-1.0);
      Tpred->Add(ZOSOFP,1.0/3);
      Tpred->Add(SBOSSFP,1.0/3);
      Tpred->Add(SBOSOFP,1.0/3);
      
      Zpred->Add(ZOSOFN,-1.0/3);
      Zpred->Add(SBOSSFN,-1.0/3);
      Zpred->Add(SBOSOFN,-1.0/3);
      
      //flipped
      flippedpred->Add(ZOSOFP,-1.0/3);
      flippedpred->Add(ZOSOFN,1.0/3);
      flippedpred->Add(SBOSSFP,-1.0/3);
      flippedpred->Add(SBOSSFN,1.0/3);
      flippedpred->Add(SBOSOFP,-1.0/3);
      flippedpred->Add(SBOSOFN,1.0/3);
      
      flippedTpred->Add(ZOSSFP,-1.0);
      flippedTpred->Add(ZOSOFN,1.0/3);
      flippedTpred->Add(SBOSSFN,1.0/3);
      flippedTpred->Add(SBOSOFN,1.0/3);
      
      flippedZpred->Add(ZOSOFP,-1.0/3);
      flippedZpred->Add(SBOSSFP,-1.0/3);
      flippedZpred->Add(SBOSOFP,-1.0/3);
    } else {
      pred->Add(ZOSOFP,1.0);
      pred->Add(ZOSOFN,-1.0);
      
      Tpred->Add(ZOSSFN,-1.0);
      Tpred->Add(ZOSOFP,1.0);
      
      Zpred->Add(ZOSOFN,-1.0);
      
      //flipped
      flippedpred->Add(ZOSOFN,1.0);
      flippedpred->Add(ZOSOFP,-1.0);
      
      flippedTpred->Add(ZOSSFP,-1.0);
      flippedTpred->Add(ZOSOFN,1.0);
      
      flippedZpred->Add(ZOSOFP,-1.0);
    }
    
    pred->SetName(predname.c_str());
    Tpred->SetName(Tpredname.c_str());
    Zpred->SetName(Zpredname.c_str());
    
    flippedpred->SetName(("flipped_"+predname).c_str());
    flippedTpred->SetName(("flipped_"+Tpredname).c_str());
    flippedZpred->SetName(("flipped_"+Zpredname).c_str());

    pred->Write();
    Tpred->Write();
    Zpred->Write();
    
    flippedpred->Write();
    flippedTpred->Write();
    flippedZpred->Write();
    
    if(identifier=="") {
      float stretchfactor=1.0;
      bool isonpeak=false;
      if(PlottingSetup::RestrictToMassPeak) {
        stretchfactor=1.0/3.0;
        isonpeak=true;
      }
      
      //--------------------------------------------------------statistical uncertainty
      
      TH1F *predstaterr = (TH1F*)ZOSSFN->Clone("predstaterr");
      predstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
      predstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
      if(isonpeak) predstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
      if(isonpeak) predstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
      if(isonpeak) predstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
      if(isonpeak) predstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
      SQRT(predstaterr);
      TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
      bgStatUp->Add(predstaterr);
      EliminateNegativeEntries(bgStatUp);
      bgStatUp->Write();
      TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDown");
      bgStatDn->Add(predstaterr,-1);
      EliminateNegativeEntries(bgStatDn);
      bgStatDn->Write();
//      delete bgStatDn;
//      delete bgStatUp;
      delete predstaterr;
      
      
      TH1F *flippedpredstaterr = (TH1F*)ZOSSFP->Clone("predstaterr");
      flippedpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
      flippedpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
      if(isonpeak) flippedpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
      if(isonpeak) flippedpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
      if(isonpeak) flippedpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
      if(isonpeak) flippedpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
      SQRT(flippedpredstaterr);
      TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
      fbgStatUp->Add(predstaterr);
      EliminateNegativeEntries(fbgStatUp);
      fbgStatUp->Write();
      TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDown");
      fbgStatDn->Add(predstaterr,-1);
      EliminateNegativeEntries(fbgStatDn);
      fbgStatDn->Write();
//      delete fbgStatDn;
//      delete fbgStatUp;
      delete predstaterr;
      
      TH1F *Tpredstaterr = (TH1F*)ZOSOFP->Clone("Tpredstaterr");
      Tpredstaterr->Scale(stretchfactor*stretchfactor);
      if(isonpeak) Tpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
      if(isonpeak) Tpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
      SQRT(Tpredstaterr);
      TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
      TpredStatUp->Add(Tpredstaterr);
      EliminateNegativeEntries(TpredStatUp);
      TpredStatUp->Write();
      TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDown");
      TpredStatDn->Add(Tpredstaterr,-1);
      EliminateNegativeEntries(TpredStatDn);
      TpredStatDn->Write();
//      delete TpredStatDn;
//      delete TpredStatUp;
      delete Tpredstaterr;
      
      TH1F *fTpredstaterr = (TH1F*)ZOSOFN->Clone("flipped_Tpredstaterr");
      fTpredstaterr->Scale(stretchfactor*stretchfactor);
      if(isonpeak) fTpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
      if(isonpeak) fTpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
      SQRT(fTpredstaterr);
      TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
      fTpredStatUp->Add(fTpredstaterr);
      EliminateNegativeEntries(fTpredStatUp);
      fTpredStatUp->Write();
      TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDown");
      fTpredStatDn->Add(fTpredstaterr,-1);
      EliminateNegativeEntries(fTpredStatDn);
      fTpredStatDn->Write();
//      delete fTpredStatDn;
//      delete fTpredStatUp;
      delete fTpredstaterr;
      
        
      TH1F *Zpredstaterr = (TH1F*)ZOSSFN->Clone("ZJetsstaterr");
      Zpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
      if(isonpeak) Zpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
      if(isonpeak) Zpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
      SQRT(Zpredstaterr);
      TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
      ZpredStatUp->Add(Zpredstaterr);
      EliminateNegativeEntries(ZpredStatUp);
      ZpredStatUp->Write();
      TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDown");
      ZpredStatDn->Add(Zpredstaterr,-1);
      EliminateNegativeEntries(ZpredStatDn);
      ZpredStatDn->Write();
//      delete ZpredStatDn;
//      delete ZpredStatUp;
      delete Zpredstaterr;
      
      TH1F *fZpredstaterr = (TH1F*)ZOSSFP->Clone("flipped_ZJetsstaterr");
      Zpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
      if(isonpeak) fZpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
      if(isonpeak) fZpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
      SQRT(fTpredstaterr);
      TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
      fZpredStatUp->Add(fZpredstaterr);
      EliminateNegativeEntries(fZpredStatUp);
      fZpredStatUp->Write();
      TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDown");
      fZpredStatDn->Add(fZpredstaterr,-1);
      EliminateNegativeEntries(fZpredStatDn);
      fZpredStatDn->Write();
//      delete fZpredStatDn;
//      delete fZpredStatUp;
      delete fZpredstaterr;
      
      //--------------------------------------------------------systematic uncertainty
      
      TH1F *predsyserr = (TH1F*)pred->Clone("SysError");
      predsyserr->Add(pred,-1);//getting everything to zero.
      predsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
      predsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
      predsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
      if(isonpeak) predsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
      if(isonpeak) predsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
      if(isonpeak) predsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
      if(isonpeak) predsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
      SQRT(predsyserr);
      TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
      bgSysUp->Add(predsyserr);
      EliminateNegativeEntries(bgSysUp);
      bgSysUp->Write();
      TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
      bgSysDn->Add(predsyserr,-1);
      EliminateNegativeEntries(bgSysDn);
      bgSysDn->Write();
      delete predsyserr;
      
      TH1F *fpredsyserr = (TH1F*)pred->Clone("SysError");
      fpredsyserr->Add(pred,-1);//getting everything to zero.
      fpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
      fpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
      fpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
      if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
      if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
      if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
      if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
      SQRT(fpredsyserr);
      TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
      fbgSysUp->Add(fpredsyserr);
      EliminateNegativeEntries(fbgSysUp);
      fbgSysUp->Write();
      TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDown");
      fbgSysDn->Add(fpredsyserr,-1);
      EliminateNegativeEntries(fbgSysDn);
      fbgSysDn->Write();
      delete fpredsyserr;
      
      
      TH1F *Tpredsyserr = (TH1F*)Tpred->Clone("SysError");
      Tpredsyserr->Add(Tpred,-1);//getting everything to zero.
      Tpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
      if(isonpeak) Tpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
      if(isonpeak) Tpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
      SQRT(Tpredsyserr);
      TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
      TpredSysUp->Add(Tpredsyserr);
      EliminateNegativeEntries(TpredSysUp);
      TpredSysUp->Write();
      TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
      TpredSysDn->Add(Tpredsyserr,-1);
      EliminateNegativeEntries(TpredSysDn);
      TpredSysDn->Write();
      delete Tpredsyserr;
      
      TH1F *fTpredsyserr = (TH1F*)Tpred->Clone("SysError");
      fTpredsyserr->Add(Tpred,-1);//getting everything to zero.
      fTpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
      if(isonpeak) fTpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
      if(isonpeak) fTpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
      SQRT(fTpredsyserr);
      TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
      fTpredSysUp->Add(fTpredsyserr);
      EliminateNegativeEntries(fTpredSysUp);
      fTpredSysUp->Write();
      TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDown");
      fTpredSysDn->Add(fTpredsyserr,-1);
      EliminateNegativeEntries(fTpredSysDn);
      fTpredSysDn->Write();
      delete fTpredsyserr;
      
      
      //------------------------------------------------------------------------------------------------------------------------
      TH1F *Zpredsyserr = (TH1F*)Zpred->Clone("SysError");
      Zpredsyserr->Add(Zpred,-1);//getting everything to zero.
      Zpredsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
      Zpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
      if(isonpeak) Zpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
      if(isonpeak) Zpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
      SQRT(Zpredsyserr);
      TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
      ZpredSysUp->Add(Zpredsyserr);
      EliminateNegativeEntries(ZpredSysUp);
      ZpredSysUp->Write();
      TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
      ZpredSysDn->Add(Zpredsyserr,-1);
      EliminateNegativeEntries(ZpredSysDn);
      ZpredSysDn->Write();
      delete Zpredsyserr;
      
      TH1F *fZpredsyserr = (TH1F*)Zpred->Clone("SysError");
      fZpredsyserr->Add(Zpred,-1);//getting everything to zero.
      fZpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
      fZpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
      if(isonpeak) fZpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
      if(isonpeak) fZpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
      SQRT(fZpredsyserr);
      TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
      fZpredSysUp->Add(fZpredsyserr);
      EliminateNegativeEntries(fZpredSysUp);
      fZpredSysUp->Write();
      TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDown");
      fZpredSysDn->Add(fZpredsyserr,-1);
      EliminateNegativeEntries(fZpredSysDn);
      fZpredSysDn->Write();
      delete fZpredsyserr;
    }
    
/*if(identifier=="") {
  for(int i=0;i<binning.size()-1;i++) {
    dout << "[ " << binning[i] << " , " << binning[i+1] << "] : O " << obs->GetBinContent(i+1) << " P " << pred->GetBinContent(i+1) << " (Z: " << Zpred->GetBinContent(i+1) << " , T: " << Tpred->GetBinContent(i+1) << ")" << endl;
  }
}*/
    delete ZOSSFP;
    delete ZOSOFP;
    delete ZOSSFN;
    delete ZOSOFN;
    
    if(PlottingSetup::RestrictToMassPeak) {
      delete SBOSSFP;
      delete SBOSOFP;
      delete SBOSSFN;
      delete SBOSOFN;
    }
  }
  
}

void DoFullCls(string RunDirectory,float firstGuess) {
  vector<float> epoints = get_expected_points(firstGuess,25,RunDirectory,true);//get 20 points;
  ofstream RealScript;
  dout << "Going to write the full CLs script to " << RunDirectory << "/LimitScript.sh" << endl;
  RealScript.open ((RunDirectory+"/LimitScript.sh").c_str());
  RealScript << "#!/bin/bash \n";
  RealScript << "\n";
  RealScript << "RunDirectory=" << RunDirectory << "\n";
  RealScript << "CMSSWDirectory=" << PlottingSetup::CMSSW_dir << "\n";
  RealScript << "\n";
  RealScript << "echo \"Going to compute limit in $RunDirectory\"\n";
  RealScript << "ORIGDIR=`pwd`\n";
  RealScript << "\n";
  RealScript << "pushd $CMSSWDirectory\n";
  RealScript << "cd ~/final_production_2011/CMSSW_4_2_8/src/HiggsAnalysis/\n";
  RealScript << "origscramarch=$SCRAM_ARCH\n";
  RealScript << "origbase=$CMSSW_BASE\n";
  RealScript << "export SCRAM_ARCH=slc5_amd64_gcc434\n";
  RealScript << "eval `scram ru -sh`\n";
  RealScript << "\n";
  RealScript << "mkdir -p $RunDirectory\n";
  RealScript << "cd $RunDirectory\n";
  RealScript << "echo `pwd`\n";
  RealScript << "mkdir -p FullCLs\n";
  RealScript << "cd FullCLs\n";
  RealScript << "\n";
  RealScript << "combineEXE=\"${CMSSW_BASE}/bin/${SCRAM_ARCH}/combine\"\n";
  RealScript << "\n";
  RealScript << "time for i in {";
  RealScript << epoints[0]/4 << ","; // adding a VERY VERY low point just to be totally safe
  RealScript << epoints[0]/2 << ","; // adding a VERY low point just to be safe
  for(int i=0;i<epoints.size();i++) RealScript << epoints[i] << ",";
  RealScript << 2*epoints[epoints.size()-1] << ","; // adding a VERY high point just to be safe
  RealScript << 4*epoints[epoints.size()-1] << ","; // adding a VERY VERY high point just to be totally safe
  RealScript << "}; do time $combineEXE ../susydatacard.txt -M HybridNew --freq -T 500 --clsAcc 0 -v 0 -n TestPoint --saveHybridResult --saveToys -s -1 -i 8 --singlePoint $i; done\n"; // 500 toys, epoints.size+2 testpoints, 8 iterations
  RealScript << "\n";
  RealScript << "hadd FullCLs.root higgsCombine*.root\n";
  RealScript << "mkdir originalFiles\n";
  RealScript << "mv higgsCombine* originalFiles\n";
  RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root\n";
  RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.5\n";
  RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.84\n";
  RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.16\n";
  RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.975\n";
  RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.025\n";
  RealScript << "\n";
  RealScript << "hadd FinalResult.root higgsCombineTest.HybridNew*\n";
  RealScript << "\n";
  RealScript << "g++ $ORIGDIR/ReadAndSave.C -o ReadAndSave.exec `root-config --glibs --cflags` \n";
  RealScript << "./ReadAndSave.exec FinalResult.root " << RunDirectory << "/FullCLsShapeDropletResult.txt \n";
  RealScript << "\n";
  RealScript << "#####\n";
  RealScript << "echo \"Finalizing ...\"\n";
  RealScript << "cd $ORIGDIR\n";
  RealScript << "echo $ORIGDIR\n";
  RealScript << "ls -ltrh ../SandBox/ | grep combine\n";
  RealScript << "export SCRAM_ARCH=$origscramarch\n";
  RealScript << "exit 0\n";
  RealScript << "\n";
  RealScript.close();
  gSystem->Exec(("bash "+RunDirectory+"/LimitScript.sh").c_str());
}
    
    
ShapeDroplet LimitsFromShapes(bool asymptotic, float firstGuess, TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
  map <  pair<float, float>, map<string, float>  >  xsec;
  if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
  
  time_t timestamp;
  tm *now;
  timestamp = time(0);
  now = localtime(&timestamp);
  stringstream RunDirectoryS;
  RunDirectoryS << PlottingSetup::cbafbasedir << "/exchange/ShapeComputation___" << name << "__" << now->tm_hour << "h_" << now->tm_min << "m_" << now->tm_sec << "s___" << rand();
  string RunDirectory = RunDirectoryS.str();
  
  ensure_directory_exists(RunDirectory);
  
  TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str(),"READ");
  if(datafile->IsZombie()) {
    write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
    assert(!datafile->IsZombie());
  }
  dout << "Run Directory: " << RunDirectory << endl;
  TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
  
  TIter nextkey(datafile->GetListOfKeys());
  TKey *key;
  while ((key = (TKey*)nextkey()))
  {
      TH1F *obj =(TH1F*) key->ReadObj();
      limfile->cd();
      obj->Write();
  }
  
  TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
  
  bool signalonly=true;
  bool dataonly=false;
  
  generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"",mcjzb,datajzb,noJES,jzbbins,limcan,addcut,xsec);
  generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,addcut,xsec);
  generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,addcut,xsec);
  
  float JZBScaleUncert=0.05; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
  float scaledown=0,scaleup=0,scalesyst=0;
  float mceff=0,mcefferr=0;
  int flipped=0;
  int Neventsinfile=10000;//doesn't matter.
  string informalname="ShapeAnalysis";
  
  bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
  
  MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
  if(mceff<0) {
    flipped=1;
    write_info(__FUNCTION__,"Doing flipping!");
  }
  doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
  float PDFuncert=0;
  int NPdfs = get_npdfs(events);
  write_warning(__FUNCTION__,"Deactivated PDFs temporarily for exp limits!!!");
//  if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);

  float JZBscale=scaledown;
  if(scaleup>scaledown) JZBscale=scaleup;
  dout << endl;
  dout << endl;
  dout << "Will use the following scalar (non-shape) systematics : " << endl;
  dout << "       JZB Scale (JSU) : " << JZBscale << endl;
  dout << "       PDF             : " << PDFuncert << endl;
  
  float SignalIntegral; 
  if(flipped) {
    TH1F *flisi = (TH1F*)(limfile->Get("flipped_signal"));
    SignalIntegral= flisi ->Integral();
    delete flisi;
  }
  else {
    TH1F *flisi = (TH1F*)(limfile->Get("signal"));
    SignalIntegral= flisi ->Integral();
    delete flisi;
  }
    
  TFile *final_limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
  
  TIter nnextkey(limfile->GetListOfKeys());
  TKey *nkey;
  
  if(SignalIntegral==0) SignalIntegral=1; //case when the integral is zero - don't try to normalize, just leave it as it is
  
  while ((nkey = (TKey*)nnextkey()))
  {
      TH1F *obj =(TH1F*) nkey->ReadObj();
      if(flipped&&!Contains(obj->GetName(),"flipped")) continue;
      if(!flipped&&Contains(obj->GetName(),"flipped")) continue;
      if(flipped) obj->SetName(strip_flip_away(obj->GetName()));
      final_limfile->cd();
      if(Contains(obj->GetName(),"signal")) { 
        obj->Scale(1.0/SignalIntegral);
      }
      obj->Write();
  }

  prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
  
  final_limfile->Close();
  limfile->Close();
  dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl;
  stringstream command;
  string DropletLocation;
  int CreatedModelFileExitCode;
  if(asymptotic) {
    if(firstGuess>0) {
      command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
    } else {
      command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
    }
    dout <<"Going to run : " << command.str() << endl;
    CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
    dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
    if(!(CreatedModelFileExitCode==0)) {
      write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
      ShapeDroplet alpha;
      alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ... 
      alpha.SignalIntegral=1;
      return alpha;
    }
    DropletLocation=RunDirectory+"/ShapeDropletResult.txt";
  }
  else {
    DoFullCls(RunDirectory,firstGuess);
    DropletLocation=RunDirectory+"/FullCLsShapeDropletResult.txt";
  }
  ShapeDroplet alpha;
  alpha.readDroplet(DropletLocation);
  alpha.PDF=PDFuncert;
  alpha.JSU=JZBscale;
  alpha.SignalIntegral=SignalIntegral;
  dout << "Done reading limit information from droplet" << endl;
  dout << alpha << endl;
  if(asymptotic) {
    dout << "Doing asymptotic limits, therefore adding an additional round!" << endl;
    command.str("");
    command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << alpha.expected; // ASYMPTOTIC LIMITS WITH FIRST GUESS
    CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
    dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
    if(!(CreatedModelFileExitCode==0)) {
      write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
      ShapeDroplet alpha;
      alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ... 
      alpha.SignalIntegral=1;
      return alpha;
    }
    alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
    alpha.PDF=PDFuncert;
    alpha.JSU=JZBscale;
    alpha.SignalIntegral=SignalIntegral;
    dout << "Done reading updated limit information from droplet" << endl;
    dout << alpha << endl;
  }
    
  dout << "Everything is saved in " << RunDirectory << endl;
  dout << "Cleaning up ... " << std::flush;
/*  dout << "   1) Make sure models directory exists ... " << std::flush;
  gSystem->Exec("mkdir -p models/");
  dout << " ok!" << endl;
  dout << "   2) Deleting any previous model files with the same name ... " << std::flush;
  stringstream delcommand;
  delcommand << "rm models/model_" << name << ".root";
  gSystem->Exec(delcommand.str().c_str());
  stringstream delcommand2;
  delcommand2 << "rm models/model_" << name << ".txt";
  gSystem->Exec(delcommand2.str().c_str());
  dout << " ok!" << endl;
  dout << "   3) Transfering the new model files ... " << std::flush;
  stringstream copycommand;
  copycommand << "cp " << RunDirectory << "/model.root models/model_" << name << ".root";
  gSystem->Exec(copycommand.str().c_str());
  stringstream copycommand2;
  copycommand2 << "cp " << RunDirectory << "/susydatacard.txt models/model_" << name << ".txt";
  gSystem->Exec(copycommand2.str().c_str());
  stringstream copycommand3;
  copycommand3 << "cp " << RunDirectory << "/limitfile.root models/model_" << name << "_histo.root";
  gSystem->Exec(copycommand3.str().c_str());
  dout << " ok!" << endl;
  dout << "   4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;*/
  gSystem->Exec(("rm -r "+RunDirectory).c_str());
  dout << " ok!" << endl;
  delete limcan;
  return alpha;
}

void PrepareDataShapes(string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
  
  map <  pair<float, float>, map<string, float>  >  xsec;
  TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
  TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
  TTree *faketree;
  bool dataonly=true;
  bool signalonly=false;
  
  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
  
  generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","",xsec);
  generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESup,"","JESUp",xsec);
  generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESdown,"","JESDown",xsec);
  generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffUp"),cutnJets,"","EffUp",xsec);
  generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffDown"),cutnJets,"","EffDown",xsec);
  generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatUp",xsec);
  generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatDown",xsec);
  /*
  datafile->Close();
  */
}
  
Revision:	1.12
Committed:	Tue May 15 19:37:47 2012 UTC (12 years, 11 months ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.11:	+125 -17 lines
Log Message:	Updated limit implementation (shape limits via grid of points)