| 1 |
#include <iostream>
|
| 2 |
#include <vector>
|
| 3 |
#include <sys/stat.h>
|
| 4 |
#include <sstream>
|
| 5 |
|
| 6 |
#include <TCut.h>
|
| 7 |
#include <TLatex.h>
|
| 8 |
#include <TROOT.h>
|
| 9 |
#include <TCanvas.h>
|
| 10 |
#include <TPad.h>
|
| 11 |
#include <TFile.h>
|
| 12 |
#include <TTree.h>
|
| 13 |
#include <TString.h>
|
| 14 |
#include <TMath.h>
|
| 15 |
#include <stdio.h>
|
| 16 |
#include <stdlib.h>
|
| 17 |
#include <TH2.h>
|
| 18 |
#include <TColor.h>
|
| 19 |
#include <TStyle.h>
|
| 20 |
#include <TText.h>
|
| 21 |
#include <TKey.h>
|
| 22 |
#include <TPaletteAxis.h>
|
| 23 |
#include <TGraph.h>
|
| 24 |
#include <TLegend.h>
|
| 25 |
#include <TLegendEntry.h>
|
| 26 |
#include "external/ExclusionPlot/Functions.C"
|
| 27 |
|
| 28 |
using namespace std;
|
| 29 |
|
| 30 |
bool draweachone=false;
|
| 31 |
bool draw2sigma=true;
|
| 32 |
|
| 33 |
float limits_lower_bound=0.1;
|
| 34 |
float limits_upper_bound=10;
|
| 35 |
float limits_ratio_lower_bound=0.005;
|
| 36 |
float limits_ratio_upper_bound=10;
|
| 37 |
|
| 38 |
float standardmargin=0.1;
|
| 39 |
float indentedmargin=0.16;
|
| 40 |
|
| 41 |
bool drawefficiencydesertline=false;
|
| 42 |
|
| 43 |
string xsecfilename;
|
| 44 |
|
| 45 |
void set_range(TH2F *histo, bool ismSUGRA, bool pushoutlabels);
|
| 46 |
void smooth_line(TGraph *gr);
|
| 47 |
void decorate_mSUGRA(TH2F *cleanhisto, TVirtualPad *cvsSys,TGraph *expected,TGraph *expected2,TGraph *observed);
|
| 48 |
TH2F* adjust_histo(TH2F *oldhist, TH2F *refhisto);
|
| 49 |
|
| 50 |
void fill_with_text(TH2F *real, TH2F *down, TH2F *up, TVirtualPad *can, bool ismSUGRA) {
|
| 51 |
can->cd();
|
| 52 |
TLegend* this_leg = new TLegend(0.18,0.6,0.38,0.75);
|
| 53 |
this_leg->SetFillColor(0);
|
| 54 |
this_leg->SetBorderSize(0);
|
| 55 |
this_leg->SetTextSize(0.035);
|
| 56 |
if(!ismSUGRA) {
|
| 57 |
this_leg->AddEntry(real,"#sigma^{prod} = #sigma^{NLO-QCD}" , "l");
|
| 58 |
this_leg->AddEntry(up,"#sigma^{prod} = 3 #times #sigma^{NLO-QCD}" , "l");
|
| 59 |
this_leg->AddEntry(down,"#sigma^{prod} = 1/3 #times #sigma^{NLO-QCD}" , "l");
|
| 60 |
} else {
|
| 61 |
write_warning(__FUNCTION__,"Not implemented yet for mSUGRA");
|
| 62 |
}
|
| 63 |
|
| 64 |
this_leg->Draw();
|
| 65 |
|
| 66 |
string legT5zz="pp #rightarrow #tilde{g} #tilde{g}, #tilde{g} #rightarrow 2j + #chi^{0}, #chi^{0} #rightarrow Z + LSP";
|
| 67 |
string legT5zzl2="m(#tilde{q}) >> m(#tilde{g})";
|
| 68 |
TText *title = write_text(0.18,0.85,legT5zz);
|
| 69 |
title->SetTextAlign(11);
|
| 70 |
title->SetTextSize(0.035);
|
| 71 |
if(!ismSUGRA) title->Draw("same");
|
| 72 |
TText *title2 = write_text(0.18,0.79,legT5zzl2);
|
| 73 |
title2->SetTextAlign(11);
|
| 74 |
title2->SetTextSize(0.035);
|
| 75 |
if(!ismSUGRA) title2->Draw("same");
|
| 76 |
DrawPrelim();
|
| 77 |
TLine *line;
|
| 78 |
line = new TLine(50.+75., 50., 1200., 1200-75.);
|
| 79 |
line->SetLineStyle(2);
|
| 80 |
line->SetLineWidth(2);
|
| 81 |
line->Draw("same");
|
| 82 |
}
|
| 83 |
|
| 84 |
TH2F* prep_histo(TH2F *oldhist, bool ismSUGRA) {///DONE
|
| 85 |
TH2F *histo = new TH2F(oldhist->GetName(),oldhist->GetName(),
|
| 86 |
oldhist->GetNbinsX(),oldhist->GetXaxis()->GetBinLowEdge(1),oldhist->GetXaxis()->GetBinLowEdge(oldhist->GetNbinsX())+oldhist->GetXaxis()->GetBinWidth(oldhist->GetNbinsX()),
|
| 87 |
oldhist->GetNbinsY(),oldhist->GetYaxis()->GetBinLowEdge(1),oldhist->GetYaxis()->GetBinLowEdge(oldhist->GetNbinsY())+oldhist->GetYaxis()->GetBinWidth(oldhist->GetNbinsY()));
|
| 88 |
|
| 89 |
for(int ix=1;ix<=oldhist->GetNbinsX();ix++) {
|
| 90 |
for(int iy=1;iy<=oldhist->GetNbinsX();iy++) {
|
| 91 |
histo->SetBinContent(ix,iy,oldhist->GetBinContent(ix,iy));
|
| 92 |
}
|
| 93 |
}
|
| 94 |
return histo;
|
| 95 |
}
|
| 96 |
|
| 97 |
TH2F* make_exclusion_shape(TH2F *excl, int isprimary) {
|
| 98 |
TH2F *exclusion = (TH2F*)excl->Clone("exclusion");
|
| 99 |
for(int i=1; i<(excl->GetNbinsX()+1); i++) {
|
| 100 |
for(int j=1; j<(excl->GetNbinsY()+1); j++) {
|
| 101 |
if(excl->GetBinContent(i,j)<1&&excl->GetBinContent(i,j)>0) exclusion->SetBinContent(i,j,0.01);
|
| 102 |
else exclusion->SetBinContent(i,j,0);
|
| 103 |
}
|
| 104 |
}
|
| 105 |
exclusion->SetLineColor(kBlue);
|
| 106 |
exclusion->SetLineWidth(2);
|
| 107 |
exclusion->SetLineStyle(isprimary);
|
| 108 |
return exclusion;
|
| 109 |
}
|
| 110 |
|
| 111 |
void produce_extensive_plots(TH2F *xsec,TH2F *limits,TH2F *exclusionshape,TH2F *exclusionshapet3,TH2F *exclusionshaped3, bool ismSUGRA) {
|
| 112 |
TCanvas *ca = new TCanvas("ca","ca",1200,1200);
|
| 113 |
ca->Divide(2,2);
|
| 114 |
ca->cd(1);
|
| 115 |
ca->cd(1)->SetLogz(1);
|
| 116 |
xsec->GetZaxis()->SetRangeUser(0.001,1000);
|
| 117 |
xsec->Draw("COLZ");
|
| 118 |
TText *title0 = write_title("Reference Cross Section");
|
| 119 |
title0->Draw("same");
|
| 120 |
ca->cd(2);
|
| 121 |
ca->cd(2)->SetLogz(1);
|
| 122 |
limits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
|
| 123 |
limits->Draw("COLZ");
|
| 124 |
TText *title = write_title("Cross Section Upper Limit");
|
| 125 |
title->Draw("same");
|
| 126 |
ca->cd(3);
|
| 127 |
ca->cd(3)->SetLogz(1);
|
| 128 |
TH2F *limit_ref = (TH2F*)limits->Clone("limit_ref");
|
| 129 |
limit_ref->Divide(xsec);
|
| 130 |
limit_ref->GetZaxis()->SetRangeUser(limits_ratio_lower_bound,limits_ratio_upper_bound);
|
| 131 |
limit_ref->Draw("COLZ");
|
| 132 |
TText *title2 = write_title("Cross Section UL / XS");
|
| 133 |
title2->Draw("same");
|
| 134 |
ca->cd(4);
|
| 135 |
ca->cd(4)->SetLogz(1);
|
| 136 |
limits->SetTitle("");
|
| 137 |
limits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
|
| 138 |
limits->SetZTitle("95% CL upper limit on #sigma [pb]");
|
| 139 |
|
| 140 |
limits->GetZaxis()->SetTitleOffset(0.95);
|
| 141 |
limits->GetZaxis()->CenterTitle();
|
| 142 |
exclusionshapet3->SetLineStyle(2);
|
| 143 |
exclusionshaped3->SetLineStyle(3);
|
| 144 |
exclusionshape->GetZaxis()->SetRangeUser(0,500*0.1);
|
| 145 |
exclusionshapet3->GetZaxis()->SetRangeUser(0,500*0.1);
|
| 146 |
exclusionshaped3->GetZaxis()->SetRangeUser(0,500*0.1);
|
| 147 |
limits->Draw("COLZ");
|
| 148 |
exclusionshape->Draw("CONT1,same");
|
| 149 |
exclusionshapet3->Draw("CONT1,same");
|
| 150 |
exclusionshaped3->Draw("CONT1,same");
|
| 151 |
stringstream partial;
|
| 152 |
partial << "Limits/exclusion__" << limits->GetName();
|
| 153 |
fill_with_text(exclusionshape,exclusionshaped3,exclusionshapet3,ca->cd(4),ismSUGRA);
|
| 154 |
CompleteSave(ca,partial.str());
|
| 155 |
delete ca;
|
| 156 |
}
|
| 157 |
|
| 158 |
|
| 159 |
|
| 160 |
void make_SMS_exclusion(TH2F *rawlimits,TH2F *xsec,bool ismSUGRA) {
|
| 161 |
TH2F *limits = prep_histo(rawlimits,ismSUGRA); // this is to be independent of the style used at creation time
|
| 162 |
//here we get some limits and the cross section; we want to make an exclusion plot!
|
| 163 |
TH2F *rellimits = (TH2F*)limits->Clone("rellimits");
|
| 164 |
TH2F *rellimitsd3 = (TH2F*)limits->Clone("rellimitsd3"); // one third the cross section ("divided by 3" -> d3)
|
| 165 |
TH2F *rellimitst3 = (TH2F*)limits->Clone("rellimitst3"); // three times the cross section ("times 3" -> t3)
|
| 166 |
|
| 167 |
if(!xsec ) {
|
| 168 |
cout << "Watch out, cross section map is invalid!" << endl;
|
| 169 |
return;
|
| 170 |
}
|
| 171 |
|
| 172 |
rellimits->Divide(xsec);
|
| 173 |
|
| 174 |
|
| 175 |
for(int i=1;i<=rellimits->GetNbinsX();i++) {
|
| 176 |
for(int j=1;j<=rellimits->GetNbinsY();j++) {
|
| 177 |
rellimitst3->SetBinContent(i,j,(rellimits->GetBinContent(i,j))/3.0);
|
| 178 |
rellimitsd3->SetBinContent(i,j,(rellimits->GetBinContent(i,j))*3.0);
|
| 179 |
}
|
| 180 |
}
|
| 181 |
|
| 182 |
TH2F *exclusionshape = make_exclusion_shape(rellimits,1);
|
| 183 |
TH2F *exclusionshapet3 = make_exclusion_shape(rellimitst3,2);
|
| 184 |
TH2F *exclusionshaped3 = make_exclusion_shape(rellimitsd3,3);
|
| 185 |
|
| 186 |
//Now let's produce the plots!
|
| 187 |
|
| 188 |
set_range(xsec,ismSUGRA,false);
|
| 189 |
set_range(limits,ismSUGRA,false);
|
| 190 |
set_range(exclusionshape,ismSUGRA,false);
|
| 191 |
set_range(exclusionshaped3,ismSUGRA,false);
|
| 192 |
set_range(exclusionshapet3,ismSUGRA,false);
|
| 193 |
|
| 194 |
produce_extensive_plots(xsec,limits,exclusionshape,exclusionshapet3,exclusionshaped3,ismSUGRA);
|
| 195 |
|
| 196 |
TCanvas *finalcanvas = new TCanvas("finalcanvas","finalcanvas");
|
| 197 |
finalcanvas->SetLogz(1);
|
| 198 |
finalcanvas->cd();
|
| 199 |
limits->Draw("COLZ");
|
| 200 |
exclusionshape->Draw("CONT1,same");
|
| 201 |
exclusionshapet3->Draw("CONT1,same");
|
| 202 |
exclusionshaped3->Draw("CONT1,same");
|
| 203 |
TLine *desertline;
|
| 204 |
if(drawefficiencydesertline) {
|
| 205 |
desertline = new TLine(375,50,1200,875);
|
| 206 |
desertline->SetLineWidth(3);
|
| 207 |
desertline->SetLineColor(kBlack);
|
| 208 |
desertline->Draw("same");
|
| 209 |
}
|
| 210 |
|
| 211 |
fill_with_text(exclusionshape,exclusionshaped3,exclusionshapet3,finalcanvas,ismSUGRA);
|
| 212 |
stringstream real;
|
| 213 |
real << "Limits/final_exclusion__" << limits->GetName();
|
| 214 |
CompleteSave(finalcanvas,real.str());
|
| 215 |
|
| 216 |
delete finalcanvas;
|
| 217 |
}
|
| 218 |
|
| 219 |
TH1* getHisto(char * filename, char* histoName, char * dirName, int nBin, double lumi)
|
| 220 |
{
|
| 221 |
TH1 * hpt_=0;
|
| 222 |
TFile *file0 = TFile::Open(filename);
|
| 223 |
if(!file0) return hpt_;
|
| 224 |
TDirectory *dir;
|
| 225 |
TH2D * hMuPt;
|
| 226 |
TH1* H;
|
| 227 |
|
| 228 |
if(dirName == "0") {
|
| 229 |
hMuPt = (TH2D*) file0->Get(histoName);
|
| 230 |
} else {
|
| 231 |
dir = (TDirectory*) file0->Get(dirName);
|
| 232 |
if(!dir) return hpt_;
|
| 233 |
hMuPt = (TH2D*) dir->Get(histoName);
|
| 234 |
}
|
| 235 |
|
| 236 |
if(hMuPt) {
|
| 237 |
hpt_ = (TH1*) hMuPt->Clone();
|
| 238 |
hpt_->Sumw2();
|
| 239 |
hpt_->Scale(1./lumi); // this take into into account the luminosity
|
| 240 |
hpt_->SetLineWidth(2);
|
| 241 |
hpt_->Rebin(nBin);
|
| 242 |
double nBinX=hpt_->GetNbinsX();
|
| 243 |
|
| 244 |
// overFlow
|
| 245 |
hpt_->SetBinContent((int)nBinX,hpt_->GetBinContent((int)nBinX)+hpt_->GetBinContent((int)nBinX+1));
|
| 246 |
hpt_->SetDirectory(0);
|
| 247 |
file0->Close();
|
| 248 |
hpt_->SetLineWidth(3);
|
| 249 |
}
|
| 250 |
return hpt_;
|
| 251 |
}
|
| 252 |
|
| 253 |
pair <float,float> find_point(TH2F* histo, int xbin, bool mSUGRA=true) {
|
| 254 |
TH1F *flathisto;
|
| 255 |
if(mSUGRA) flathisto = new TH1F("flat","flat",histo->GetNbinsY(),histo->GetYaxis()->GetBinLowEdge(1),histo->GetYaxis()->GetBinLowEdge(histo->GetNbinsY())+histo->GetYaxis()->GetBinWidth(histo->GetNbinsY()));
|
| 256 |
else flathisto = new TH1F("flat","flat",histo->GetNbinsX(),histo->GetXaxis()->GetBinLowEdge(1),histo->GetXaxis()->GetBinLowEdge(histo->GetNbinsX())+histo->GetXaxis()->GetBinWidth(histo->GetNbinsX()));
|
| 257 |
|
| 258 |
int nbins=histo->GetNbinsY();
|
| 259 |
if(!mSUGRA) nbins=histo->GetNbinsX();
|
| 260 |
for(int i=1;i<nbins;i++) {
|
| 261 |
float value=(histo->GetBinContent(xbin,i));
|
| 262 |
if(value<20&&value>0) flathisto->SetBinContent(i,value);
|
| 263 |
}
|
| 264 |
|
| 265 |
float pointone=-100;
|
| 266 |
nbins=flathisto->GetNbinsX();
|
| 267 |
if(!mSUGRA) nbins=flathisto->GetNbinsY();
|
| 268 |
for(int i=nbins;i>=1;i--) {
|
| 269 |
if(pointone<0&&flathisto->GetBinContent(i)<1&&flathisto->GetBinContent(i)>0) pointone=flathisto->GetBinLowEdge(i)+flathisto->GetBinWidth(i);
|
| 270 |
}
|
| 271 |
pair <float,float> anything;
|
| 272 |
if(mSUGRA) anything.first=histo->GetXaxis()->GetBinCenter(xbin);
|
| 273 |
else anything.first=histo->GetYaxis()->GetBinCenter(xbin);
|
| 274 |
anything.second=pointone;
|
| 275 |
stringstream flatname;
|
| 276 |
delete flathisto;
|
| 277 |
return anything;
|
| 278 |
}
|
| 279 |
|
| 280 |
pair <float,float> find_interpolated_point(TH2F* histo, int xbin, bool mSUGRA=true) {
|
| 281 |
|
| 282 |
int minaccept=4;
|
| 283 |
TCanvas *flatcan = new TCanvas("flatcan","flatcan");
|
| 284 |
stringstream histoname;
|
| 285 |
histoname << "exclusion shape for x bin " << xbin;
|
| 286 |
TH1F *flathisto;
|
| 287 |
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()));
|
| 288 |
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()));
|
| 289 |
|
| 290 |
int acceptedpoints=0;
|
| 291 |
int nbins=histo->GetNbinsY();
|
| 292 |
if(!mSUGRA) histo->GetNbinsX();
|
| 293 |
for(int i=1;i<nbins;i++) {
|
| 294 |
float value=0;
|
| 295 |
if(i<=nbins-2) value=((1/3.0)*(histo->GetBinContent(xbin,i)+histo->GetBinContent(xbin+1,i)+histo->GetBinContent(xbin+2,i)));
|
| 296 |
if(i==nbins-1) value=((1/2.0)*(histo->GetBinContent(xbin,i)+histo->GetBinContent(xbin,i+1)));
|
| 297 |
if(i==nbins) value=(histo->GetBinContent(xbin,i));
|
| 298 |
if(value<20&&value>0) {
|
| 299 |
flathisto->SetBinContent(i,value);
|
| 300 |
flathisto->SetBinError(i,TMath::Sqrt(value));
|
| 301 |
acceptedpoints++;
|
| 302 |
}
|
| 303 |
}
|
| 304 |
|
| 305 |
float pointone=-100;
|
| 306 |
TLine *excluded;
|
| 307 |
if(acceptedpoints>minaccept) {
|
| 308 |
flathisto->Fit("expo","Q");
|
| 309 |
TF1 *fitfunc = (TF1*)flathisto->GetFunction("expo");
|
| 310 |
pointone=-(fitfunc->GetParameter(0)/fitfunc->GetParameter(1));
|
| 311 |
excluded=new TLine(pointone,0,pointone,10);
|
| 312 |
}
|
| 313 |
|
| 314 |
pair <float,float> anything;
|
| 315 |
anything.first=histo->GetXaxis()->GetBinCenter(xbin);
|
| 316 |
anything.second=pointone;
|
| 317 |
stringstream flatname;
|
| 318 |
flathisto->GetYaxis()->SetRangeUser(0,10);
|
| 319 |
flathisto->Draw();
|
| 320 |
if(acceptedpoints>minaccept) excluded->SetLineColor(kGreen);
|
| 321 |
if(acceptedpoints>minaccept) excluded->SetLineStyle(2);
|
| 322 |
if(acceptedpoints>minaccept) excluded->SetLineWidth(2);
|
| 323 |
if(acceptedpoints>minaccept) excluded->Draw("same");
|
| 324 |
flatname << "Limits/partials/partial_" << xbin << "___" << histo->GetName() << ".png";
|
| 325 |
if(draweachone) CompleteSave(flatcan,flatname.str());
|
| 326 |
delete flatcan;
|
| 327 |
delete flathisto;
|
| 328 |
return anything;
|
| 329 |
}
|
| 330 |
|
| 331 |
TGraph* get_exclusion_line(TH2F *exclusionhisto) {
|
| 332 |
exclusionhisto->SetStats(0);
|
| 333 |
exclusionhisto->GetZaxis()->SetRangeUser(0,2);
|
| 334 |
|
| 335 |
vector<pair <float,float> > points;
|
| 336 |
|
| 337 |
for(int i=1;i<exclusionhisto->GetNbinsX();i++) {
|
| 338 |
pair<float,float> anything = find_point(exclusionhisto,i);
|
| 339 |
pair<float,float> intthing = find_interpolated_point(exclusionhisto,i);
|
| 340 |
float value=anything.second;
|
| 341 |
if(intthing.second>anything.second) anything.second=intthing.second;
|
| 342 |
if(anything.second>100&&anything.second<1000) points.push_back(anything);
|
| 343 |
}
|
| 344 |
|
| 345 |
Double_t xpoints[points.size()];
|
| 346 |
Double_t spoints[points.size()];
|
| 347 |
|
| 348 |
TGraph *graph = new TGraph(points.size());
|
| 349 |
graph->SetLineColor(kBlack);
|
| 350 |
graph->SetLineWidth(2);
|
| 351 |
|
| 352 |
for(int i=0;i<points.size();i++) {
|
| 353 |
xpoints[i]=points[i].first;
|
| 354 |
if(i>1&&i<points.size()-1) spoints[i]=(1.0/3.0)*(points[i-1].second+points[i].second+points[i+1].second);
|
| 355 |
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);
|
| 356 |
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);
|
| 357 |
else spoints[i]=points[i].second;
|
| 358 |
graph->SetPoint(i,points[i].first,spoints[i]);
|
| 359 |
}
|
| 360 |
return graph;
|
| 361 |
}
|
| 362 |
|
| 363 |
void draw_mSUGRA_exclusion(TH2F *crosssection, TH2F *limitmap, TH2F *expmap, TH2F *expplusmap, TH2F *expminusmap, TH2F *exp2plusmap, TH2F *exp2minusmap) {
|
| 364 |
TH2F *cleanhisto = (TH2F*)limitmap->Clone("clean");
|
| 365 |
for(int ix=1;ix<=cleanhisto->GetNbinsX();ix++) {
|
| 366 |
for(int iy=1;iy<=cleanhisto->GetNbinsY();iy++) {
|
| 367 |
cleanhisto->SetBinContent(ix,iy,0);
|
| 368 |
}
|
| 369 |
}
|
| 370 |
|
| 371 |
TH2F *limits = (TH2F*)limitmap->Clone("limits");
|
| 372 |
set_range(limits,true,false);
|
| 373 |
limitmap->Divide(crosssection);
|
| 374 |
expminusmap->Divide(crosssection);
|
| 375 |
expplusmap->Divide(crosssection);
|
| 376 |
exp2minusmap->Divide(crosssection);
|
| 377 |
exp2plusmap->Divide(crosssection);
|
| 378 |
expmap->Divide(crosssection);
|
| 379 |
TGraph *observed = get_exclusion_line(limitmap);
|
| 380 |
observed->SetLineColor(kRed);
|
| 381 |
TGraph *expminus = get_exclusion_line(expminusmap);
|
| 382 |
TGraph *expplus = get_exclusion_line(expplusmap);
|
| 383 |
TGraph *exp2minus;
|
| 384 |
if(draw2sigma) exp2minus = get_exclusion_line(exp2minusmap);
|
| 385 |
TGraph *exp2plus;
|
| 386 |
if(draw2sigma) exp2plus = get_exclusion_line(exp2plusmap);
|
| 387 |
|
| 388 |
TGraph *expected = new TGraph(expminus->GetN()+expplus->GetN());
|
| 389 |
TGraph *expected2;
|
| 390 |
if(draw2sigma) expected2 = new TGraph(exp2minus->GetN()+exp2plus->GetN());
|
| 391 |
for(int i=0;i<=expminus->GetN();i++) {
|
| 392 |
Double_t x,y;
|
| 393 |
expminus->GetPoint(i,x,y);
|
| 394 |
expected->SetPoint(i,x,y);
|
| 395 |
}
|
| 396 |
|
| 397 |
for(int i=0;i<=exp2minus->GetN();i++) {
|
| 398 |
Double_t x,y;
|
| 399 |
exp2minus->GetPoint(i,x,y);
|
| 400 |
expected2->SetPoint(i,x,y);
|
| 401 |
}
|
| 402 |
for(int i=exp2plus->GetN()-1;i>=0;i--) {
|
| 403 |
Double_t x,y;
|
| 404 |
exp2plus->GetPoint(i,x,y);
|
| 405 |
expected2->SetPoint(exp2minus->GetN()+(exp2plus->GetN()-i),x,y);
|
| 406 |
}
|
| 407 |
for(int i=expplus->GetN()-1;i>=0;i--) {
|
| 408 |
Double_t x,y;
|
| 409 |
expplus->GetPoint(i,x,y);
|
| 410 |
expected->SetPoint(expminus->GetN()+(expplus->GetN()-i),x,y);
|
| 411 |
}
|
| 412 |
expected->SetFillColor(TColor::GetColor("#9FF781"));
|
| 413 |
if(draw2sigma) expected2->SetFillColor(TColor::GetColor("#F3F781"));
|
| 414 |
|
| 415 |
smooth_line(observed);
|
| 416 |
smooth_line(expected);
|
| 417 |
if(draw2sigma) smooth_line(expected2);
|
| 418 |
|
| 419 |
TCanvas *te = new TCanvas("te","te");
|
| 420 |
decorate_mSUGRA(cleanhisto,te,expected,expected2,observed);
|
| 421 |
stringstream saveas;
|
| 422 |
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);
|
| 423 |
else saveas << "Limits/final_exclusion_for_bestlimits";
|
| 424 |
CompleteSave(te,saveas.str());
|
| 425 |
delete te;
|
| 426 |
|
| 427 |
TCanvas *overview = new TCanvas("overview","overview",1000,1000);
|
| 428 |
set_range(crosssection,true,false);
|
| 429 |
set_range(limits,true,false);
|
| 430 |
set_range(limitmap,true,false);
|
| 431 |
|
| 432 |
overview->Divide(2,2);
|
| 433 |
overview->cd(1);
|
| 434 |
overview->cd(1)->SetLogz(1);
|
| 435 |
crosssection->GetZaxis()->SetRangeUser(0.0001,100);
|
| 436 |
crosssection->Draw("COLZ");
|
| 437 |
TText *title0 = write_title("Cross Section");
|
| 438 |
title0->Draw("same");
|
| 439 |
overview->cd(2);
|
| 440 |
overview->cd(2)->SetLogz(1);
|
| 441 |
limits->GetZaxis()->SetRangeUser(0.1,100);
|
| 442 |
limits->Draw("COLZ");
|
| 443 |
TText *title1 = write_title("Cross Section Upper Limit");
|
| 444 |
title1->Draw("same");
|
| 445 |
overview->cd(3);
|
| 446 |
limitmap->Draw("COLZ");
|
| 447 |
TText *title2 = write_title("UL/XS");
|
| 448 |
title2->Draw("same");
|
| 449 |
observed->Draw("c");
|
| 450 |
overview->cd(4);
|
| 451 |
decorate_mSUGRA(cleanhisto,overview->cd(4),expected,expected2,observed);
|
| 452 |
stringstream saveas2;
|
| 453 |
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);
|
| 454 |
else saveas2 << "Limits/exclusion_overview_for_bestlimits";
|
| 455 |
CompleteSave(overview,saveas2.str());
|
| 456 |
delete overview;
|
| 457 |
|
| 458 |
|
| 459 |
}
|
| 460 |
|
| 461 |
TH2F* get_XS(string filename, string mass, TH2F *histo) {///DONE
|
| 462 |
TH1 *hRef=getHisto((char*) filename.c_str(), (char*) mass.c_str(), (char*) "0", 1,1);
|
| 463 |
TH2F * xsec= new TH2F("xsec","xsec",60, 0., 1500., 60, 0., 1500.);
|
| 464 |
for(int i=1; i<(histo->GetNbinsX()+1); i++) {
|
| 465 |
for(int j=1; j<(histo->GetNbinsY()+1); j++) {
|
| 466 |
if(!(histo->GetBinContent(i,j)==0)) xsec->SetBinContent(i,j,hRef->GetBinContent(hRef->FindBin(xsec->GetXaxis()->GetBinCenter(i))));
|
| 467 |
}
|
| 468 |
}
|
| 469 |
return xsec;
|
| 470 |
}
|
| 471 |
|
| 472 |
string give_nice_source_label(string name) {
|
| 473 |
int mappoint=name.find("map");
|
| 474 |
if(mappoint<0||mappoint>500) return name; // this mean that something weird is happening
|
| 475 |
stringstream nice_label;
|
| 476 |
nice_label << "JZB > " << name.substr(mappoint+3,name.size()) << " GeV";
|
| 477 |
return nice_label.str();
|
| 478 |
}
|
| 479 |
|
| 480 |
TH2F *make_best_limits(vector<TH2F*> explimits,vector<TH2F*> obslimits, bool ismSUGRA, vector<TH2F*> exp1mlimits, vector<TH2F*> exp1plimits, vector<TH2F*> exp2mlimits, vector<TH2F*> exp2plimits, vector<TH2F*> &allbestexplimits) {
|
| 481 |
if(obslimits.size()==0) {
|
| 482 |
write_warning(__FUNCTION__,"There are no observed limits! Cannot continue!");
|
| 483 |
TH2F *err;
|
| 484 |
return err;
|
| 485 |
}
|
| 486 |
if(explimits.size()==0) {
|
| 487 |
write_warning(__FUNCTION__,"There are no expected limits! Will compute best limits based on observed limits! (WATCH OUT THIS IS DISCOURAGED!");
|
| 488 |
for(int i=0;i<obslimits.size();i++) explimits.push_back(obslimits[i]);
|
| 489 |
}
|
| 490 |
TH2F *bestlimit=(TH2F*)obslimits[0]->Clone("bestlimits");
|
| 491 |
TH2F *bestexplimit=(TH2F*)obslimits[0]->Clone("bestexplimits");
|
| 492 |
TH2F *bestlimitsource=(TH2F*)obslimits[0]->Clone("bestlimitsource");
|
| 493 |
TH2F *bestexp1plimit=(TH2F*)obslimits[0]->Clone("bestexp1plimit");
|
| 494 |
TH2F *bestexp1mlimit=(TH2F*)obslimits[0]->Clone("bestexp1mlimit");
|
| 495 |
TH2F *bestexp2plimit=(TH2F*)obslimits[0]->Clone("bestexp2plimit");
|
| 496 |
TH2F *bestexp2mlimit=(TH2F*)obslimits[0]->Clone("bestexp2mlimit");
|
| 497 |
|
| 498 |
for(int i=1;i<=explimits[0]->GetNbinsX();i++) {
|
| 499 |
for(int j=1;j<=explimits[0]->GetNbinsY();j++) {
|
| 500 |
float min=explimits[0]->GetBinContent(i,j);
|
| 501 |
float omin=obslimits[0]->GetBinContent(i,j);
|
| 502 |
int source=1;
|
| 503 |
for(int k=0;k<explimits.size();k++) {
|
| 504 |
float currlim=explimits[k]->GetBinContent(i,j);
|
| 505 |
if(currlim<min&&currlim>0) {
|
| 506 |
min=currlim;
|
| 507 |
omin=obslimits[k]->GetBinContent(i,j);
|
| 508 |
source=k+1;
|
| 509 |
}
|
| 510 |
}
|
| 511 |
if(min>0) {
|
| 512 |
bestlimitsource->SetBinContent(i,j,source);
|
| 513 |
bestlimit->SetBinContent(i,j,omin);
|
| 514 |
bestexplimit->SetBinContent(i,j,min);
|
| 515 |
if(ismSUGRA) {
|
| 516 |
bestexp1plimit->SetBinContent(i,j,exp1plimits[source-1]->GetBinContent(i,j));
|
| 517 |
bestexp1mlimit->SetBinContent(i,j,exp1mlimits[source-1]->GetBinContent(i,j));
|
| 518 |
bestexp2plimit->SetBinContent(i,j,exp2plimits[source-1]->GetBinContent(i,j));
|
| 519 |
bestexp2mlimit->SetBinContent(i,j,exp2mlimits[source-1]->GetBinContent(i,j));
|
| 520 |
}
|
| 521 |
}
|
| 522 |
}
|
| 523 |
}
|
| 524 |
gStyle->SetPadRightMargin(0.15);
|
| 525 |
TCanvas *canlimsource = new TCanvas("limsource","Source of best limits");
|
| 526 |
set_range(bestlimitsource,ismSUGRA,false);
|
| 527 |
bestlimitsource->SetTitle("Source of limit for best limits");
|
| 528 |
bestlimitsource->GetZaxis()->SetRangeUser(0,explimits.size()+1);
|
| 529 |
bestlimitsource->Draw("COLZ");
|
| 530 |
if(!ismSUGRA) bestlimitsource->Draw("TEXT,same");
|
| 531 |
TLegend *sourceleg = new TLegend(0.15,0.6,0.55,0.85);
|
| 532 |
for(int i=0;i<explimits.size();i++) {
|
| 533 |
stringstream legendentry;
|
| 534 |
legendentry << i+1 << " = " << give_nice_source_label(explimits[i]->GetName());
|
| 535 |
sourceleg->AddEntry(explimits[i], legendentry.str().c_str(),"");
|
| 536 |
}
|
| 537 |
sourceleg->SetLineColor(kWhite);sourceleg->SetFillColor(kWhite);
|
| 538 |
sourceleg->SetTextSize(0.03);
|
| 539 |
sourceleg->Draw("same");
|
| 540 |
DrawPrelim();
|
| 541 |
|
| 542 |
CompleteSave(canlimsource,"Limits/SourceOfBestLimits");
|
| 543 |
allbestexplimits.push_back(bestexplimit);
|
| 544 |
allbestexplimits.push_back(bestexp1plimit);
|
| 545 |
allbestexplimits.push_back(bestexp1mlimit);
|
| 546 |
allbestexplimits.push_back(bestexp2plimit);
|
| 547 |
allbestexplimits.push_back(bestexp2mlimit);
|
| 548 |
|
| 549 |
delete canlimsource;
|
| 550 |
delete bestexplimit;
|
| 551 |
delete bestlimitsource;
|
| 552 |
return bestlimit;
|
| 553 |
}
|
| 554 |
|
| 555 |
|
| 556 |
void create_exclusion_plots(vector<TH2F*> limits, bool ismSUGRA) {
|
| 557 |
TFile *xsecfile;
|
| 558 |
if(!ismSUGRA) {
|
| 559 |
xsecfile = new TFile(xsecfilename.c_str());
|
| 560 |
if(xsecfile->IsZombie()&&!ismSUGRA) {
|
| 561 |
write_error(__FUNCTION__,"Cross section file is invalid!!!!");
|
| 562 |
return;
|
| 563 |
}
|
| 564 |
xsecfile->Close();
|
| 565 |
}
|
| 566 |
if(!ismSUGRA) for(int i=0;i<limits.size();i++) limits[i]->Scale(1./0.19); // because for T5zz we forced one z to decay leptonically
|
| 567 |
|
| 568 |
vector<TH2F*> explimits;
|
| 569 |
vector<TH2F*> exp1mlimits;
|
| 570 |
vector<TH2F*> exp1plimits;
|
| 571 |
vector<TH2F*> exp2mlimits;
|
| 572 |
vector<TH2F*> exp2plimits;
|
| 573 |
vector<TH2F*> obslimits;
|
| 574 |
vector<TH2F*> flipmaps;
|
| 575 |
vector<TH2F*> crosssections;
|
| 576 |
|
| 577 |
cout << __LINE__ << endl;
|
| 578 |
for(int ilim=0;ilim<limits.size();ilim++) {
|
| 579 |
if(TString(limits[ilim]->GetName()).Contains("_explimitmap")) explimits.push_back(limits[ilim]);
|
| 580 |
if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_exp1plimitmap")) exp1plimits.push_back(limits[ilim]);
|
| 581 |
if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_exp1mlimitmap")) exp1mlimits.push_back(limits[ilim]);
|
| 582 |
if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_exp2plimitmap")) exp2plimits.push_back(limits[ilim]);
|
| 583 |
if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_exp2mlimitmap")) exp2mlimits.push_back(limits[ilim]);
|
| 584 |
if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_absolutecrosssectionmap")) crosssections.push_back(limits[ilim]);
|
| 585 |
if(TString(limits[ilim]->GetName()).Contains("_limitmap")) obslimits.push_back(limits[ilim]);
|
| 586 |
// if(TString(limits[ilim]->GetName()).Contains("_limitflipmap")) flipmaps.push_back(limits[ilim]);
|
| 587 |
}
|
| 588 |
|
| 589 |
TH2F *xsec;
|
| 590 |
if(!ismSUGRA) xsec = adjust_histo(get_XS(xsecfilename,"gluino",limits[0]),limits[0]);
|
| 591 |
vector<TH2F*> bestexplimits;
|
| 592 |
TH2F *bestlimits = make_best_limits(explimits,obslimits,ismSUGRA, exp1mlimits, exp1plimits, exp2mlimits, exp2plimits, bestexplimits);
|
| 593 |
bestlimits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
|
| 594 |
|
| 595 |
for(int ilim=0;ilim<limits.size();ilim++) {
|
| 596 |
limits[ilim]->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
|
| 597 |
}
|
| 598 |
|
| 599 |
if(!ismSUGRA) {
|
| 600 |
for(int ilim=0;ilim<obslimits.size();ilim++) make_SMS_exclusion(obslimits[ilim],xsec,ismSUGRA);
|
| 601 |
make_SMS_exclusion(bestlimits,xsec,ismSUGRA);
|
| 602 |
} else {
|
| 603 |
for(int ilim=0;ilim<obslimits.size();ilim++) {
|
| 604 |
draw_mSUGRA_exclusion(crosssections[0], obslimits[ilim], explimits[ilim], exp1mlimits[ilim], exp1plimits[ilim], exp2mlimits[ilim], exp2plimits[ilim]);
|
| 605 |
}
|
| 606 |
draw_mSUGRA_exclusion(crosssections[0], bestlimits, bestexplimits[0], bestexplimits[1], bestexplimits[2], bestexplimits[3], bestexplimits[4]);
|
| 607 |
}
|
| 608 |
delete bestlimits;
|
| 609 |
}
|
| 610 |
|
| 611 |
void decorate_mSUGRA(TH2F *cleanhisto, TVirtualPad *cvsSys,TGraph *expected,TGraph *expected2,TGraph *observed) {
|
| 612 |
cvsSys->SetRightMargin(standardmargin);
|
| 613 |
Int_t tanBeta_ = 10;
|
| 614 |
Bool_t plotLO_ = false;
|
| 615 |
|
| 616 |
//convert tanb value to string
|
| 617 |
std::stringstream tmp;
|
| 618 |
tmp << tanBeta_;
|
| 619 |
TString tanb( tmp.str() );
|
| 620 |
|
| 621 |
//set old exclusion Limits
|
| 622 |
TGraph* LEP_ch = set_lep_ch(tanBeta_);
|
| 623 |
TGraph* LEP_sl = set_lep_sl(tanBeta_);//slepton curve
|
| 624 |
TGraph* TEV_sg_cdf = set_tev_sg_cdf(tanBeta_);//squark gluino cdf
|
| 625 |
TGraph* TEV_sg_d0 = set_tev_sg_d0(tanBeta_);//squark gluino d0
|
| 626 |
TGraph* stau = set_tev_stau(tanBeta_);//stau
|
| 627 |
TGraph* NoEWSB = set_NoEWSB(tanBeta_);
|
| 628 |
|
| 629 |
TGraph* TEV_sn_d0_1 = set_sneutrino_d0_1(tanBeta_);
|
| 630 |
TGraph* TEV_sn_d0_2 = set_sneutrino_d0_2(tanBeta_);
|
| 631 |
|
| 632 |
//constant ssqquark and gluino lines
|
| 633 |
TF1* lnsq[10];
|
| 634 |
TF1* lngl[10];
|
| 635 |
|
| 636 |
TLatex* sq_text[10];
|
| 637 |
TLatex* gl_text[10];
|
| 638 |
|
| 639 |
for(int i = 0; i < 6; i++){
|
| 640 |
lnsq[i] = constant_squark(tanBeta_,i);
|
| 641 |
sq_text[i] = constant_squark_text(i,*lnsq[i],tanBeta_);
|
| 642 |
lngl[i] = constant_gluino(tanBeta_,i);
|
| 643 |
gl_text[i] = constant_gluino_text(i,*lngl[i]);
|
| 644 |
}
|
| 645 |
|
| 646 |
//Legends
|
| 647 |
TLegend* legst = makeStauLegend(0.05,tanBeta_);
|
| 648 |
TLegend* legNoEWSB = makeNoEWSBLegend(0.05,tanBeta_);
|
| 649 |
TLegend* legexp = makeExpLegend( *TEV_sg_cdf,*TEV_sg_d0,*LEP_ch,*LEP_sl,*TEV_sn_d0_1,0.035,tanBeta_);
|
| 650 |
|
| 651 |
TEV_sn_d0_1->SetLineWidth(1);
|
| 652 |
TEV_sn_d0_2->SetLineWidth(1);
|
| 653 |
TEV_sg_d0->SetLineWidth(1);
|
| 654 |
|
| 655 |
double m0min = 0;
|
| 656 |
if (tanBeta_ == 50) m0min=200;
|
| 657 |
TH2F* hist = new TH2F("h","h",100,m0min,1000,100,120,700);
|
| 658 |
hist->Draw();
|
| 659 |
hist->GetXaxis()->SetTitle("m_{0} [GeV]");
|
| 660 |
hist->GetXaxis()->CenterTitle();
|
| 661 |
hist->GetYaxis()->SetTitle("m_{1/2} [GeV]");
|
| 662 |
hist->GetYaxis()->CenterTitle();
|
| 663 |
hist->GetXaxis()->SetTitleSize(0.04);
|
| 664 |
hist->GetYaxis()->SetTitleSize(0.04);
|
| 665 |
hist->GetXaxis()->SetTitleOffset(1.2);
|
| 666 |
hist->GetYaxis()->SetTitleOffset(1.5);
|
| 667 |
hist->GetXaxis()->SetNdivisions(506);
|
| 668 |
hist->GetYaxis()->SetNdivisions(506);
|
| 669 |
|
| 670 |
int col[]={2,3,4};
|
| 671 |
|
| 672 |
TLegend* myleg;
|
| 673 |
|
| 674 |
if( plotLO_ ) myleg = new TLegend(0.3,0.65,0.65,0.8,NULL,"brNDC");
|
| 675 |
else myleg = new TLegend(0.25,0.76,0.44,0.91,NULL,"brNDC");
|
| 676 |
|
| 677 |
myleg->SetFillColor(0);
|
| 678 |
myleg->SetShadowColor(0);
|
| 679 |
myleg->SetTextSize(0.04);
|
| 680 |
myleg->SetBorderSize(0);
|
| 681 |
|
| 682 |
//constant squark and gluino mass contours
|
| 683 |
for (int it=0;it<5;it++) {
|
| 684 |
lngl[it]->Draw("same");
|
| 685 |
lnsq[it]->Draw("same");
|
| 686 |
sq_text[it]->Draw();
|
| 687 |
gl_text[it]->Draw();
|
| 688 |
}
|
| 689 |
|
| 690 |
//exclusion limits previous experiments
|
| 691 |
if(tanBeta_ == 3){
|
| 692 |
TEV_sn_d0_1->Draw("fsame");
|
| 693 |
TEV_sn_d0_2->Draw("fsame");
|
| 694 |
}
|
| 695 |
LEP_ch->Draw("fsame");
|
| 696 |
if (tanBeta_ != 50) LEP_sl->Draw("fsame");
|
| 697 |
|
| 698 |
//remove CDF/D0 excluded regions
|
| 699 |
TEV_sg_cdf->Draw("fsame");
|
| 700 |
TEV_sg_d0->Draw("same");
|
| 701 |
TEV_sg_d0->Draw("fsame");
|
| 702 |
|
| 703 |
//other labels
|
| 704 |
Double_t xpos = 0;
|
| 705 |
Double_t xposi = 0;
|
| 706 |
Double_t ypos = 0;
|
| 707 |
if(tanBeta_ == 50) xposi = 100;
|
| 708 |
if(tanBeta_ == 50) xpos = 200;
|
| 709 |
if(tanBeta_ == 50) ypos = -10;
|
| 710 |
|
| 711 |
TString text_tanBeta;
|
| 712 |
text_tanBeta = "tan#beta = "+tanb+", A_{0} = 0, #mu > 0";
|
| 713 |
TLatex* cmssmpars = new TLatex(/*530.+xpos,690.+ypos-130*/150,650,text_tanBeta);
|
| 714 |
|
| 715 |
cmssmpars->SetTextSize(0.03);
|
| 716 |
cmssmpars->Draw("same");
|
| 717 |
|
| 718 |
TLatex* lep_chargino = new TLatex(250,135,"LEP2 #tilde{#chi}_{1}^{#pm}");
|
| 719 |
lep_chargino->SetTextSize(0.03);
|
| 720 |
lep_chargino->SetTextFont(42);
|
| 721 |
|
| 722 |
TLatex* lep_slepton = new TLatex(26,190,"LEP2 #tilde{#font[12]{l}}^{#pm}");
|
| 723 |
lep_slepton->SetTextSize(0.03);
|
| 724 |
lep_slepton->SetTextAngle(-83);
|
| 725 |
lep_slepton->SetTextFont(42);
|
| 726 |
|
| 727 |
if(draw2sigma) expected2->Draw("f");
|
| 728 |
expected->SetLineColor(expected->GetFillColor());
|
| 729 |
expected2->SetLineColor(expected2->GetFillColor());
|
| 730 |
expected->Draw("f");
|
| 731 |
observed->Draw("c");
|
| 732 |
|
| 733 |
stau->Draw("fsame");
|
| 734 |
NoEWSB->Draw("fsame");
|
| 735 |
|
| 736 |
legexp->AddEntry(observed,"Observed limit","l");
|
| 737 |
legexp->AddEntry(expected,"Expected 1#sigma limit","f");
|
| 738 |
if(draw2sigma) legexp->AddEntry(expected2,"Expected 2#sigma limit","f");
|
| 739 |
legexp->SetY1(0.60);
|
| 740 |
legexp->SetX1(0.55);
|
| 741 |
legexp->Draw();
|
| 742 |
legst->Draw();
|
| 743 |
|
| 744 |
hist->Draw("sameaxis");
|
| 745 |
cvsSys->RedrawAxis();
|
| 746 |
cvsSys->Update();
|
| 747 |
DrawPrelim();
|
| 748 |
}
|
| 749 |
|
| 750 |
void smooth_line_once(TGraph *gr) {
|
| 751 |
//going to smooth graph
|
| 752 |
//need to take into account the DIRECTION we're heading so we noticed for expected shape when we turn around and go back
|
| 753 |
float previousX=-100;
|
| 754 |
int sign=1;
|
| 755 |
for(int i=0;i<gr->GetN();i++) {
|
| 756 |
Double_t x,y,x1,y1,x2,y2;
|
| 757 |
bool turning=false;
|
| 758 |
gr->GetPoint(i,x,y);
|
| 759 |
gr->GetPoint(i+1,x1,y1);//need to handle exception here
|
| 760 |
gr->GetPoint(i-1,x2,y2);//need to handle exception here
|
| 761 |
if(i!=gr->GetN()-1&& (sign*x<=sign*previousX||sign*x1<=sign*x)) {
|
| 762 |
//turned around!
|
| 763 |
sign=(-1)*sign;
|
| 764 |
// do NOT do any smoothing on this point!
|
| 765 |
} else {
|
| 766 |
float newval=y;
|
| 767 |
if(i>0&&i<(gr->GetN())-1) newval=(1.0/3.0)*(y+y1+y2);
|
| 768 |
if(i==0) newval=(1.0/2.0)*(y+y1);
|
| 769 |
if(i==gr->GetN()-1) newval=(1.0/2.0)*(y+y2);
|
| 770 |
gr->SetPoint(i,x,newval);
|
| 771 |
}
|
| 772 |
previousX=x;
|
| 773 |
}
|
| 774 |
}
|
| 775 |
|
| 776 |
void smooth_line(TGraph *gr) {
|
| 777 |
smooth_line_once(gr);
|
| 778 |
smooth_line_once(gr);
|
| 779 |
smooth_line_once(gr);
|
| 780 |
}
|
| 781 |
void set_range(TH2F *histo, bool ismSUGRA, bool pushoutyz=false) {
|
| 782 |
if(ismSUGRA) {
|
| 783 |
histo->GetXaxis()->SetRangeUser(0,2000);
|
| 784 |
histo->GetYaxis()->SetRangeUser(0,800);
|
| 785 |
histo->GetXaxis()->SetTitle("m_{0} [GeV]");
|
| 786 |
histo->GetYaxis()->SetTitle("m_{1/2} [GeV]");
|
| 787 |
histo->GetXaxis()->SetRangeUser(0,2000);
|
| 788 |
histo->GetXaxis()->SetNdivisions(504,true);
|
| 789 |
} else {
|
| 790 |
histo->GetXaxis()->SetRangeUser(50,1200);
|
| 791 |
histo->GetYaxis()->SetRangeUser(50,1200);
|
| 792 |
histo->GetXaxis()->SetTitle("m_{#tilde{g}} [GeV]");
|
| 793 |
histo->GetYaxis()->SetTitle("m_{LSP} [GeV]");
|
| 794 |
histo->GetXaxis()->SetTitleSize(0.04);
|
| 795 |
histo->GetYaxis()->SetTitleSize(0.04);
|
| 796 |
histo->GetZaxis()->SetTitleSize(0.04);
|
| 797 |
histo->GetXaxis()->SetTitleOffset(1.2);
|
| 798 |
histo->GetYaxis()->SetTitleOffset(1.5);
|
| 799 |
// if(pushoutyz) histo->GetYaxis()->SetTitleOffset(1.6);
|
| 800 |
if(pushoutyz) histo->GetZaxis()->SetTitleOffset(1.6);
|
| 801 |
// histo->GetYaxis()->SetTitleOffset(1.2);
|
| 802 |
// histo->GetZaxis()->SetTitleOffset(1.35);
|
| 803 |
|
| 804 |
|
| 805 |
}
|
| 806 |
histo->GetXaxis()->CenterTitle();
|
| 807 |
histo->GetYaxis()->CenterTitle();
|
| 808 |
histo->SetStats(0);
|
| 809 |
}
|
| 810 |
|
| 811 |
TH2F* adjust_histo(TH2F *oldhist, TH2F *refhisto) {
|
| 812 |
TH2F *histo = new TH2F(refhisto->GetName(),refhisto->GetName(),
|
| 813 |
refhisto->GetNbinsX(),refhisto->GetXaxis()->GetBinLowEdge(1),refhisto->GetXaxis()->GetBinLowEdge(refhisto->GetNbinsX())+refhisto->GetXaxis()->GetBinWidth(refhisto->GetNbinsX()),
|
| 814 |
refhisto->GetNbinsY(),refhisto->GetYaxis()->GetBinLowEdge(1),refhisto->GetYaxis()->GetBinLowEdge(refhisto->GetNbinsY())+refhisto->GetYaxis()->GetBinWidth(refhisto->GetNbinsY()));
|
| 815 |
|
| 816 |
for(int ix=1;ix<=refhisto->GetNbinsX();ix++) {
|
| 817 |
for(int iy=1;iy<=refhisto->GetNbinsX();iy++) {
|
| 818 |
int binnum=refhisto->FindBin(refhisto->GetXaxis()->GetBinCenter(ix),refhisto->GetYaxis()->GetBinCenter(iy));
|
| 819 |
histo->SetBinContent(binnum,oldhist->GetBinContent(ix,iy));
|
| 820 |
}
|
| 821 |
}
|
| 822 |
return histo;
|
| 823 |
}
|
| 824 |
|
| 825 |
void process_syst_plot(TH2F *rhisto,string saveto,bool ismSUGRA) {
|
| 826 |
TH2F *histo = prep_histo(rhisto,ismSUGRA); // this is to be independent of the style used at creation time
|
| 827 |
float rightmargin=gStyle->GetPadRightMargin();
|
| 828 |
gStyle->SetPadRightMargin(0.20);
|
| 829 |
TString name = rhisto->GetName();
|
| 830 |
if(name.Contains("Nevents")) gStyle->SetPadRightMargin(0.22);
|
| 831 |
TCanvas *can = new TCanvas("syst_plot","Systematics Plot");
|
| 832 |
set_range(histo,ismSUGRA,true);
|
| 833 |
|
| 834 |
|
| 835 |
if(name.Contains("efficiency")) {
|
| 836 |
histo->GetZaxis()->SetTitle("A #times #varepsilon (#geq 1 Z(ll))");
|
| 837 |
// histo->GetZaxis()->SetRangeUser(0.0,0.2);
|
| 838 |
}
|
| 839 |
if(name.Contains("Nevents")) {
|
| 840 |
histo->GetZaxis()->SetTitle("N(events)");
|
| 841 |
histo->GetZaxis()->SetTitleOffset(histo->GetZaxis()->GetTitleOffset()+0.4);
|
| 842 |
}
|
| 843 |
if(name.Contains("sysjes")) {
|
| 844 |
histo->GetZaxis()->SetTitle("Jet Energy Scale");
|
| 845 |
histo->GetZaxis()->SetRangeUser(0.0,0.2);
|
| 846 |
}
|
| 847 |
if(name.Contains("sysjsu")) {
|
| 848 |
histo->GetZaxis()->SetTitle("JZB Scale Uncertainty");
|
| 849 |
histo->GetZaxis()->SetRangeUser(0.0,0.5);
|
| 850 |
}
|
| 851 |
if(name.Contains("sysresmap")) {
|
| 852 |
histo->GetZaxis()->SetTitle("Resulution");
|
| 853 |
histo->GetZaxis()->SetRangeUser(0.0,0.5);
|
| 854 |
}
|
| 855 |
if(name.Contains("sysstatmap")) {
|
| 856 |
histo->GetZaxis()->SetTitle("Statistical Error");
|
| 857 |
histo->GetZaxis()->SetRangeUser(0.0,0.01);
|
| 858 |
}
|
| 859 |
if(name.Contains("systotmap")) {
|
| 860 |
histo->GetZaxis()->SetTitle("All Systematic Errors");
|
| 861 |
histo->GetZaxis()->SetRangeUser(0.0,0.5);
|
| 862 |
}
|
| 863 |
|
| 864 |
histo->GetZaxis()->CenterTitle();
|
| 865 |
gStyle->SetStripDecimals(false);
|
| 866 |
histo->GetXaxis()->SetDecimals(true);
|
| 867 |
|
| 868 |
histo->Draw("COLZ");
|
| 869 |
DrawPrelim();
|
| 870 |
|
| 871 |
CompleteSave(can,(saveto+(string)histo->GetName()));
|
| 872 |
|
| 873 |
gStyle->SetPadRightMargin(rightmargin);
|
| 874 |
|
| 875 |
delete can;
|
| 876 |
}
|
| 877 |
|
| 878 |
void make_all_syst_plots(vector<TH2F*> all_histos,bool ismSUGRA) {
|
| 879 |
string saveto="Systematics/";
|
| 880 |
for(int iplot=0;iplot<all_histos.size();iplot++) {
|
| 881 |
process_syst_plot(all_histos[iplot],saveto,ismSUGRA);
|
| 882 |
}
|
| 883 |
}
|
| 884 |
|
| 885 |
void process_file(TFile* file, float stdmargin) {
|
| 886 |
standardmargin=stdmargin;
|
| 887 |
xsecfilename="reference_xSec_SMS.root";
|
| 888 |
|
| 889 |
// can receive a file with systematics and limits mixed, or a file with systematics only , or a file with limits only.
|
| 890 |
TIter nextkey(file->GetListOfKeys());
|
| 891 |
TKey *key;
|
| 892 |
|
| 893 |
bool ismSUGRA=false;
|
| 894 |
|
| 895 |
vector<TH2F*> systematics_histos;
|
| 896 |
vector<TH2F*> limits_histos;
|
| 897 |
while ((key = (TKey*)nextkey()))
|
| 898 |
{
|
| 899 |
TObject *obj = key->ReadObj();
|
| 900 |
if (!(obj->IsA()->InheritsFrom("TH2"))) continue;
|
| 901 |
TString name=(TString)(obj->GetName());
|
| 902 |
bool is_limit=false;
|
| 903 |
bool is_systematic=false;
|
| 904 |
|
| 905 |
if(name.Contains("limitmap")) is_limit=true;
|
| 906 |
if(name.Contains("explimitmap")) is_limit=true;
|
| 907 |
if(name.Contains("exp1plimitmap")) is_limit=true;
|
| 908 |
if(name.Contains("exp2plimitmap")) is_limit=true;
|
| 909 |
if(name.Contains("exp1mlimitmap")) is_limit=true;
|
| 910 |
if(name.Contains("exp2mlimitmap")) is_limit=true;
|
| 911 |
if(name.Contains("exclusionmap")) is_limit=true;
|
| 912 |
if(name.Contains("crosssectionmap")) is_limit=true;
|
| 913 |
if(name.Contains("absolutecrosssectionmap")) is_limit=true;
|
| 914 |
if(name.Contains("limitflipmap")) is_limit=true;
|
| 915 |
|
| 916 |
if(name.Contains("sysjes")) is_systematic=true;
|
| 917 |
if(name.Contains("sysjsu")) is_systematic=true;
|
| 918 |
if(name.Contains("sysresmap")) is_systematic=true;
|
| 919 |
if(name.Contains("efficiencymap")) is_systematic=true;
|
| 920 |
if(name.Contains("efficiencymapAcc")) is_systematic=true;
|
| 921 |
if(name.Contains("efficiencymapID")) is_systematic=true;
|
| 922 |
if(name.Contains("efficiencymapJets")) is_systematic=true;
|
| 923 |
if(name.Contains("efficiencymapSignal")) is_systematic=true;
|
| 924 |
if(name.Contains("noscefficiencymap")) is_systematic=true;
|
| 925 |
if(name.Contains("Neventsmap")) is_systematic=true;
|
| 926 |
if(name.Contains("ipointmap")) is_systematic=true;
|
| 927 |
if(name.Contains("syspdfmap")) is_systematic=true;
|
| 928 |
if(name.Contains("systotmap")) is_systematic=true;
|
| 929 |
if(name.Contains("sysstatmap")) is_systematic=true;
|
| 930 |
if(name.Contains("timemap")) is_systematic=true;
|
| 931 |
if(name.Contains("flipmap")) is_systematic=true;
|
| 932 |
|
| 933 |
if(is_limit) limits_histos.push_back((TH2F*) obj);
|
| 934 |
else if(is_systematic) systematics_histos.push_back((TH2F*) obj);
|
| 935 |
if(name.Contains("mSUGRA")) ismSUGRA=true;
|
| 936 |
}
|
| 937 |
if(systematics_histos.size()>0) make_all_syst_plots(systematics_histos,ismSUGRA);
|
| 938 |
if(limits_histos.size()>0) create_exclusion_plots(limits_histos,ismSUGRA);
|
| 939 |
}
|
| 940 |
|
