ra4_macros  bede988c286599a3a84b77a4d788ac0a971e89f9
plot_trigger.cxx
Go to the documentation of this file.
1 // plot_trigger: Plots for the trigger section of the RA4 note
2 
3 #include <iostream>
4 #include <fstream>
5 #include <vector>
6 #include <cmath>
7 
8 #include "TChain.h"
9 #include "TH1D.h"
10 #include "TF1.h"
11 #include "TMath.h"
12 #include "TCanvas.h"
13 #include "TPad.h"
14 #include "TLegend.h"
15 #include "TLine.h"
16 #include "TGaxis.h"
17 #include "TBox.h"
18 #include "TLatex.h"
19 #include "TPaveText.h"
20 #include "TString.h"
21 #include "TStyle.h"
22 #include "TGraphAsymmErrors.h"
23 #include "TEfficiency.h"
24 
25 #include "styles.hpp"
26 #include "utilities.hpp"
27 #include "utilities_macros.hpp"
28 
29 namespace {
30  bool do_rund = false;
31  bool do_dps = false;
32  bool do_rpv = true;
33  TString plot_type = ".pdf";
34 }
35 
36 using namespace std;
37 
38 Double_t errorFun(Double_t *x, Double_t *par) {
39  double value(0.5*par[0]*(1. + TMath::Erf( (x[0] - par[1]) / (sqrt(2.)*par[2]) )));
40  return value;
41 }
42 
43 void PlotTurnOn(TChain *data, TString var, int nbins, double minx, double maxx, TString xtitle,
44  TString den, TString num, TString title="", TString ytitle="", float minfit=-1., bool isData=true);
45 TString Efficiency(TChain *data, TString den, TString num, float &effic, float &errup, float &errdown);
46 
47 int main(){
48 
49  styles style("HLTStyle"); style.setDefaultStyle();
50  gStyle->SetPadTickY(0);
51 
52  TString folder("/cms2r0/babymaker/babies/2015_09_28/skim_json/");
53 
54  TChain c_htmht("tree"); c_htmht.Add(folder+"*HTMHT*root");
55  TChain c_jetht("tree"); c_jetht.Add(folder+"*JetHT*root");
56  TChain c_met("tree"); c_met.Add(folder+"*MET*root");
57  TChain c_mu("tree"); c_mu.Add(folder+"*SingleMuon*");
58  TChain c_el("tree"); c_el.Add(folder+"*SingleElectron*");
59  TChain c_lep("tree"); c_lep.Add(folder+"singlelep/*Single*");
60  TChain c_all("tree"); c_all.Add(folder+"alldata/*root");
61  TChain c_had("tree"); c_had.Add(folder+"hadronic/*root");
62  TChain c_jetht_rpv("tree");
63  c_jetht_rpv.Add("/homes/cawest/links/JetHT_Run2015C_25ns-05Oct2015-v1/*root");
64  c_jetht_rpv.Add("/homes/cawest/links/JetHT_Run2015D-05Oct2015-v1/*root");
65  c_jetht_rpv.Add("/homes/cawest/links/JetHT_Run2015D-PromptReco-v4/*root");
66  TChain c_mu_rpv("tree");
67  c_mu_rpv.Add("/homes/cawest/links/SingleMuon_Run2015C-23Sep2015-v1/*root");
68  c_mu_rpv.Add("/homes/cawest/links/SingleMuon_Run2015D-05Oct2015-v1/*root");
69  c_mu_rpv.Add("/homes/cawest/links/SingleMuon_Run2015D-PromptReco-v4/*root");
70 
71  if(do_rpv){
72  float lmin(500), lmax(2000);
73  int lbins(static_cast<int>((lmax-lmin)/25));
74  PlotTurnOn(&c_mu_rpv, "ht", lbins,lmin,lmax, "H_{T}",
75  "(trig[20]&&nmus==1&&mus_pt[0]>35)", "trig[12]","HLT_IsoMu27 && p_{T,#mu}>35 GeV", "PFHT800",700);
76  PlotTurnOn(&c_jetht_rpv, "ht", lbins,lmin,lmax, "H_{T}",
77  "(trig[11])", "trig[12]","PFHT450", "PFHT800",800);
78  }
79 
80  if(do_rund){
81  float lmin(25), lmax(300);
82  int lbins(static_cast<int>((lmax-lmin)/12.5));
83 
84  lmin = 10; lmax = 80; lbins = static_cast<int>((lmax-lmin)/2.5);
85  PlotTurnOn(&c_htmht, "Max$(mus_pt*(mus_sigid&&mus_miniso<0.2))-0.1", lbins,lmin,lmax, "#mu_{medium} p_{T}",
86  "(trig[0])&&met>150", "trig[4]","HT350_MET100, MET > 150", "Mu15_HT350");
87  PlotTurnOn(&c_htmht, "Max$(mus_pt*(mus_sigid&&mus_miniso<0.2))-0.1", lbins,lmin,lmax, "#mu_{medium} p_{T}",
88  "(trig[0])&&met>150", "trig[17]","HT350_MET100, MET > 150", "IsoMu18");
89 
90  PlotTurnOn(&c_htmht, "Max$(els_pt*(els_sigid&&els_miniso<0.1))-0.1", lbins,lmin,lmax, "e_{medium} p_{T}",
91  "(trig[0])&&met>150", "trig[8]","HT350_MET100, MET > 150", "Ele15_HT350");
92  PlotTurnOn(&c_htmht, "Max$(els_pt*(els_sigid&&els_miniso<0.1))-0.1", lbins,lmin,lmax, "e_{medium} p_{T}",
93  "(trig[0])&&met>150", "trig[23]","HT350_MET100, MET > 150", "Ele23_WPLoose");
94 
95  float htmin(175), htmax(850);
96  int htbins(static_cast<int>((htmax-htmin)/12.5));
97  htmin = 175; htmax = 850; htbins = static_cast<int>((htmax-htmin)/12.5);
98  PlotTurnOn(&c_met, "ht", htbins,htmin,htmax, "H_{T}",
99  "trig[14]&&nvleps>=1", "trig[0]","MET170, n_{l} #geq 1", "HT350_MET100", 300);
100  htmin = 100; htmax = 850; htbins = static_cast<int>((htmax-htmin)/25);
101  PlotTurnOn(&c_met, "ht", htbins,htmin,htmax, "H_{T}",
102  "(trig[13])&&(nmus)>=1&&onmaxmu>15", "trig[1]","MET120_Mu5, n_{#mu} #geq 1",
103  "Mu15_HT350_MET50", 200);
104 
105  float metmin(0), metmax(460);
106  int metbins(static_cast<int>((metmax-metmin)/20));
107  PlotTurnOn(&c_el, "met", metbins,metmin,metmax, "E_{T}^{miss}",
108  "(trig[8])&&njets_ra2>=4&&ht_ra2>450&&ht_ra2<750", "trig[5]",
109  "Ele15_HT350, n_{j} #geq 4, 450 < H_{T} < 750", "Ele15_HT350_MET50");
110  PlotTurnOn(&c_el, "met", metbins,metmin,metmax, "E_{T}^{miss}",
111  "(trig[8])&&njets_ra2>=4&&ht_ra2>750", "trig[5]",
112  "Ele15_HT350, n_{j} #geq 4, H_{T} > 750", "Ele15_HT350_MET50");
113  PlotTurnOn(&c_el, "met", metbins,metmin,metmax, "E_{T}^{miss}",
114  "(trig[8])&&njets_ra2>=4&&ht_ra2>450", "trig[5]",
115  "Ele15_HT350, n_{j} #geq 4, H_{T} > 450", "Ele15_HT350_MET50");
116 
117  } // if(do_rund)
118 
119  if(do_dps){
120  float htmin(175), htmax(850);
121  int htbins(static_cast<int>((htmax-htmin)/12.5));
122  htmin = 175; htmax = 850; htbins = static_cast<int>((htmax-htmin)/12.5);
123  PlotTurnOn(&c_met, "ht", htbins,htmin,htmax, "H_{T}",
124  "trig[14]&&nvleps==0", "trig[0]","MET170, no leptons", "HT350_MET100", 300);
125  htmin = 100; htmax = 850; htbins = static_cast<int>((htmax-htmin)/25);
126  PlotTurnOn(&c_met, "ht", htbins,htmin,htmax, "H_{T}",
127  "(trig[13])&&(nmus)>=1&&onmaxmu>15", "trig[1]","MET120_Mu5, n_{#mu} #geq 1",
128  "Mu15_HT350_MET50", 200);
129 
130  htmin = 575; htmax = 1450; htbins = static_cast<int>((htmax-htmin)/25);
131  PlotTurnOn(&c_jetht, "ht", htbins,htmin,htmax, "H_{T}",
132  "trig[11]&&nvleps==0", "trig[12]","HT475, no leptons", "HT800");
133  htmin = 375; htmax = 1250; htbins = static_cast<int>((htmax-htmin)/25);
134  PlotTurnOn(&c_el, "ht", htbins,htmin,htmax, "H_{T}",
135  "(trig[5])", "trig[6]","Ele15_HT350_MET50", "Ele15_HT600", 550);
136  PlotTurnOn(&c_el, "ht", htbins,htmin,htmax, "H_{T}",
137  "(trig[23])&&nels>=1&&onmaxel>15", "trig[6]","Ele23", "Ele15_HT600", 550);
138 
139  htmin = 175; htmax = 850; htbins = static_cast<int>((htmax-htmin)/12.5);
140  PlotTurnOn(&c_met, "ht_hlt", htbins,htmin,htmax, "H_{T}",
141  "trig[14]&&nvleps==0", "trig[0]","MET170, no leptons", "HT350_MET100", 300);
142  PlotTurnOn(&c_met, "ht_hlt", htbins,htmin,htmax, "H_{T}",
143  "(trig[13])&&(nmus)>=1&&onmaxmu>15", "trig[1]","MET120_Mu5, n_{#mu} #geq 1",
144  "Mu15_HT350_MET50", 300);
145 
146  htmin = 575; htmax = 1450; htbins = static_cast<int>((htmax-htmin)/25);
147  PlotTurnOn(&c_jetht, "ht_hlt", htbins,htmin,htmax, "H_{T}",
148  "trig[11]&&nvleps==0", "trig[12]","HT475, no leptons", "HT800");
149  htmin = 375; htmax = 1250; htbins = static_cast<int>((htmax-htmin)/25);
150  PlotTurnOn(&c_el, "ht_hlt", htbins,htmin,htmax, "H_{T}",
151  "(trig[5])", "trig[6]","Ele15_HT350_MET50", "Ele15_HT600", 550);
152  PlotTurnOn(&c_el, "ht_hlt", htbins,htmin,htmax, "H_{T}",
153  "(trig[23])&&nels>=1&&onmaxel>15", "trig[6]","Ele23", "Ele15_HT600", 550);
154 
155  float metmin(0), metmax(460);
156  int metbins(static_cast<int>((metmax-metmin)/20));
157  PlotTurnOn(&c_el, "met", metbins,metmin,metmax, "E_{T}^{miss}",
158  "(trig[8])&&njets_ra2>=4&&ht_ra2>450&&ht_ra2<750", "trig[5]",
159  "Ele15_HT350, n_{j} #geq 4, 450 < H_{T} < 750", "Ele15_HT350_MET50");
160  PlotTurnOn(&c_el, "met", metbins,metmin,metmax, "E_{T}^{miss}",
161  "(trig[8])&&njets_ra2>=4&&ht_ra2>750", "trig[5]",
162  "Ele15_HT350, n_{j} #geq 4, H_{T} > 750", "Ele15_HT350_MET50");
163  PlotTurnOn(&c_el, "met", metbins,metmin,metmax, "E_{T}^{miss}",
164  "(trig[8])&&njets_ra2>=4&&ht_ra2>450", "trig[5]",
165  "Ele15_HT350, n_{j} #geq 4, H_{T} > 450", "Ele15_HT350_MET50");
166  PlotTurnOn(&c_el, "met", metbins,metmin,metmax, "E_{T}^{miss}",
167  "(trig[8])&&njets_ra2>=4&&ht_ra2>450", "trig[0]",
168  "Ele15_HT350, n_{j} #geq 4, H_{T} > 450", "HT350_MET100");
169  PlotTurnOn(&c_mu, "met", metbins,metmin,metmax, "E_{T}^{miss}",
170  "(trig[4])&&njets_ra2>=4&&ht_ra2>450", "trig[0]",
171  "Mu15_HT350, n_{j} #geq 4, H_{T} > 450", "HT350_MET100");
172 
173  PlotTurnOn(&c_el, "mht_ra2", metbins,metmin,metmax, "H_{T}^{miss}",
174  "(trig[8])&&njets_ra2>=4&&ht_ra2>450", "trig[0]",
175  "Ele15_HT350, n_{j} #geq 4, H_{T} > 450", "HT350_MET100");
176  PlotTurnOn(&c_el, "mht_ra2", metbins,metmin,metmax, "H_{T}^{miss}",
177  "(trig[8])&&njets_ra2>=4&&ht_ra2>450&&mht_ra2/met>0.5&&mht_ra2/met<2", "trig[0]",
178  "Ele15_HT350, n_{j} #geq 4, H_{T} > 450", "HT350_MET100");
179 
180 
181  }
182 
183 }
184 
185 void PlotTurnOn(TChain *data, TString var, int nbins, double minx, double maxx, TString xtitle,
186  TString den, TString num, TString title, TString ytitle, float minfit, bool isData){
187  styles style("HLTStyle"); gStyle->SetPadTickY(0);
188  bool dofit(minfit>=-1);
189  TCanvas can;
190  //can.SetGrid();
191  TH1D* histo[2];
192  TString hname, totCut, pname;
193  // den = "("+den+")&&json&&pass";
194  hname = "den"; totCut = den;
195  histo[0] = new TH1D(hname, "", nbins, minx, maxx);
196  data->Project(hname, var, totCut);
197 
198  hname = "num"; totCut = "("+den+")&&("+num+")";
199  histo[1] = new TH1D(hname, "", nbins, minx, maxx);
200  data->Project(hname, var, totCut);
201 
202  // Adding overflow bins
203  for(unsigned his(0); his<2; his++)
204  histo[his]->SetBinContent(nbins, histo[his]->GetBinContent(nbins)+histo[his]->GetBinContent(nbins+1));
205 
206  TGraphAsymmErrors heff(histo[1], histo[0], "cp");
207  //TEfficiency heff(*histo[1], *histo[0]);
208  heff.SetMarkerStyle(20); heff.SetMarkerSize(0.9);
209  heff.SetMarkerColor(1); heff.SetLineColor(1);
210 
211  TString epsi("#scale[1.3]{#font[122]{e}}");
212  //epsi = "Efficiency";
213  // Ploting denominator
214  float hscaled(0.3), maxeff(1.42);
215  float hfactor(hscaled/histo[1]->GetMaximum()), hmax(histo[1]->GetMaximum());
216  float axismax(hmax*maxeff/hscaled);
217  histo[1]->Scale(hfactor);
218  //histo[1]->SetFillStyle(3344);
219  histo[1]->SetFillColor(kBlue-9);
220  histo[1]->SetLineStyle(0);
221  //histo[1]->SetTitle("Denom: "+title);
222  histo[1]->GetXaxis()->SetTitle("Offline "+xtitle+" [GeV]");
223  histo[1]->GetYaxis()->SetTitle(epsi+" ["+ytitle+"]");
224  histo[1]->GetYaxis()->SetTitle("Efficiency ["+ytitle+"]");
225  histo[1]->GetYaxis()->SetRangeUser(0,maxeff);
226  histo[1]->GetYaxis()->CenterTitle(true);
227  histo[1]->Draw();
228 
229  TLine line; line.SetLineStyle(3);
230  line.DrawLine(minx, 1, maxx, 1);
231 
232  histo[0]->Scale(hfactor);
233  histo[0]->SetLineColor(kBlue+2);
234  histo[0]->SetLineStyle(2);
235  histo[0]->SetLineWidth(2);
236  histo[0]->Draw("same");
237 
238  pname = "plots/turnon_"+format_tag(var)+"_";
239  pname += maxx; pname += "_"+format_tag(num)+"_"+format_tag(den);
240  if(minfit>0) {pname += "_min"; pname += minfit; }
241  if(isData) pname += "_data";
242  else pname += "_mc";
243  pname += plot_type;
244  pname.ReplaceAll("_json_pass","");
245 
246 
247  // Fitting turn on curve
248  if(minfit<0 && dofit) minfit = minx;
249  TF1 *fitCurve = new TF1("fitCurve",errorFun,minx,maxx,3);
250  Double_t params[] = {0.99,minx+(maxx-minx)/4.,50., 50.};
251  //Double_t params[] = {0.99,880,50.};
252  fitCurve->SetParameters(params);
253  fitCurve->SetParNames("#epsilon","#mu","#sigma");
254  fitCurve->SetLineColor(2);
255  fitCurve->SetLineWidth(2);
256  if(dofit) heff.Fit("fitCurve","QEM+","",minfit, maxx);
257 
258  heff.Draw("p");
259  histo[1]->Draw("axis same");
260 
261  float binw((maxx-minx)/static_cast<float>(nbins));
262  int digits((binw-floor(binw))>0?1:0);
263 
264  TString ntitle("Events/("+RoundNumber(maxx-minx,digits,nbins)+" GeV)");
265  TGaxis *axis = new TGaxis(maxx,0, maxx, maxeff,0,axismax,508,"+L");
266  axis->SetLineColor(kBlue+2);
267  axis->SetTextColor(kBlue+2); axis->SetLabelColor(kBlue+2);
268  axis->SetTextFont(style.nFont); axis->SetLabelFont(style.nFont);
269  axis->SetTitleSize(style.LabelSize); axis->SetLabelSize(style.LabelSize);
270  if(axismax>=1000) axis->SetTitleOffset(style.yTitleOffset+0.34);
271  else axis->SetTitleOffset(style.yTitleOffset+0.22);
272  axis->SetTitle(ntitle); axis->CenterTitle(true);
273  axis->Draw();
274 
275 
276  float effic, errup, errdown;
277  float mu(fitCurve->GetParameter(1)), sigma(fitCurve->GetParameter(2));
278  float rgev(mu>200?10:5);
279  float var_plateau_f(floor((mu+3*sigma+5)/rgev)*rgev);
280  if(!dofit) var_plateau_f = fabs(minfit);
281  TString den_plateau(den), var_plateau(RoundNumber(var_plateau_f, 0));
282  den_plateau += ("&&"+var+">"+var_plateau);
283  Efficiency(data, den_plateau, num, effic, errup, errdown);
284 
285  // 98th percentile of Gaussian from Wolfram Alpha
286  float p98(fitCurve->GetParameter(1)+2.05*fitCurve->GetParameter(2));
287  float eplateau(fitCurve->GetParError(0)*100);
288  if(eplateau<0.1 && dofit) cout<<"Error on plateau "<<eplateau<<"%"<<endl;
289  epsi.ToLower();
290  // TString fitpar("Plateau "+epsi+" = ("+RoundNumber(fitCurve->GetParameter(0)*100,1)+" #pm "+
291  // RoundNumber(eplateau,1)+")%");
292  // TString fitpar("Plateau "+epsi+" = "+RoundNumber(effic*100,1)+"^{+"+RoundNumber(errup*100,1)+
293  // "}_{-"+RoundNumber(errdown*100,1)+"} %");
294  TString fitpar(epsi+"("+xtitle+" > "+var_plateau+" GeV) = "+RoundNumber(effic*100,1)+"^{+"+RoundNumber(errup*100,1)+
295  "}_{-"+RoundNumber(errdown*100,1)+"} %");
296  TLatex label; label.SetTextSize(0.042);
297  label.SetTextAlign(33); //label.SetNDC();
298  float range(maxx-minx);
299  float x2(maxx-0.04*range), y2(maxeff-0.07), ysingle(0.1);
300  label.DrawLatex(x2, y2, "Denom: #font[52]{"+title+"}");
301  label.DrawLatex(x2, y2-ysingle,fitpar);
302  fitpar = "98% of plateau at "+RoundNumber(p98,0)+" GeV";
303  if(dofit) label.DrawLatex(x2, y2-2.3*ysingle,fitpar);
304 
305  // Drawing CMS preliminary
306  label.SetNDC(); label.SetTextAlign(11); label.SetTextSize(0.045);
307  if(isData) label.DrawLatex(0.13, 0.93, "#font[61]{CMS} #scale[0.8]{#font[52]{Preliminary}}");
308  else label.DrawLatex(0.13, 0.93, "#font[61]{CMS} #scale[0.8]{#font[52]{Simulation}}");
309  // Drawing luminosity
310  label.SetTextAlign(31);
311  if(isData) {
312  if(do_rpv) label.DrawLatex(0.85, 0.93, "2.69 fb^{-1} (13 TeV)");
313  else label.DrawLatex(0.85, 0.93, "116 pb^{-1} (13 TeV)");
314  }
315  else label.DrawLatex(0.85, 0.93, "Spring15 t#bar{t}");
316 
317  can.SaveAs(pname);
318 
319  for(unsigned his(0); his<2; his++)
320  histo[his]->Delete();
321  fitCurve->Delete();
322 }
323 
324 TString Efficiency(TChain *data, TString den, TString num, float &effic, float &errup, float &errdown){
325  TH1D* histo[2];
326  TString hname, totCut, pname;
327  // den = "("+den+")&&json&&pass";
328  hname = "eden"; totCut = den;
329  histo[0] = new TH1D(hname, "", 1, 0, 1);
330  float denom(data->GetEntries(totCut));
331  histo[0]->SetBinContent(1,denom);
332 
333  hname = "enum"; totCut = "("+den+")&&("+num+")";
334  histo[1] = new TH1D(hname, "", 1, 0, 1);
335  float numer(data->GetEntries(totCut));
336  histo[1]->SetBinContent(1,numer);
337 
338  TGraphAsymmErrors heff(histo[1], histo[0], "cp");
339  //TEfficiency heff(*histo[1], *histo[0]);
340 
341  effic = numer/denom;
342  errup = heff.GetErrorYhigh(0); errdown = heff.GetErrorYlow(0);
343  //float errup(heff.GetEfficiencyErrorUp(0)), errdown(heff.GetEfficiencyErrorLow(0));
344 
345  den.ReplaceAll("&&json","");
346  if(denom) cout<<endl<<"Eff = "<<RoundNumber(numer*100,1,denom)<<"+"<<RoundNumber(errup*100,1)
347  <<"-"<<RoundNumber(errdown*100,1)<<" for num "<<num<<" and "<<den<<" with "<<denom<<" entries"<<endl;
348  else cout<<"Denominator is zero"<<endl;
349 
350  TString efficiency(RoundNumber(numer*100,1,denom)+"^{+"+RoundNumber(errup*100,1)+
351  "}_{-"+RoundNumber(errdown*100,1)+"}");
352 
353  for(unsigned his(0); his<2; his++)
354  histo[his]->Delete();
355 
356  return efficiency;
357 
358 }
359 
styles::setDefaultStyle
void setDefaultStyle()
Definition: styles.cpp:36
styles::LabelSize
float LabelSize
Definition: styles.hpp:35
std
STL namespace.
anonymous_namespace{plot_trigger.cxx}::plot_type
TString plot_type
Definition: plot_trigger.cxx:33
anonymous_namespace{plot_trigger.cxx}::do_rpv
bool do_rpv
Definition: plot_trigger.cxx:32
errorFun
Double_t errorFun(Double_t *x, Double_t *par)
Definition: plot_trigger.cxx:38
format_tag
TString format_tag(TString tag)
Definition: utilities_macros.cpp:870
anonymous_namespace{plot_trigger.cxx}::do_rund
bool do_rund
Definition: plot_trigger.cxx:30
styles::yTitleOffset
float yTitleOffset
Definition: styles.hpp:35
RoundNumber
TString RoundNumber(double num, int decimals, double denom=1.)
Definition: utilities.cpp:191
main
int main()
Definition: plot_trigger.cxx:47
styles::nFont
int nFont
Definition: styles.hpp:33
Efficiency
TString Efficiency(TChain *data, TString den, TString num, float &effic, float &errup, float &errdown)
Definition: plot_trigger.cxx:324
PlotTurnOn
void PlotTurnOn(TChain *data, TString var, int nbins, double minx, double maxx, TString xtitle, TString den, TString num, TString title="", TString ytitle="", float minfit=-1., bool isData=true)
Definition: plot_trigger.cxx:185
parse_card.data
tuple data
Definition: parse_card.py:260
efficiency_map.num
tuple num
Definition: efficiency_map.py:180
utilities_macros.hpp
anonymous_namespace{plot_trigger.cxx}::do_dps
bool do_dps
Definition: plot_trigger.cxx:31
anonymous_namespace{make_rpv_datacard.cxx}::nbins
unsigned int nbins
Definition: make_rpv_datacard.cxx:21
Generated by   doxygen 1.8.11