Plotting/Modules/ExclusionPlot.C

#include <iostream>
#include <vector>
#include <sys/stat.h>
#include <sstream>
#include <assert.h>

#include <TCut.h>
#include <TLatex.h>
#include <TROOT.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include <TPRegexp.h>
#include <TMath.h>
#include <stdio.h>
#include <stdlib.h>
#include <TH2.h>
#include <TColor.h>
#include <TStyle.h>
#include <TText.h>
#include <TKey.h>
#include <TPaletteAxis.h>
#include <TGraph.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include "external/ExclusionPlot/Functions.C"

using namespace std;

bool draweachone=false;
bool draw2sigma=true;

float limits_lower_bound=0.05;
float limits_upper_bound=10;
float limits_ratio_lower_bound=0.005;
float limits_ratio_upper_bound=10;

float standardmargin=0.1;
float indentedmargin=0.16;

bool drawefficiencydesertline=false;

bool wrongwaytodothis=true;

string xsecfilename;


void set_range(TH2F *histo, int scantype, bool pushoutlabels);
void smooth_line(TGraph *gr);
void draw_diagonal_xchange(int scantype, std::string scanx);
TH2F* adjust_histo(TH2F *oldhist, TH2F *refhisto);
TGraph* get_mSUGRA_exclusion_line(TH2F *exclusionhisto, int scantype);
TGraph* thin_line(TGraph *gr);
TGraph *MarcosExclusionLine(TH2F *exclusionshape, int scantype);
TH2F* get_XS(string filename, string mass, TH2F *histo);

float getSMSxs(float mlsp,float mglu) {
  TH2F *refh  =  new TH2F("ReferenceHisto","ReferenceHisto",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
  refh->SetBinContent(refh->FindBin(mglu,mlsp),1);//only compute the cross section for our point
  TFile *xsecfile = new TFile((PlottingSetup::cbafbasedir+"/"+PlottingSetup::SMSReferenceXSFile).c_str());
  if(xsecfile->IsZombie()) {
      write_error(__FUNCTION__,"Cross section file is invalid!!!!");
      return -1;
  }
  xsecfile->Close();
  delete xsecfile;
  TH2F *xsec = adjust_histo(get_XS(PlottingSetup::cbafbasedir+"/"+PlottingSetup::SMSReferenceXSFile,"gluino",refh),refh);
  int GlobalBin = xsec->FindBin(mglu,mlsp);
  float refxs=xsec->GetBinContent(GlobalBin);
  delete refh;
  delete xsec;
  return refxs;
}

void write_SMS_text(int scantype, std::string& scanx, float xpos = 0.22 ) {
  string legT5zz="pp #rightarrow  #tilde{g} #tilde{g}, #tilde{g} #rightarrow 2j + #chi^{0}_{2}, #chi^{0}_{2} #rightarrow Z #chi^{0}_{1}";
  if ( scantype==PlottingSetup::GMSB) legT5zz="pp #rightarrow  #tilde{g} #tilde{g}, #tilde{g} #rightarrow 2j + #chi^{0}_{1}, #chi^{0}_{1} #rightarrow Z #tilde{G}";

  string legT5zzl2="m(#tilde{q}) >> m(#tilde{g})";
  if (scantype==PlottingSetup::SMS) legT5zzl2 += ", x = "+scanx;
  string legT5zzl3=" GMSB";

  TText *title = write_text(xpos,0.85,legT5zz);
  title->SetTextAlign(11);
  title->SetTextSize(0.035);
  if(scantype!=PlottingSetup::mSUGRA) title->Draw("same");
  TText *title2 = write_text(xpos,0.79,legT5zzl2);
  title2->SetTextAlign(11);
  title2->SetTextSize(0.035);
  if(scantype!=PlottingSetup::mSUGRA) title2->Draw("same");
  TText *title3 = write_text(0.40,0.79,legT5zzl3);
  title3->SetTextAlign(11);
  title3->SetTextSize(0.035);
  title3->SetTextColor(kRed);
  //  if(scantype==PlottingSetup::GMSB) title3->Draw("same");
}

void draw_diagonal_xchange(int scantype, std::string scanx = "" ) {
  // Line marking the diagonal
  TLine *line;
  float verticaloffset=0.0;
  if(scantype==PlottingSetup::GMSB) verticaloffset=75.0;
  line = new TLine(50.+75.0, 50.0+verticaloffset, 1200., 1200.0-75.0+verticaloffset);
  line->SetLineStyle(7);
  line->SetLineWidth(4);
  //line->Draw("same"); // Do not draw: draw the other one instead
  
  // Add a dashed line to indicate where x changes
  float offset = 0.;
  if ( 0 == scanx.compare("0.5") ) { offset = 91.2/0.5; }
  else if ( 0 == scanx.compare("0.25") ) { offset = 91.2/0.25; }
  else if ( 0 == scanx.compare("0.75") ) { offset = 91.2/0.75; }
  else if ( scantype==PlottingSetup::GMSB) { offset = 0; };
  
  if ( offset>0. ) {
    line->DrawLine(50+offset, 50.0, 1175., 1175.0-offset);
  }else if ( scantype==PlottingSetup::GMSB ) {
    line->DrawLine(100, 100, 1175., 1175.0-offset);
 }
}

void fill_with_text(TGraph *real, TGraph *down, TGraph *up, TVirtualPad *can, int scantype, std::string scanx = "") {
  can->cd();
  float xpos_of_text = 0.22;
  TLegend* this_leg = new TLegend(xpos_of_text,0.6,0.38,0.75);
  //TLegend* this_leg = new TLegend(xpos_of_text,0.55,0.45,0.75,"n_{jets} #geq 3"); // this was the style of the paper. 
  this_leg->SetFillColor(0);
  this_leg->SetBorderSize(0);
  this_leg->SetTextSize(0.035);
  //this_leg->SetTextSize(0.04); // paper style. 
  if(scantype==PlottingSetup::SMS||scantype==PlottingSetup::GMSB) {
    //this_leg->AddEntry(real,"#sigma^{prod} = #sigma^{NLO-QCD}" , "l");
    //this_leg->AddEntry(up,"#sigma^{prod} = 3 #times #sigma^{NLO-QCD}" , "l");
    //this_leg->AddEntry(down,"#sigma^{prod} = 1/3 #times #sigma^{NLO-QCD}" , "l");
    this_leg->AddEntry(real,"#sigma^{NLO-QCD}" , "l");
    this_leg->AddEntry(up,"3 #times #sigma^{NLO-QCD}" , "l");
    this_leg->AddEntry(down,"1/3 #times #sigma^{NLO-QCD}" , "l");
  } else {
    write_warning(__FUNCTION__,"Not implemented yet for mSUGRA");
  }
    
  this_leg->Draw();
  TText *title = write_text(xpos_of_text+0.005,0.52,"JZB");
  title->SetTextSize(0.04);
  title->SetTextAlign(13);
  title->SetTextFont(62);
  title->Draw();
  
  write_SMS_text( scantype, scanx, xpos_of_text );
  
//  //string legT5zz="pp #rightarrow  #tilde{g} #tilde{g}, #tilde{g} #rightarrow 2j + #chi^{0}_{1}, #chi^{0}_{1} #rightarrow Z + #tilde{G}";
//  string legT5zz="pp #rightarrow  #tilde{g} #tilde{g}, #tilde{g} #rightarrow 2j + #chi^{0}_{2}, #chi^{0}_{2} #rightarrow Z #chi^{0}_{1}";
//  string legT5zzl2="m(#tilde{q}) >> m(#tilde{g}), x = "+scanx;
//  string legT5zzl3=" GMSB";
//  TText *title = write_text(xpos_of_text,0.85,legT5zz);
//  title->SetTextAlign(11);
//  title->SetTextSize(0.035);
//  if(scantype!=PlottingSetup::mSUGRA) title->Draw("same");
//  TText *title2 = write_text(xpos_of_text,0.79,legT5zzl2);
//  title2->SetTextAlign(11);
//  title2->SetTextSize(0.035);
//  if(scantype!=PlottingSetup::mSUGRA) title2->Draw("same");
//  TText *title3 = write_text(0.40,0.79,legT5zzl3);
//  title3->SetTextAlign(11);
//  title3->SetTextSize(0.035);
//  title3->SetTextColor(kRed);
////  if(scantype==PlottingSetup::GMSB) title3->Draw("same");
  DrawPrelim();
  draw_diagonal_xchange( scantype, scanx );
}

TH2F* prep_histo(TH2F *oldhist, int scantype) {///DONE
  string name=oldhist->GetName();
  oldhist->SetName(((string)"unprepped"+oldhist->GetName()).c_str());
  TH2F *histo = new TH2F(name.c_str(),name.c_str(),
                         oldhist->GetNbinsX(),oldhist->GetXaxis()->GetBinLowEdge(1),oldhist->GetXaxis()->GetBinLowEdge(oldhist->GetNbinsX())+oldhist->GetXaxis()->GetBinWidth(oldhist->GetNbinsX()),
                         oldhist->GetNbinsY(),oldhist->GetYaxis()->GetBinLowEdge(1),oldhist->GetYaxis()->GetBinLowEdge(oldhist->GetNbinsY())+oldhist->GetYaxis()->GetBinWidth(oldhist->GetNbinsY()));
  
  for(int ix=1;ix<=oldhist->GetNbinsX();ix++) {
    for(int iy=1;iy<=oldhist->GetNbinsX();iy++) {
      histo->SetBinContent(ix,iy,oldhist->GetBinContent(ix,iy));
    }
  }
  return histo;
}

TH2F *hardlimit(TH2F *limit) {
  TH2F *co = (TH2F*)limit->Clone("hcopy");
  for(int i=1;i<limit->GetNbinsX();i++) {
    for(int j=1;j<limit->GetNbinsY();j++) {
       if(limit->GetBinContent(i,j)<1&&limit->GetBinContent(i,j)>0) co->SetBinContent(i,j,1);
       else co->SetBinContent(i,j,0);
    }
  }
  return co;
}

TH2F* make_exclusion_shape(TH2F *excl, int isprimary) {
  TH2F *exclusion = (TH2F*)excl->Clone("exclusion");
  for(int i=1; i<(excl->GetNbinsX()+1); i++) {
    for(int j=1; j<(excl->GetNbinsY()+1); j++) {
      if(excl->GetBinContent(i,j)<1&&excl->GetBinContent(i,j)>0) exclusion->SetBinContent(i,j,0.01);
      else exclusion->SetBinContent(i,j,0);
    }
  }
  exclusion->SetLineColor(kBlue);
  exclusion->SetLineWidth(2);
  //exclusion->SetLineWidth(4); // paper style
  exclusion->SetLineStyle(isprimary);
  return exclusion;
}

void produce_extensive_plots(TH2F *xsec,TH2F *limits,TH2F *rellimits, TH2F *rellimitsd3, TH2F *rellimitst3, int scantype) {
  TCanvas *ca = new TCanvas("ca","ca",2400,1200);
  ca->Divide(4,2);
  ca->cd(1);
  ca->cd(1)->SetLogz(1);
  xsec->GetZaxis()->SetRangeUser(0.001,1000);
  xsec->Draw("COLZ");
  TText *title0 = write_title("Reference Cross Section");
  title0->Draw("same");
  ca->cd(2);
  ca->cd(2)->SetLogz(1);
  limits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
  limits->Draw("COLZ");
  TText *title = write_title("Cross Section Upper Limit");
  title->Draw("same");
  ca->cd(3);
  ca->cd(3)->SetLogz(1);
  TH2F *limit_ref = (TH2F*)limits->Clone("limit_ref");
  limit_ref->Divide(xsec);
  limit_ref->GetZaxis()->SetRangeUser(limits_ratio_lower_bound,limits_ratio_upper_bound);
  limit_ref->Draw("COLZ");
  TText *title2 = write_title("Cross Section UL / XS");
  title2->Draw("same");
  ca->cd(4);
  ca->cd(4)->SetLogz(1);
  limits->SetTitle("");
  limits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
  limits->SetZTitle("95% CL upper limit on #sigma [pb]");

//  limits->GetZaxis()->SetTitleOffset(1.1); // This is set in set_range
//  limits->GetZaxis()->CenterTitle(); // This is set in set_range
  limits->Draw("COLZ");

  TGraph *exclline = MarcosExclusionLine(rellimits, scantype);
  TGraph *thinexclline = thin_line(exclline);

  TGraph *excllinet3 = MarcosExclusionLine(rellimitst3, scantype);
  excllinet3->SetLineStyle(2);
  TGraph *thinexcllinet3 = thin_line(excllinet3);

  TGraph *excllined3 = MarcosExclusionLine(rellimitsd3, scantype);
  excllined3->SetLineStyle(3);
  TGraph *thinexcllined3 = thin_line(excllined3);

  ca->cd(4);
  exclline->Draw();
  thinexclline->Draw();
  excllinet3->Draw();
  thinexcllinet3->Draw();
  excllined3->Draw();
  thinexcllined3->Draw();
  stringstream partial;
  partial << "Limits/exclusion__" << limits->GetName();
  fill_with_text(exclline,excllined3,excllinet3,ca->cd(4),scantype);
  CompleteSave(ca,partial.str());

  ca->cd(5);
  (hardlimit(rellimits))->Draw("COL");
  exclline->Draw("same");
  ca->cd(6);
  (hardlimit(rellimitst3))->Draw("COL");
  excllinet3->Draw("same");
  ca->cd(7);
  (hardlimit(rellimitsd3))->Draw("COLZ");
  excllined3->Draw("same");
  
  CompleteSave(ca,partial.str()+"__PlusInfo");
  delete ca;
}

bool fail(double xs, double xsLimit) {return xsLimit>1 or !xsLimit;}

void get_Marias_exclusion_line(TH2F *limit_ref, float pointsX[200], float pointsY[200], int &ixNew, int &counter, int &foundDiag) {
// part of Mariarosaria's getRefXsecGraph function
  bool enough = false;
  Int_t  nBinX= limit_ref->GetXaxis()->GetNbins();
  Int_t  nBinY= limit_ref->GetYaxis()->GetNbins();
  for(int i=1; i<(nBinX+1); i++) {
    for(int j=1; j<(nBinY+1); j++) {
      if(limit_ref->GetBinContent(i,j-1)==0) continue;
      if(limit_ref->GetBinContent(i,j+1)==0) continue;
      if(limit_ref->GetBinContent(i,j)==0) continue;
      if( limit_ref->GetBinContent(i,j)>0. && limit_ref->GetBinContent(i,j)<=1.) {
        double xsLimitAbove = limit_ref->GetBinContent(i, j+1);

        double xsLimit = limit_ref->GetBinContent(i, j);
        if((fail(1.0,xsLimitAbove)) && (!fail(1.0,xsLimit)) ) {
          pointsX[counter]=limit_ref->GetXaxis()->GetBinCenter(i);
          pointsY[counter]=limit_ref->GetYaxis()->GetBinCenter(j);
          counter++;
          enough = !xsLimitAbove;          
          if(!enough && !foundDiag) foundDiag=counter;
//                cout << "enough " << enough << " foundDiag " << foundDiag << endl;
        }
      }
    }
  }
  
}

TH2F* flipth2f(TH2F *histo) {
   if(histo->GetNbinsX()!=histo->GetNbinsY()) {
        write_error(__FUNCTION__,"number of bins is not the same for x & y !!!");
        return NULL;
   }
   TH2F *rethisto = (TH2F*)histo->Clone();
   for(int i=1;i<histo->GetNbinsX();i++) {
        for(int j=1;j<histo->GetNbinsY();j++) {
                rethisto->SetBinContent(j,i,histo->GetBinContent(i,j));
        }
   }
   return rethisto;
}

TGraph* get_exclusion_line(TH2F *limit_ref) {
  
 float pointsX[200];
 float pointsY[200];
 int ixNew=0;
 int counter=0;
 int foundDiag=0;

 TH2F *fakehisto = flipth2f(limit_ref);
 get_Marias_exclusion_line(fakehisto,pointsY,pointsX,ixNew,counter,foundDiag);// x&y deliberately switched!
  if(counter>1) {
    pointsX[counter]=pointsX[counter-1];
    pointsY[counter]=50;  
  }

  const int newCounter=counter;
  Double_t newPointsX[newCounter+1];
  Double_t newPointsY[newCounter+1];
  float lastx=-100;
  float lasty=-100;
  
  
  for (int ix = 0; ix < counter+1; ++ix) {
    if(ix<(foundDiag-2)) continue;
    if(ix!=counter && pointsX[ix+1]==pointsX[ix]) continue;
    if(pointsX[ix]<PlottingSetup::mglustart) continue;
    if(lasty>-100&&pointsY[ix]<51) continue;
//    cout << pointsX[ix] << " , " << pointsY[ix] << endl;
   newPointsX[ixNew]=pointsX[ix]; 
   newPointsY[ixNew]=pointsY[ix];
   ixNew++;
   lastx=pointsX[ix];
   lasty=pointsY[ix];
  }
  string titleHisto="tester";
//  sprintf(titleHisto,"graph_%s_%f",limit_ref->GetName(),refMult);
//  cout << "HERE " << titleHisto << endl;

 TGraph * gr = new TGraph(ixNew,newPointsX,newPointsY);
 gr->SetName(titleHisto.c_str());
 gr->Draw("same");
 return gr;
}

TGraph *MarcosExclusionLine(TH2F *exclusionshape, int scantype) {
//  write_warning(__FUNCTION__,"USING MARIAS ALGORITHM...");
//  return get_exclusion_line(exclusionshape);
  
  TH2F *fakehisto = flipth2f(exclusionshape);
  TGraph *fakegraph = get_mSUGRA_exclusion_line(fakehisto, scantype);
  TGraph *realgraph = new TGraph(fakegraph->GetN());
  double x,y;
  float last_x=0;
  float last_y=0;
  int counter=0;
  for(int i=0;i<fakegraph->GetN();i++) {
        fakegraph->GetPoint(i,x,y);
        if(scantype==PlottingSetup::SMS) {
          if(y-x<75) {
            realgraph->SetPoint(counter,last_x,last_y);
            counter++;
            continue;
          }
        }
        realgraph->SetPoint(counter,y,x);
        last_x=y;
        last_y=x;
        counter++;
  }
  realgraph->SetLineColor(TColor::GetColor("#151515")); //nice black
  realgraph->SetLineWidth(2);
  //realgraph->SetLineWidth(4);//paper style

  return realgraph;
}
  
TGraph* thin_line(TGraph *gr) {
  TGraph *thin = (TGraph*) gr->Clone(concatenate("thin",gr->GetTitle()));
  thin->SetLineWidth(1);
  thin->SetLineColor(kWhite);
  return thin;
}

void make_SMS_exclusion(TH2F *rawlimits,TH2F *xsec,int scantype,std::string& scanx) {
  TH2F *limits = prep_histo(rawlimits,scantype); // this is to be independent of the style used at creation time
  //here we get some limits and the cross section; we want to make an exclusion plot!
  TH2F *rellimits = (TH2F*)limits->Clone("rellimits");
  TH2F *rellimitsd3 = (TH2F*)limits->Clone("rellimitsd3"); // one third the cross section ("divided by 3" -> d3)
  TH2F *rellimitst3 = (TH2F*)limits->Clone("rellimitst3"); // three times the cross section ("times 3" -> t3)
  
  if(!xsec ) {
    cout << "Watch out, cross section map is invalid!" << endl;
    return;
  }
  
  rellimits->Divide(xsec);
  
  
  for(int i=1;i<=rellimits->GetNbinsX();i++) {
    for(int j=1;j<=rellimits->GetNbinsY();j++) {
      rellimitst3->SetBinContent(i,j,(rellimits->GetBinContent(i,j))/3.0);
      rellimitsd3->SetBinContent(i,j,(rellimits->GetBinContent(i,j))*3.0);
    }
  }
  

//  TH2F *exclusionshape = make_exclusion_shape(rellimits,1);
//  TH2F *exclusionshapet3 = make_exclusion_shape(rellimitst3,2);
//  TH2F *exclusionshaped3 = make_exclusion_shape(rellimitsd3,3);
  
  //Now let's produce the plots!
  
  set_range(xsec,scantype,false);
  set_range(limits,scantype,false);
  
  bool drawdoubleline=false; //draw nice thin line on top of thick outer line
  TGraph *exclline = MarcosExclusionLine(rellimits, scantype);
  TGraph *thinexcline = thin_line(exclline);

  TGraph *excllinet3 = MarcosExclusionLine(rellimitst3, scantype);
  excllinet3->SetLineStyle(2);
  TGraph *thinexclinet3 = thin_line(excllinet3);

  TGraph *excllined3 = MarcosExclusionLine(rellimitsd3, scantype);
  excllined3->SetLineStyle(3);
  TGraph *thinexclined3 = thin_line(excllined3);

  produce_extensive_plots(xsec,limits,rellimits, rellimitst3, rellimitsd3,scantype);
  
  TCanvas *finalcanvas = new TCanvas("finalcanvas","finalcanvas");
  finalcanvas->SetLogz(1);
  finalcanvas->cd();
  limits->Draw("COLZ");


  TLine *desertline;
  if(drawefficiencydesertline) {
        desertline = new TLine(375,50,1200,875);
        desertline->SetLineWidth(3);
        //desertline->SetLineWidth(4); // paper style
        desertline->SetLineColor(kBlack);
        desertline->Draw("same");
  }

  exclline->Draw("c");

//  fill_with_text(exclline,excllined3,excllinet3,finalcanvas,scantype,scanx);
  stringstream real;
  real << "Limits/final_exclusion__" << limits->GetName();
  
  if(Contains(limits->GetName(),"bestlimits")) {
    cout << "----------> " << limits->GetName() << endl;
    TFile *f = new TFile("limits.root","RECREATE");
    thinexcline->SetName("ExclusionLine");
    limits->SetName("UpperLimits");
    limits->Write();
    thinexcline->Write();
    f->Close();
  }
  
  
  exclline->Draw("l");
  if(drawdoubleline) thinexcline->Draw("");
  excllinet3->Draw("");
  if(drawdoubleline) thinexclinet3->Draw("");
  excllined3->Draw("");
  if(drawdoubleline) thinexclined3->Draw("");
  CompleteSave(finalcanvas,real.str());

  delete finalcanvas;
}

TH1* getHisto(char * filename, char* histoName, char * dirName, int nBin, double lumi)
{
  TH1 * hpt_=0;
  TFile *file0 = TFile::Open(filename);
  if(!file0) return hpt_;
  TDirectory *dir;
  TH2D * hMuPt;

  if(dirName == "0") {
    hMuPt = (TH2D*) file0->Get(histoName);
  } else {
    dir = (TDirectory*) file0->Get(dirName);
    if(!dir) return hpt_;
    hMuPt = (TH2D*) dir->Get(histoName);
  }

  if(hMuPt) {
  hpt_ = (TH1*) hMuPt->Clone();
  hpt_->Sumw2();
  hpt_->Scale(1./lumi); // this take into into account the luminosity
  hpt_->SetLineWidth(2);
  hpt_->Rebin(nBin);
  double nBinX=hpt_->GetNbinsX();

  // overFlow
  hpt_->SetBinContent((int)nBinX,hpt_->GetBinContent((int)nBinX)+hpt_->GetBinContent((int)nBinX+1));
  hpt_->SetDirectory(0);
  file0->Close();
  hpt_->SetLineWidth(3);
  }
  return hpt_;
}

pair <float,float> find_point(TH2F* histo, int xbin, bool mSUGRA=true) {
  TH1F *flathisto;
  if(mSUGRA) flathisto = new TH1F("flat","flat",histo->GetNbinsY(),histo->GetYaxis()->GetBinLowEdge(1),histo->GetYaxis()->GetBinLowEdge(histo->GetNbinsY())+histo->GetYaxis()->GetBinWidth(histo->GetNbinsY()));
  else flathisto = new TH1F("flat","flat",histo->GetNbinsX(),histo->GetXaxis()->GetBinLowEdge(1),histo->GetXaxis()->GetBinLowEdge(histo->GetNbinsX())+histo->GetXaxis()->GetBinWidth(histo->GetNbinsX()));
  
  int nbins=histo->GetNbinsY();
  if(!mSUGRA) nbins=histo->GetNbinsX();
  for(int i=1;i<nbins;i++) {
    float value=(histo->GetBinContent(xbin,i));
    if(value<20&&value>0) flathisto->SetBinContent(i,value);
  }
  
  float pointone=-100;
  nbins=flathisto->GetNbinsX();
  if(!mSUGRA) nbins=flathisto->GetNbinsY();
  for(int i=nbins;i>=1;i--) {
    if(pointone<0&&flathisto->GetBinContent(i)<1&&flathisto->GetBinContent(i)>0) pointone=flathisto->GetBinLowEdge(i)+flathisto->GetBinWidth(i);
  }
  pair <float,float> anything;
  if(mSUGRA) anything.first=histo->GetXaxis()->GetBinCenter(xbin);
  else anything.first=histo->GetYaxis()->GetBinCenter(xbin);
  anything.second=pointone;
  delete flathisto;
  return anything;
}

pair <float,float> find_interpolated_point(TH2F* histo, int xbin, bool mSUGRA=true) {
  
  int minaccept=4;
  TCanvas *flatcan = new TCanvas("flatcan","flatcan");
  stringstream histoname;
  histoname << "exclusion shape for x bin " << xbin;
  TH1F *flathisto;
  if(mSUGRA) flathisto = new TH1F("flat",histoname.str().c_str(),histo->GetNbinsY(),histo->GetYaxis()->GetBinLowEdge(1),histo->GetYaxis()->GetBinLowEdge(histo->GetNbinsY())+histo->GetYaxis()->GetBinWidth(histo->GetNbinsY()));
  else flathisto = new TH1F("flat",histoname.str().c_str(),histo->GetNbinsX(),histo->GetXaxis()->GetBinLowEdge(1),histo->GetXaxis()->GetBinLowEdge(histo->GetNbinsX())+histo->GetXaxis()->GetBinWidth(histo->GetNbinsX()));
  
  int acceptedpoints=0;
  int nbins=histo->GetNbinsY();
  if(!mSUGRA) histo->GetNbinsX();
  for(int i=1;i<nbins;i++) {
    float value=0;
    if(i<=nbins-2) value=((1/3.0)*(histo->GetBinContent(xbin,i)+histo->GetBinContent(xbin+1,i)+histo->GetBinContent(xbin+2,i)));
    if(i==nbins-1) value=((1/2.0)*(histo->GetBinContent(xbin,i)+histo->GetBinContent(xbin,i+1)));
    if(i==nbins) value=(histo->GetBinContent(xbin,i));
    if(value<20&&value>0) {
      flathisto->SetBinContent(i,value);
      flathisto->SetBinError(i,TMath::Sqrt(value));
      acceptedpoints++;
    }
  }
  
  float pointone=-100;
  TLine *excluded;
  if(acceptedpoints>minaccept) {
    flathisto->Fit("expo","Q");
    TF1 *fitfunc = (TF1*)flathisto->GetFunction("expo");
    pointone=-(fitfunc->GetParameter(0)/fitfunc->GetParameter(1));
    excluded=new TLine(pointone,0,pointone,10);
  }
  
  pair <float,float> anything;
  anything.first=histo->GetXaxis()->GetBinCenter(xbin);
  anything.second=pointone;
  stringstream flatname;
  flathisto->GetYaxis()->SetRangeUser(0,10);
  flathisto->Draw();
  if(acceptedpoints>minaccept) excluded->SetLineColor(kGreen);
  if(acceptedpoints>minaccept) excluded->SetLineStyle(2);
  if(acceptedpoints>minaccept) excluded->SetLineWidth(2);
  if(acceptedpoints>minaccept) excluded->Draw("same");
  flatname << "Limits/partials/partial_" << xbin << "___" << histo->GetName() << ".png";
  if(draweachone) CompleteSave(flatcan,flatname.str());
  delete flatcan;
  delete flathisto;
  return anything;
}

void project_to_last_point_on_line(float x2, float y2, float x1, float y1, float &x, float &y, int scantype) {
// original code: continue the same way to the diagonal (i.e. conserve the angle to the diagonal)
/*
  float deltax=x2-x1;
  float deltay=y2-y1;
  float b = (y1-x1)/(deltax-deltay);
  if(scantype==PlottingSetup::SMS) b = (y1-75.-x1)/(deltax-deltay);
  x = x1 + b * deltax;
  y = y1 + b * deltay;
  if(x<0||y<0) x=-1;
  if((scantype==PlottingSetup::GMSB||scantype==PlottingSetup::SMS) && (x>PlottingSetup::mgluend||y>PlottingSetup::mLSPend)) x=-1;
*/
//new code: just connect directly with the diagonal
float b = 0; 
if(scantype==PlottingSetup::SMS) b = -75;
//float alpha = (b+x2-y2)/2.0;
x           = floor((x2+y2-b)/2.0);
y           = floor(x + b);
}

TGraph* get_mSUGRA_exclusion_line(TH2F *exclusionhisto, int scantype) {
  exclusionhisto->SetStats(0);
  exclusionhisto->GetZaxis()->SetRangeUser(0,2);

  vector<pair <float,float> > points;
  
  for(int i=1;i<exclusionhisto->GetNbinsX();i++) {
    pair<float,float> anything = find_point(exclusionhisto,i);
    pair<float,float> intthing = find_interpolated_point(exclusionhisto,i);
    if(intthing.second>anything.second) anything.second=intthing.second;
    if(anything.second>100&&anything.second<1500) points.push_back(anything);
  }
  
  Double_t xpoints[points.size()];
  Double_t spoints[points.size()];
  
  TGraph *graph = new TGraph(points.size()+1);
  graph->SetLineColor(kBlack);
  graph->SetLineWidth(2);
  
  int pointcounter=0;
  float lastx=0;
  float lasty=0;
  float lastx2=0;
  float lasty2=0;
  
  for(int i=0;i<(int)points.size();i++) {
    xpoints[i]=points[i].first;
    spoints[i]=points[i].second;
    if(scantype==PlottingSetup::mSUGRA) {
        if(i>1&&i<(int)points.size()-1) spoints[i]=(1.0/3.0)*(points[i-1].second+points[i].second+points[i+1].second);
        if(i>2&&i<(int)points.size()-2) spoints[i]=(1.0/5.0)*(points[i-2].second+points[i-1].second+points[i].second+points[i+1].second+points[i+2].second);
        if(i>3&&i<(int)points.size()-3) spoints[i]=(1.0/7.0)*(points[i-3].second+points[i-2].second+points[i-1].second+points[i].second+points[i+1].second+points[i+2].second+points[i+3].second);
    }
 
    bool usethispoint=true;
    if(scantype==PlottingSetup::GMSB) if(spoints[i]<points[i].first) usethispoint=false;
    if(scantype==PlottingSetup::SMS)  if(spoints[i]-75.<points[i].first) usethispoint=false;

    if((scantype==PlottingSetup::GMSB||scantype==PlottingSetup::SMS) && (spoints[i]>PlottingSetup::mgluend||points[i].first>PlottingSetup::mLSPend)) usethispoint=false; 

    if(usethispoint) {
        graph->SetPoint(pointcounter,points[i].first,spoints[i]);
        lastx2=lastx;
        lasty2=lasty;
        lastx=points[i].first;
        lasty=spoints[i];
        pointcounter++;
    }
  }
  for(int i=pointcounter;i<=(int)points.size();i++) graph->SetPoint(i,lastx,lasty);
  if(scantype==PlottingSetup::GMSB||scantype==PlottingSetup::SMS) {
        //The final point will be a continuation of the last one, towards the diagonal
        float x,y;
        project_to_last_point_on_line(lastx,lasty,lastx2,lasty2,x,y,scantype);
        if(x>0&&y>0) graph->SetPoint(points.size(),y,x);
        else graph->SetPoint(points.size(),lastx,lasty);
  }

  return graph;
}

TH2F* cast_into_shape(TH2F *origin, TH2F *reference) {
  TH2F *newh = (TH2F*)reference->Clone(origin->GetName());
  for(int ix=1;ix<=reference->GetNbinsX();ix++) {
    for(int iy=1;iy<=reference->GetNbinsY();iy++) {
//      reference->SetBinContent(ix,iy,origin->GetBinContent(ix,iy));
      newh->SetBinContent(ix,iy,origin->GetBinContent(ix,iy));
    }
  }
  return newh;
}

void draw_mSUGRA_exclusion(TH2F *ocrosssection, TH2F *oFilterEfficiency, TH2F *oabsXS, TH2F *limitmap, TH2F *expmap, TH2F *expplusmap, TH2F *expminusmap, TH2F *exp2plusmap, TH2F *exp2minusmap) {
  TH2F *crosssection = (TH2F*)ocrosssection->Clone("crosssection");
//  TH2F *limitmap = (TH2F*)olimitmap->Clone(((string)olimitmap->GetName()+"clone").c_str());
  TH2F *cleanhisto = (TH2F*)limitmap->Clone("clean");
  for(int ix=1;ix<=cleanhisto->GetNbinsX();ix++) {
    for(int iy=1;iy<=cleanhisto->GetNbinsY();iy++) {
      cleanhisto->SetBinContent(ix,iy,0);
    }
  }
  
  
  TH2F *FilterEfficiency;
  TH2F *absXS;

  
  write_warning(__FUNCTION__,"You'll want to switch off 'wrongwaytodothis')");
  
  if(wrongwaytodothis) {
    //this part is the one you want to remove.
    TFile *Efficiencies = new TFile("FilterEfficiencyv3.root");
    FilterEfficiency = cast_into_shape((TH2F*) Efficiencies->Get("FilterEfficiency"),limitmap);
    assert(FilterEfficiency);
    assert(crosssection);
    absXS=(TH2F*)crosssection->Clone("absXS");
    crosssection->Multiply(FilterEfficiency);
  } else {
    //this part is the one you want to keep!
    FilterEfficiency=(TH2F*)oFilterEfficiency->Clone("FilterEfficiency");
    absXS=(TH2F*)oabsXS->Clone("absXS");
  }
    
  
  TH2F *limits = (TH2F*)limitmap->Clone("limits");
  set_range(limits,true,false);
  limitmap->Divide(crosssection);
  expminusmap->Divide(crosssection);
  expplusmap->Divide(crosssection);
  exp2minusmap->Divide(crosssection);
  exp2plusmap->Divide(crosssection);
  expmap->Divide(crosssection);
  TGraph *observed = get_mSUGRA_exclusion_line(limitmap, PlottingSetup::mSUGRA);
  observed->SetLineColor(kRed);
  TGraph *expminus = get_mSUGRA_exclusion_line(expminusmap, PlottingSetup::mSUGRA);
  TGraph *expplus  = get_mSUGRA_exclusion_line(expplusmap, PlottingSetup::mSUGRA);
  TGraph *exp2minus;
  if(draw2sigma) exp2minus = get_mSUGRA_exclusion_line(exp2minusmap, PlottingSetup::mSUGRA);
  TGraph *exp2plus;
  if(draw2sigma) exp2plus = get_mSUGRA_exclusion_line(exp2plusmap, PlottingSetup::mSUGRA);

  TGraph *expected = new TGraph(expminus->GetN()+expplus->GetN());
  TGraph *expected2; 
  if(draw2sigma) expected2 = new TGraph(exp2minus->GetN()+exp2plus->GetN());
  for(int i=0;i<=expminus->GetN();i++) {
    Double_t x,y;
    expminus->GetPoint(i,x,y);
    expected->SetPoint(i,x,y);
  }
    
  for(int i=0;i<=exp2minus->GetN();i++) {
    Double_t x,y;
    exp2minus->GetPoint(i,x,y);
    expected2->SetPoint(i,x,y);
  }
  for(int i=exp2plus->GetN()-1;i>=0;i--) {
    Double_t x,y;
    exp2plus->GetPoint(i,x,y);
    expected2->SetPoint(exp2minus->GetN()+(exp2plus->GetN()-i),x,y);
  }
  for(int i=expplus->GetN()-1;i>=0;i--) {
    Double_t x,y;
    expplus->GetPoint(i,x,y);
    expected->SetPoint(expminus->GetN()+(expplus->GetN()-i),x,y);
  }
  expected->SetFillColor(TColor::GetColor("#9FF781"));
  if(draw2sigma) expected2->SetFillColor(TColor::GetColor("#F3F781"));
  
  smooth_line(observed);
  smooth_line(expected);
  if(draw2sigma) smooth_line(expected2);
  
  TCanvas *te = new TCanvas("te","te");
  te->SetRightMargin(standardmargin);
//  decorate_mSUGRA(cleanhisto,te,expected,expected2,observed);
  TH2F *noh = new TH2F("noh","noh",1,1,2,1,1,2);
  SugarCoatThis(te,10,noh,observed);
//  expected->Draw("c");
//  observed->Draw("c");
  stringstream saveas;
  if((int)((string)limitmap->GetName()).find("limitmap")>0) saveas << "Limits/final_exclusion_for_JZB_geq_" << ((string)limitmap->GetName()).substr(((string)limitmap->GetName()).find("limitmap")+8,10);
  else saveas << "Limits/final_exclusion_for_bestlimits";
  CompleteSave(te,saveas.str());
  delete te;
  
  TCanvas *overview = new TCanvas("overview","overview",1500,1000);
  
  set_range(crosssection,true,false);
  set_range(limits,true,false);
  set_range(limitmap,true,false);
  
  overview->Divide(3,2);
  overview->cd(1);
  overview->cd(1)->SetLogz(1);
  absXS->GetZaxis()->SetRangeUser(0.0001,100);
  absXS->Draw("COLZ");
  TText *title0 = write_title("Cross Section");
  title0->Draw("same");
  overview->cd(2);
  FilterEfficiency->GetZaxis()->SetRangeUser(0.01,0.7);
  FilterEfficiency->Draw("COLZ");
  TText *title0aa = write_title("Filter #epsilon");
  title0aa->Draw("same");
  overview->cd(3);
  overview->cd(3)->SetLogz(1);
  crosssection->GetZaxis()->SetRangeUser(0.0001,100);
  crosssection->Draw("COLZ");
  TText *title0a = write_title("Filter #epsilon x Cross Section");
  title0a->Draw("same");
  
  overview->cd(4);
  overview->cd(4)->SetLogz(1);
  limits->GetZaxis()->SetRangeUser(0.01,100);
  limits->Draw("COLZ");
  TText *title1 = write_title("Cross Section Upper Limit");
  title1->Draw("same");
  overview->cd(5);
  limitmap->Draw("COLZ");
  TText *title2 = write_title("UL/XS");
  title2->Draw("same");
  observed->Draw("c");
  overview->cd(6);
  overview->cd(6)->SetRightMargin(standardmargin);
//  decorate_mSUGRA(cleanhisto,overview->cd(4),expected,expected2,observed);
  SugarCoatThis(overview->cd(6),10,noh,observed);
//  observed->Draw("c");
  stringstream saveas2;
  if((int)((string)limitmap->GetName()).find("limitmap")>0) saveas2 << "Limits/exclusion_overview_for_JZB_geq_" << ((string)limitmap->GetName()).substr(((string)limitmap->GetName()).find("limitmap")+8,10);
  else saveas2 << "Limits/exclusion_overview_for_bestlimits";
  CompleteSave(overview,saveas2.str());
  delete overview;
  delete noh;
  delete crosssection;
  delete absXS;
  delete FilterEfficiency;
  
}

TH2F* get_XS(string filename, string mass, TH2F *histo) {///DONE
  TH1 *hRef=getHisto((char*) filename.c_str(), (char*) mass.c_str(), (char*) "0", 1,1);
  TH2F * xsec= new TH2F("xsec","xsec",60, 0., 1500., 60, 0., 1500.);
  for(int i=1; i<(histo->GetNbinsX()+1); i++) {
    for(int j=1; j<(histo->GetNbinsY()+1); j++) {
      if(!(histo->GetBinContent(i,j)==0)) xsec->SetBinContent(i,j,hRef->GetBinContent(hRef->FindBin(xsec->GetXaxis()->GetBinCenter(i))));       
    }
  }
  return xsec;
}

string give_nice_source_label(string name) {
  int mappoint=name.find("map");
  if(mappoint<0||mappoint>500) return name; // this mean that something weird is happening
  stringstream nice_label;
  nice_label << "JZB > " << name.substr(mappoint+3,name.size()) << " GeV";
  return nice_label.str();
}

Int_t get_exclusion_region_color(double value, TH2F *histo) {
   Double_t wmin  = histo->GetMinimum();
   Double_t wmax  = histo->GetMaximum();
   Double_t wlmin = wmin;
   Double_t wlmax = wmax;


   Int_t ncolors = gStyle->GetNumberOfColors();
   Int_t ndivz = histo->GetContour();
   if (ndivz == 0) return 0;
   ndivz = TMath::Abs(ndivz);
   Int_t theColor, color;
   Double_t scale = ndivz / (wlmax - wlmin);

   if (value < wlmin) value = wlmin;

   color = Int_t(0.01 + (value - wlmin) * scale);

   theColor = Int_t((color + 0.99) * Double_t(ncolors) / Double_t(ndivz));
   return gStyle->GetColorPalette(theColor);
}  

float findmaxentry(vector<TH2F*> histos, int i, int j) {
  float max = histos[0]->GetBinContent(i,j);
  for(int k=0;k<(int)histos.size();k++) {
    float entry = histos[k]->GetBinContent(i,j);
    if(entry>=max) max=entry;
  }

  return max;
}

TH2F *make_best_limits(vector<TH2F*> explimits,vector<TH2F*> obslimits, int scantype, std::string& scanx, vector<TH2F*> exp1mlimits, vector<TH2F*> exp1plimits, vector<TH2F*> exp2mlimits, vector<TH2F*> exp2plimits, vector<TH2F*> &allbestexplimits) {
  if(obslimits.size()==0) {
    write_warning(__FUNCTION__,"There are no observed limits! Cannot continue!");
    TH2F *err;
    return err;
  }
  if(explimits.size()==0) {
    write_warning(__FUNCTION__,"There are no expected limits! Will compute best limits based on observed limits! (WATCH OUT THIS IS DISCOURAGED!");
    for(int i=0;i<(int)obslimits.size();i++) explimits.push_back(obslimits[i]);
  }
  TH2F *bestlimit=(TH2F*)obslimits[0]->Clone("bestlimits");
  TH2F *bestexplimit=(TH2F*)obslimits[0]->Clone("bestexplimits");
  TH2F *bestlimitsource=(TH2F*)obslimits[0]->Clone("bestlimitsource");
  TH2F *bestexp1plimit=(TH2F*)obslimits[0]->Clone("bestexp1plimit");
  TH2F *bestexp1mlimit=(TH2F*)obslimits[0]->Clone("bestexp1mlimit");
  TH2F *bestexp2plimit=(TH2F*)obslimits[0]->Clone("bestexp2plimit");
  TH2F *bestexp2mlimit=(TH2F*)obslimits[0]->Clone("bestexp2mlimit");
  
  for(int i=1;i<=explimits[0]->GetNbinsX();i++) {
    for(int j=1;j<=explimits[0]->GetNbinsY();j++) {
      float min=findmaxentry(explimits,i,j);
      float omin=obslimits[0]->GetBinContent(i,j);
      int source=0;
      for(int k=0;k<(int)explimits.size();k++) {
        float currlim=explimits[k]->GetBinContent(i,j);
        if(currlim<=min&&currlim>0) {
          min=currlim;
          omin=obslimits[k]->GetBinContent(i,j);
          source=k+1;
        }
      }
      if(min>0) {
        bestlimitsource->SetBinContent(i,j,source);
        bestlimit->SetBinContent(i,j,omin);
        bestexplimit->SetBinContent(i,j,min);
        if(scantype==PlottingSetup::mSUGRA) {
          bestexp1plimit->SetBinContent(i,j,exp1plimits[source-1]->GetBinContent(i,j));
          bestexp1mlimit->SetBinContent(i,j,exp1mlimits[source-1]->GetBinContent(i,j));
          bestexp2plimit->SetBinContent(i,j,exp2plimits[source-1]->GetBinContent(i,j));
          bestexp2mlimit->SetBinContent(i,j,exp2mlimits[source-1]->GetBinContent(i,j));
        }
      }
    }
  }
  gStyle->SetPadRightMargin(standardmargin);
  gStyle->SetPadLeftMargin(0.18);

  TCanvas *canlimsource = new TCanvas("limsource","Source of best limits");
  set_range(bestlimitsource,scantype,false);
  bestlimitsource->SetTitle("Source of limit for best limits");
  bestlimitsource->GetZaxis()->SetRangeUser(0,explimits.size()+1);
  bestlimitsource->GetXaxis()->CenterTitle(0);
  bestlimitsource->GetYaxis()->CenterTitle(0);
  bestlimitsource->Draw("COL");
  gPad->Update();
  TLegend *sourceleg = new TLegend(0.2,0.5,0.55,0.75);
  for(int i=0;i<(int)explimits.size();i++) {
    stringstream legendentry;
    legendentry << i+1 << " = " << give_nice_source_label(explimits[i]->GetName());
    Int_t currcol=get_exclusion_region_color(i+1,bestlimitsource);
    explimits[i]->SetFillColor(currcol);
    explimits[i]->SetLineColor(currcol);
    sourceleg->AddEntry(explimits[i], legendentry.str().c_str(),"f");
  }
  sourceleg->SetLineColor(kWhite);sourceleg->SetFillColor(kWhite);
  sourceleg->SetTextSize(0.03);
  sourceleg->Draw("same");
  DrawPrelim();
  write_SMS_text(scantype,scanx);

  CompleteSave(canlimsource,"Limits/SourceOfBestLimits");
  gStyle->SetPadRightMargin(indentedmargin);
  allbestexplimits.push_back(bestexplimit);
  allbestexplimits.push_back(bestexp1plimit);
  allbestexplimits.push_back(bestexp1mlimit);
  allbestexplimits.push_back(bestexp2plimit);
  allbestexplimits.push_back(bestexp2mlimit);
  
  delete canlimsource;
  delete bestlimitsource;
  return bestlimit;
}


void create_exclusion_plots(vector<TH2F*> limits, int scantype, std::string& scanx) {
  TFile *xsecfile;
  if(scantype!=PlottingSetup::mSUGRA) {
    xsecfile = new TFile(xsecfilename.c_str());
    if(xsecfile->IsZombie()&&(scantype!=PlottingSetup::mSUGRA)) {
      write_error(__FUNCTION__,"Cross section file is invalid!!!!");
      return;
    }
    xsecfile->Close();
  }
  if(scantype!=PlottingSetup::mSUGRA) for(int i=0;i<(int)limits.size();i++) limits[i]->Scale(1./0.19); // because for T5zz we forced one z to decay leptonically
    
  vector<TH2F*> explimits;
  vector<TH2F*> exp1mlimits;
  vector<TH2F*> exp1plimits;
  vector<TH2F*> exp2mlimits;
  vector<TH2F*> exp2plimits;
  vector<TH2F*> obslimits;
  vector<TH2F*> flipmaps;
  vector<TH2F*> crosssections;
  vector<TH2F*> AbsCrossSection;
  vector<TH2F*> FilterEfficiencies;
  for(int ilim=0;ilim<(int)limits.size();ilim++) {
    if(TString(limits[ilim]->GetName()).Contains("_explimitmap")) explimits.push_back(limits[ilim]);
    if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_exp1plimitmap")) exp1plimits.push_back(limits[ilim]);
    if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_exp1mlimitmap")) exp1mlimits.push_back(limits[ilim]);
    if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_exp2plimitmap")) exp2plimits.push_back(limits[ilim]);
    if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_exp2mlimitmap")) exp2mlimits.push_back(limits[ilim]);
//    if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_XS")) crosssections.push_back(limits[ilim]);
    if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_absXS")) AbsCrossSection.push_back(limits[ilim]);
    if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_FilterEfficiency")) FilterEfficiencies.push_back(limits[ilim]);
    if(wrongwaytodothis) {
      if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_XS")) AbsCrossSection.push_back(limits[ilim]);
      if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_absXS")) crosssections  .push_back(limits[ilim]);
      if(scantype==PlottingSetup::mSUGRA && TString(limits[ilim]->GetName()).Contains("_XS")) FilterEfficiencies.push_back(limits[ilim]);
    }
    if(TString(limits[ilim]->GetName()).Contains("_limitmap")) obslimits.push_back(limits[ilim]);
//    if(TString(limits[ilim]->GetName()).Contains("_limitflipmap")) flipmaps.push_back(limits[ilim]);
  }
    
cout << "Size: " << AbsCrossSection.size() << endl;
  TH2F *xsec;
  if(scantype!=PlottingSetup::mSUGRA) xsec = adjust_histo(get_XS(xsecfilename,"gluino",limits[0]),limits[0]);
  vector<TH2F*> bestexplimits;
  TH2F *bestlimits = make_best_limits(explimits,obslimits,scantype, scanx, exp1mlimits, exp1plimits, exp2mlimits, exp2plimits, bestexplimits);
  bestlimits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);

  for(int ilim=0;ilim<(int)limits.size();ilim++) {
    limits[ilim]->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
  }
  
  if(scantype!=PlottingSetup::mSUGRA) {
    for(int ilim=0;ilim<(int)obslimits.size();ilim++) make_SMS_exclusion(obslimits[ilim],xsec,scantype,scanx);
    make_SMS_exclusion(bestlimits,xsec,scantype,scanx);
  } else {
    for(int ilim=0;ilim<(int)obslimits.size();ilim++) {
      draw_mSUGRA_exclusion(crosssections[0],FilterEfficiencies[0],AbsCrossSection[0],obslimits[ilim], explimits[ilim], exp1mlimits[ilim], exp1plimits[ilim], exp2mlimits[ilim], exp2plimits[ilim]);
    }
    draw_mSUGRA_exclusion(crosssections[0],FilterEfficiencies[0],AbsCrossSection[0],bestlimits, bestexplimits[0], bestexplimits[1], bestexplimits[2], bestexplimits[3], bestexplimits[4]);
  }
  delete bestlimits;
}


void smooth_line_once(TGraph *gr) {
  //going to smooth graph
  //need to take into account the DIRECTION we're heading so we noticed for expected shape when we turn around and go back
  float previousX=-100;
  int sign=1;
  for(int i=0;i<gr->GetN();i++) {
    Double_t x,y,x1,y1,x2,y2;
    gr->GetPoint(i,x,y);
    gr->GetPoint(i+1,x1,y1);//need to handle exception here
    gr->GetPoint(i-1,x2,y2);//need to handle exception here
    if(i!=gr->GetN()-1&& (sign*x<=sign*previousX||sign*x1<=sign*x)) {
      //turned around!
      sign=(-1)*sign;
      // do NOT do any smoothing on this point!
    } else {
      float newval=y;
      if(i>0&&i<(gr->GetN())-1) newval=(1.0/3.0)*(y+y1+y2);
      if(i==0) newval=(1.0/2.0)*(y+y1);
      if(i==gr->GetN()-1) newval=(1.0/2.0)*(y+y2);
      gr->SetPoint(i,x,newval);
    }
    previousX=x;
  }
}

void smooth_line(TGraph *gr) {
  smooth_line_once(gr);
  smooth_line_once(gr);
  smooth_line_once(gr);
}
void set_range(TH2F *histo, int scantype, bool pushoutyz=false) {
  gStyle->SetPadLeftMargin(0.18);
  gStyle->SetPadRightMargin(0.19);
  // histo->GetXaxis()->SetLabelSize(0.035); //paper style
  // histo->GetYaxis()->SetLabelSize(0.035); //paper style
  if(scantype==PlottingSetup::mSUGRA) {
    histo->GetXaxis()->SetRangeUser(0,PlottingSetup::m0end);
    histo->GetYaxis()->SetRangeUser(0,PlottingSetup::m12end);
    histo->GetXaxis()->SetTitle("m_{0} [GeV]");
    histo->GetYaxis()->SetTitle("m_{1/2} [GeV]");
    histo->GetXaxis()->SetRangeUser(0,PlottingSetup::m0end);
    histo->GetXaxis()->SetNdivisions(504,true);
  }
  if(scantype==PlottingSetup::SMS) {
    histo->GetXaxis()->SetRangeUser(50,1200);
    histo->GetYaxis()->SetRangeUser(50,1200);
    histo->GetXaxis()->SetTitle("gluino mass [GeV]");
    histo->GetYaxis()->SetTitle("#chi_{1}^{0} mass [GeV]");
    histo->GetXaxis()->SetTitleSize(0.055);
    histo->GetYaxis()->SetTitleSize(0.055);
    histo->GetZaxis()->SetTitleSize(0.055);
    histo->GetXaxis()->SetTitleOffset(1.0);
    histo->GetYaxis()->SetTitleOffset(1.4);
    histo->GetZaxis()->SetTitleOffset(1.2);
    if(pushoutyz) histo->GetZaxis()->SetTitleOffset(1.2);
  }
  if(scantype==PlottingSetup::GMSB) {
    histo->GetXaxis()->SetRangeUser(100,1200);
    histo->GetYaxis()->SetRangeUser(100,1200);
    histo->GetXaxis()->SetTitle("gluino mass [GeV]");
    histo->GetYaxis()->SetTitle("#chi_{1}^{0} mass [GeV]");
    histo->GetXaxis()->SetTitleSize(0.055);
    histo->GetYaxis()->SetTitleSize(0.055);
    histo->GetZaxis()->SetTitleSize(0.055);
    histo->GetXaxis()->SetTitleOffset(1.0);
    histo->GetYaxis()->SetTitleOffset(1.4);
    histo->GetZaxis()->SetTitleOffset(1.2); 
    if(pushoutyz) histo->GetZaxis()->SetTitleOffset(1.2);
  }

  //histo->GetXaxis()->CenterTitle();
  //histo->GetYaxis()->CenterTitle();
  histo->SetStats(0);
}

TH2F* adjust_histo(TH2F *oldhist, TH2F *refhisto) {
  TH2F *histo = new TH2F(refhisto->GetName(),refhisto->GetName(),
                         refhisto->GetNbinsX(),refhisto->GetXaxis()->GetBinLowEdge(1),refhisto->GetXaxis()->GetBinLowEdge(refhisto->GetNbinsX())+refhisto->GetXaxis()->GetBinWidth(refhisto->GetNbinsX()),
                         refhisto->GetNbinsY(),refhisto->GetYaxis()->GetBinLowEdge(1),refhisto->GetYaxis()->GetBinLowEdge(refhisto->GetNbinsY())+refhisto->GetYaxis()->GetBinWidth(refhisto->GetNbinsY()));
                         
  for(int ix=1;ix<=refhisto->GetNbinsX();ix++) {
    for(int iy=1;iy<=refhisto->GetNbinsX();iy++) {
      int binnum=refhisto->FindBin(refhisto->GetXaxis()->GetBinCenter(ix),refhisto->GetYaxis()->GetBinCenter(iy));
      histo->SetBinContent(binnum,oldhist->GetBinContent(ix,iy));
    }
  }
  return histo;
}

void process_syst_plot(TH2F *rhisto,string saveto,int scantype, std::string scanx = "") {
  TH2F *histo = prep_histo(rhisto,scantype); // this is to be independent of the style used at creation time
  float rightmargin=gStyle->GetPadRightMargin();
  gStyle->SetPadRightMargin(0.20);
  TString name = rhisto->GetName();
  if(name.Contains("Nevents")) gStyle->SetPadRightMargin(0.22);
  TCanvas *can = new TCanvas("syst_plot","Systematics Plot");
  set_range(histo,scantype,true);
  
  histo->GetZaxis()->CenterTitle();
  gStyle->SetStripDecimals(false);
  histo->GetXaxis()->SetDecimals(true);
  
  if(name.Contains("efficiency")) {
    histo->GetZaxis()->SetTitle("A #times #varepsilon (#geq 1 Z(ll))");
    histo->GetZaxis()->CenterTitle(0);
    //histo->GetZaxis()->SetRangeUser(0.0,0.15);
  }
  if(name.Contains("Nevents")) {
     histo->GetZaxis()->SetTitle("N(events)");
     histo->GetZaxis()->SetTitleOffset(histo->GetZaxis()->GetTitleOffset()+0.4);
  }
  if(name.Contains("sysjes")) {
    histo->GetZaxis()->SetTitle("Jet Energy Scale");
    histo->GetZaxis()->SetRangeUser(0.0,0.2);
  }
  if(name.Contains("sysjsu")) {
    histo->GetZaxis()->SetTitle("JZB Scale Uncertainty");
    histo->GetZaxis()->SetRangeUser(0.0,0.5);
  }
  if(name.Contains("sysresmap")) {
    histo->GetZaxis()->SetTitle("Resulution");
    histo->GetZaxis()->SetRangeUser(0.0,0.5);
  }
  if(name.Contains("sysstatmap")) {
    histo->GetZaxis()->SetTitle("Statistical Error");
    histo->GetZaxis()->SetRangeUser(0.0,0.01);
  }
  if(name.Contains("systotmap")) {
    histo->GetZaxis()->SetTitle("All Systematic Errors");
    histo->GetZaxis()->SetRangeUser(0.0,0.5);
  }
  
  
  histo->Draw("COLZ");
  DrawPrelim();
  float xpos_of_text = 0.22;
  if(name.Contains("_noscefficiencymap"))  {
    // Add some more decorations on these publication plots
    write_SMS_text( scantype, scanx, xpos_of_text );
    // Decode name and get JZB cut
    TPRegexp pat("\\d+$");
    size_t index = name.Index(pat,0);
    string cut = string("JZB > ")+(name(index,name.Length()-index).Data())+" GeV";
    //string cut = string("#splitline{JZB > ")+(name(index,name.Length()-index).Data())+" GeV}{n_{jets} #geq 3}"; //paper style
    TText *text = write_text(xpos_of_text,0.73,cut);
    text->SetTextAlign(11);
    text->SetTextSize(0.035);
    text->Draw();
    draw_diagonal_xchange( scantype, scanx );
  }
  
  CompleteSave(can,(saveto+(string)histo->GetName()));

  gStyle->SetPadRightMargin(rightmargin);
  
  delete can;
}

void make_all_syst_plots(vector<TH2F*> all_histos,int scantype,std::string& scanx) {
  string saveto="Systematics/";
  for(int iplot=0;iplot<(int)all_histos.size();iplot++) {
      process_syst_plot(all_histos[iplot],saveto,scantype,scanx);
  }
}

void process_file(TFile* file, float stdmargin) {
  standardmargin=stdmargin;
  xsecfilename=PlottingSetup::cbafbasedir+"/"+PlottingSetup::SMSReferenceXSFile;

  // can receive a file with systematics and limits mixed, or a file with systematics only , or a file with limits only. 
  TIter nextkey(file->GetListOfKeys());
  TKey *key;
  
  int scantype=PlottingSetup::SMS;
  std::string scanx = "0"; // Just for the legend
  if (TString(file->GetName()).Contains("T5zzl")) scanx = "0.75";
  else  if(TString(file->GetName()).Contains("T5zzh")) scanx = "0.25";
  else if(TString(file->GetName()).Contains("T5zz")) scanx = "0.5";
  
  vector<TH2F*> systematics_histos;
  vector<TH2F*> limits_histos;
  while ((key = (TKey*)nextkey()))
    {
      TObject *obj = key->ReadObj();
      if (!(obj->IsA()->InheritsFrom("TH2"))) continue;
      TString name=(TString)(obj->GetName());
      bool is_limit=false;
      bool is_systematic=false;
      
      if(name.Contains("limitmap")) is_limit=true;
      if(name.Contains("explimitmap")) is_limit=true;
      if(name.Contains("exp1plimitmap")) is_limit=true;
      if(name.Contains("exp2plimitmap")) is_limit=true;
      if(name.Contains("exp1mlimitmap")) is_limit=true;
      if(name.Contains("exp2mlimitmap")) is_limit=true;
      if(name.Contains("exclusionmap")) is_limit=true;
      if(name.Contains("crosssectionmap")) is_limit=true;
      if(name.Contains("XS")) is_limit=true;
      if(name.Contains("absXS")) is_limit=true;
      if(name.Contains("limitflipmap")) is_limit=true;

      if(name.Contains("sysjes")) is_systematic=true;
      if(name.Contains("sysjsu")) is_systematic=true;
      if(name.Contains("sysresmap")) is_systematic=true;
      if(name.Contains("efficiencymap")) is_systematic=true;
      if(name.Contains("efficiencymapAcc")) is_systematic=true;
      if(name.Contains("efficiencymapID")) is_systematic=true;
      if(name.Contains("efficiencymapJets")) is_systematic=true;
      if(name.Contains("efficiencymapSignal")) is_systematic=true;
      if(name.Contains("noscefficiencymap")) is_systematic=true;
      if(name.Contains("Neventsmap")) is_systematic=true;
      if(name.Contains("ipointmap")) is_systematic=true;
      if(name.Contains("syspdfmap")) is_systematic=true;
      if(name.Contains("systotmap")) is_systematic=true;
      if(name.Contains("sysstatmap")) is_systematic=true;
      if(name.Contains("timemap")) is_systematic=true;
      if(name.Contains("flipmap")) is_systematic=true;

      if(is_limit) limits_histos.push_back((TH2F*) obj);
      else if(is_systematic) systematics_histos.push_back((TH2F*) obj);
      if(name.Contains("mSUGRA")) scantype=PlottingSetup::mSUGRA;
      if(name.Contains("GMSB")) scantype=PlottingSetup::GMSB;
    }
    write_warning(__FUNCTION__,"Deactivated systematics plots");
//  if(systematics_histos.size()>0) make_all_syst_plots(systematics_histos,scantype,scanx);
  if(limits_histos.size()>0) create_exclusion_plots(limits_histos,scantype,scanx);
}

Revision:	1.9
Committed:	Mon Apr 30 08:39:44 2012 UTC (13 years ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.8:	+15 -25 lines
Log Message:	fixed some uint vs int comparisons; commented out/removed unused variables; made Wall happy