Plotting/Modules/ExclusionPlot.C

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

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

string xsecfilename;

void set_range(TH2F *histo, int scantype, bool pushoutlabels);
void smooth_line(TGraph *gr);
void decorate_mSUGRA(TH2F *cleanhisto, TVirtualPad *cvsSys,TGraph *expected,TGraph *expected2,TGraph *observed);
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);

void fill_with_text(TGraph *real, TGraph *down, TGraph *up, TVirtualPad *can, int scantype) {
  can->cd();
  TLegend* this_leg = new TLegend(0.18,0.6,0.38,0.75);
  this_leg->SetFillColor(0);
  this_leg->SetBorderSize(0);
  this_leg->SetTextSize(0.035);
  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");
  } else {
    write_warning(__FUNCTION__,"Not implemented yet for mSUGRA");
  }
    
  this_leg->Draw();
  
  string legT5zz="pp #rightarrow  #tilde{g} #tilde{g}, #tilde{g} #rightarrow 2j + #chi^{0}_{2}, #chi^{0}_{2} #rightarrow Z + LSP";
  string legT5zzl2="m(#tilde{q}) >> m(#tilde{g})";
  string legT5zzl3=" GMSB";
  TText *title = write_text(0.18,0.85,legT5zz);
  title->SetTextAlign(11);
  title->SetTextSize(0.035);
  if(scantype!=PlottingSetup::mSUGRA) title->Draw("same");
  TText *title2 = write_text(0.18,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();
  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(2);
  line->SetLineWidth(2);
  line->Draw("same");
}

TH2F* prep_histo(TH2F *oldhist, int scantype) {///DONE
  TH2F *histo = new TH2F(oldhist->GetName(),oldhist->GetName(),
                         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* 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->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",1200,1200);
  ca->Divide(2,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(0.95);
  limits->GetZaxis()->CenterTitle();
  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());
  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
  double refMult = 3.0; 
  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 xsLimitBelow = 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(3);

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

  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->SetLineColor(kBlack);
        desertline->Draw("same");
  }

  fill_with_text(exclline,excllined3,excllinet3,finalcanvas,scantype);
  stringstream real;
  real << "Limits/final_exclusion__" << limits->GetName();
  exclline->Draw("l");
  thinexcline->Draw("");
  excllinet3->Draw("");
  thinexclinet3->Draw("");
  excllined3->Draw("");
  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;
  TH1* H;

  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;
  stringstream flatname;
  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);
    float value=anything.second;
    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<points.size();i++) {
    xpoints[i]=points[i].first;
    spoints[i]=points[i].second;
    if(scantype==PlottingSetup::mSUGRA) {
        if(i>1&&i<points.size()-1) spoints[i]=(1.0/3.0)*(points[i-1].second+points[i].second+points[i+1].second);
        if(i>2&&i<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<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<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;
}

void draw_mSUGRA_exclusion(TH2F *crosssection, TH2F *limitmap, TH2F *expmap, TH2F *expplusmap, TH2F *expminusmap, TH2F *exp2plusmap, TH2F *exp2minusmap) {
  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 *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");
  decorate_mSUGRA(cleanhisto,te,expected,expected2,observed);
  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",1000,1000);
  set_range(crosssection,true,false);
  set_range(limits,true,false);
  set_range(limitmap,true,false);
  
  overview->Divide(2,2);
  overview->cd(1);
  overview->cd(1)->SetLogz(1);
  crosssection->GetZaxis()->SetRangeUser(0.0001,100);
  crosssection->Draw("COLZ");
  TText *title0 = write_title("Cross Section");
  title0->Draw("same");
  overview->cd(2);
  overview->cd(2)->SetLogz(1);
  limits->GetZaxis()->SetRangeUser(0.1,100);
  limits->Draw("COLZ");
  TText *title1 = write_title("Cross Section Upper Limit");
  title1->Draw("same");
  overview->cd(3);
  limitmap->Draw("COLZ");
  TText *title2 = write_title("UL/XS");
  title2->Draw("same");
  observed->Draw("c");
  overview->cd(4);
  decorate_mSUGRA(cleanhisto,overview->cd(4),expected,expected2,observed);
  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;

  
}

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<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, 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<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<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);
  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->Draw("COL");
  gPad->Update();
  if(scantype!=PlottingSetup::mSUGRA) bestlimitsource->Draw("TEXT,same");
  TLegend *sourceleg = new TLegend(0.2,0.6,0.55,0.85);
  for(int i=0;i<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();

  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 bestexplimit;
  delete bestlimitsource;
  return bestlimit;
}


void create_exclusion_plots(vector<TH2F*> limits, int scantype) {
  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<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;
  
  for(int ilim=0;ilim<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("_absolutecrosssectionmap")) crosssections.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]);
  }
    
  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, exp1mlimits, exp1plimits, exp2mlimits, exp2plimits, bestexplimits);
  bestlimits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
  
  for(int ilim=0;ilim<limits.size();ilim++) {
    limits[ilim]->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
  }
  
  if(scantype!=PlottingSetup::mSUGRA) {
    for(int ilim=0;ilim<obslimits.size();ilim++) make_SMS_exclusion(obslimits[ilim],xsec,scantype);
    make_SMS_exclusion(bestlimits,xsec,scantype);
  } else {
    for(int ilim=0;ilim<obslimits.size();ilim++) {
      draw_mSUGRA_exclusion(crosssections[0], obslimits[ilim], explimits[ilim], exp1mlimits[ilim], exp1plimits[ilim], exp2mlimits[ilim], exp2plimits[ilim]);
    }
    draw_mSUGRA_exclusion(crosssections[0], bestlimits, bestexplimits[0], bestexplimits[1], bestexplimits[2], bestexplimits[3], bestexplimits[4]);
  }
  delete bestlimits;
}

void decorate_mSUGRA(TH2F *cleanhisto, TVirtualPad *cvsSys,TGraph *expected,TGraph *expected2,TGraph *observed) {
  cvsSys->SetRightMargin(standardmargin);
  Int_t tanBeta_ = 10;
  Bool_t plotLO_ = false;
  
  //convert tanb value to string
  std::stringstream tmp;
  tmp << tanBeta_;
  TString tanb( tmp.str() );

  //set old exclusion Limits
  TGraph* LEP_ch = set_lep_ch(tanBeta_);
  TGraph* LEP_sl = set_lep_sl(tanBeta_);//slepton curve
  TGraph* TEV_sg_cdf = set_tev_sg_cdf(tanBeta_);//squark gluino cdf
  TGraph* TEV_sg_d0 = set_tev_sg_d0(tanBeta_);//squark gluino d0
  TGraph* stau   = set_tev_stau(tanBeta_);//stau 
  TGraph* NoEWSB = set_NoEWSB(tanBeta_); 

  TGraph* TEV_sn_d0_1 = set_sneutrino_d0_1(tanBeta_);
  TGraph* TEV_sn_d0_2 = set_sneutrino_d0_2(tanBeta_);

  //constant ssqquark and gluino lines
  TF1* lnsq[10];
  TF1* lngl[10];
  
  TLatex* sq_text[10];
  TLatex* gl_text[10];
  
  for(int i = 0; i < 6; i++){
    lnsq[i] = constant_squark(tanBeta_,i);
    sq_text[i] = constant_squark_text(i,*lnsq[i],tanBeta_);
    lngl[i] = constant_gluino(tanBeta_,i);
    gl_text[i] = constant_gluino_text(i,*lngl[i]);
  }

  //Legends
  TLegend* legst  = makeStauLegend(0.05,tanBeta_);
  TLegend* legNoEWSB  = makeNoEWSBLegend(0.05,tanBeta_);
  TLegend* legexp = makeExpLegend( *TEV_sg_cdf,*TEV_sg_d0,*LEP_ch,*LEP_sl,*TEV_sn_d0_1,0.035,tanBeta_);
 
  TEV_sn_d0_1->SetLineWidth(1);
  TEV_sn_d0_2->SetLineWidth(1);
  TEV_sg_d0->SetLineWidth(1);
  
  double m0min = 0;
  if (tanBeta_ == 50) m0min=200;
  TH2F* hist = new TH2F("h","h",100,m0min,1000,100,120,700);
  hist->Draw();  
  hist->GetXaxis()->SetTitle("m_{0} [GeV]");
  hist->GetXaxis()->CenterTitle();
  hist->GetYaxis()->SetTitle("m_{1/2} [GeV]");
  hist->GetYaxis()->CenterTitle();
  hist->GetXaxis()->SetTitleSize(0.04);
  hist->GetYaxis()->SetTitleSize(0.04);
  hist->GetXaxis()->SetTitleOffset(1.2);
  hist->GetYaxis()->SetTitleOffset(1.5);
  hist->GetXaxis()->SetNdivisions(506);
  hist->GetYaxis()->SetNdivisions(506);

  int col[]={2,3,4};

  TLegend* myleg;

  if( plotLO_ ) myleg = new TLegend(0.3,0.65,0.65,0.8,NULL,"brNDC");
  else          myleg = new TLegend(0.25,0.76,0.44,0.91,NULL,"brNDC");

  myleg->SetFillColor(0); 
  myleg->SetShadowColor(0);
  myleg->SetTextSize(0.04);
  myleg->SetBorderSize(0);

  //constant squark and gluino mass contours
  for (int it=0;it<5;it++) {   
    lngl[it]->Draw("same");   
    lnsq[it]->Draw("same");
    sq_text[it]->Draw();
    gl_text[it]->Draw();
  }
  
//exclusion limits previous experiments
  if(tanBeta_ == 3){
    TEV_sn_d0_1->Draw("fsame");
    TEV_sn_d0_2->Draw("fsame");
  }
  LEP_ch->Draw("fsame");
  if (tanBeta_ != 50) LEP_sl->Draw("fsame");

  //remove CDF/D0 excluded regions
  TEV_sg_cdf->Draw("fsame");
  TEV_sg_d0->Draw("same");  
  TEV_sg_d0->Draw("fsame");

  //other labels
  Double_t xpos = 0;
  Double_t xposi = 0;
  Double_t ypos = 0;
  if(tanBeta_ == 50) xposi = 100;
  if(tanBeta_ == 50) xpos = 200;
  if(tanBeta_ == 50) ypos = -10;

  TString text_tanBeta;
  text_tanBeta =  "tan#beta = "+tanb+",  A_{0} = 0,  #mu > 0";
  TLatex* cmssmpars = new TLatex(/*530.+xpos,690.+ypos-130*/150,650,text_tanBeta);

  cmssmpars->SetTextSize(0.03);
  cmssmpars->Draw("same");

  TLatex* lep_chargino = new TLatex(250,135,"LEP2 #tilde{#chi}_{1}^{#pm}");
  lep_chargino->SetTextSize(0.03);
  lep_chargino->SetTextFont(42);

  TLatex* lep_slepton = new TLatex(26,190,"LEP2 #tilde{#font[12]{l}}^{#pm}");
  lep_slepton->SetTextSize(0.03);
  lep_slepton->SetTextAngle(-83);
  lep_slepton->SetTextFont(42);

  if(draw2sigma) expected2->Draw("f");
  expected->SetLineColor(expected->GetFillColor());
  expected2->SetLineColor(expected2->GetFillColor());
  expected->Draw("f");
  observed->Draw("c");
  
  stau->Draw("fsame");
  NoEWSB->Draw("fsame");
 
  legexp->AddEntry(observed,"Observed limit","l");
  legexp->AddEntry(expected,"Expected 1#sigma limit","f");
  if(draw2sigma) legexp->AddEntry(expected2,"Expected 2#sigma limit","f");
  legexp->SetY1(0.60);
  legexp->SetX1(0.55);
  legexp->Draw();
  legst->Draw();

  hist->Draw("sameaxis");
  cvsSys->RedrawAxis();
  cvsSys->Update();
  DrawPrelim();
}

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;
    bool turning=false;
    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) {
  if(scantype==PlottingSetup::mSUGRA) {
    histo->GetXaxis()->SetRangeUser(0,2000);
    histo->GetYaxis()->SetRangeUser(0,800);
    histo->GetXaxis()->SetTitle("m_{0} [GeV]");
    histo->GetYaxis()->SetTitle("m_{1/2} [GeV]");
    histo->GetXaxis()->SetRangeUser(0,2000);
    histo->GetXaxis()->SetNdivisions(504,true);
  }
  if(scantype==PlottingSetup::SMS) {
    histo->GetXaxis()->SetRangeUser(50,1200);
    histo->GetYaxis()->SetRangeUser(50,1200);
    histo->GetXaxis()->SetTitle("m_{#tilde{g}} [GeV]");
    histo->GetYaxis()->SetTitle("m_{LSP} [GeV]");
    histo->GetXaxis()->SetTitleSize(0.04);
    histo->GetYaxis()->SetTitleSize(0.04);
    histo->GetZaxis()->SetTitleSize(0.04);
    histo->GetXaxis()->SetTitleOffset(1.2);
    histo->GetYaxis()->SetTitleOffset(1.5);
    if(pushoutyz) histo->GetZaxis()->SetTitleOffset(1.6);
  }
  if(scantype==PlottingSetup::GMSB) {
    histo->GetXaxis()->SetRangeUser(100,1200);
    histo->GetYaxis()->SetRangeUser(100,1200);
    histo->GetXaxis()->SetTitle("m_{#tilde{g}} [GeV]");
    histo->GetYaxis()->SetTitle("m_{#chi} [GeV]");
    histo->GetXaxis()->SetTitleSize(0.04);
    histo->GetYaxis()->SetTitleSize(0.04);
    histo->GetZaxis()->SetTitleSize(0.04);
    histo->GetXaxis()->SetTitleOffset(1.2);
    histo->GetYaxis()->SetTitleOffset(1.5);
    if(pushoutyz) histo->GetZaxis()->SetTitleOffset(1.6);
  }

  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) {
  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);
  
  
  if(name.Contains("efficiency")) {
    histo->GetZaxis()->SetTitle("A #times #varepsilon (#geq 1 Z(ll))");
    //histo->GetZaxis()->SetRangeUser(0.0,0.3);
  }
  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->GetZaxis()->CenterTitle();
  gStyle->SetStripDecimals(false);
  histo->GetXaxis()->SetDecimals(true);
  
  histo->Draw("COLZ");
  DrawPrelim();
  
  CompleteSave(can,(saveto+(string)histo->GetName()));

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

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

void process_file(TFile* file, float stdmargin) {
  standardmargin=stdmargin;
  xsecfilename="reference_xSec_SMS.root";

  // 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;
  
  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("absolutecrosssectionmap")) 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;
    }
  if(systematics_histos.size()>0) make_all_syst_plots(systematics_histos,scantype);
  if(limits_histos.size()>0) create_exclusion_plots(limits_histos,scantype);
}

Revision:	1.17
Committed:	Tue Dec 13 10:44:44 2011 UTC (13 years, 4 months ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.16:	+11 -2 lines
Log Message:	Fix special case: if the first histo of expected limits has an empty bin (i.e. bin content 0), then there's no chance any of the following histos will ever beat that and fill in something meaningful ...