CMS 3D CMS Logo

BTagDifferentialPlot.cc
Go to the documentation of this file.
5 #include "TF1.h"
6 
10 
12 
13 #include <algorithm>
14 #include <iostream>
15 #include <sstream>
16 using namespace std;
17 
18 
19 
20 BTagDifferentialPlot::BTagDifferentialPlot(double bEff, const ConstVarType& constVariable,
21  const std::string & tagName, unsigned int mc) :
22  fixedBEfficiency (bEff),
23  noProcessing (false), processed(false), constVar(constVariable),
24  constVariableName ("") , diffVariableName ("") ,
25  constVariableValue (999.9, 999.9), commonName("MisidForBEff_" + tagName+"_"),
26  theDifferentialHistoB_d (nullptr),
27  theDifferentialHistoB_u (nullptr),
28  theDifferentialHistoB_s (nullptr),
29  theDifferentialHistoB_c (nullptr),
30  theDifferentialHistoB_b (nullptr),
31  theDifferentialHistoB_g (nullptr),
32  theDifferentialHistoB_ni (nullptr),
33  theDifferentialHistoB_dus(nullptr),
34  theDifferentialHistoB_dusg(nullptr),
35  theDifferentialHistoB_pu (nullptr),
36  mcPlots_(mc)
37 {}
38 
39 
41 }
42 
43 
44 
45 
46 
47 
48 void BTagDifferentialPlot::plot(TCanvas & thePlotCanvas) {
49 
50 // thePlotCanvas = new TCanvas( commonName,
51 // commonName,
52 // btppXCanvas, btppYCanvas);
53 //
54 // if (!btppTitle) gStyle->SetOptTitle(0);
55 
56  if (!processed) return;
57 //fixme:
58  bool btppNI = false;
59  bool btppColour = true;
60 
61  thePlotCanvas.SetFillColor(0);
62  thePlotCanvas.cd(1);
63  gPad->SetLogy(1);
64  gPad->SetGridx(1);
65  gPad->SetGridy(1);
66 
67 // int col_b;
68  int col_c;
69  int col_g;
70  int col_dus;
71  int col_ni;
72 
73 // int mStyle_b;
74  int mStyle_c;
75  int mStyle_g;
76  int mStyle_dus;
77  int mStyle_ni;
78 
79  // marker size(same for all)
80  float mSize = 1.5;
81 
82  if (btppColour) {
83 // col_b = 2;
84  col_c = 6;
85  col_g = 3;
86  col_dus = 4;
87  col_ni = 5;
88 // mStyle_b = 20;
89  mStyle_c = 20;
90  mStyle_g = 20;
91  mStyle_dus = 20;
92  mStyle_ni = 20;
93  }
94  else {
95 // col_b = 1;
96  col_c = 1;
97  col_g = 1;
98  col_dus = 1;
99  col_ni = 1;
100 // mStyle_b = 12;
101  mStyle_c = 22;
102  mStyle_g = 29;
103  mStyle_dus = 20;
104  mStyle_ni = 27;
105  }
106 
107  // for the moment: plot b(to see what the constant b-efficiency is), c, g, uds
108  // b in red
109  // No, do not plot b(because only visible for the soft leptons)
110  // theDifferentialHistoB_b -> GetXaxis()->SetTitle(diffVariableName);
111  // theDifferentialHistoB_b -> GetYaxis()->SetTitle("non b-jet efficiency");
112  // theDifferentialHistoB_b -> GetYaxis()->SetTitleOffset(1.25);
113  // theDifferentialHistoB_b -> SetMaximum(0.4);
114  // theDifferentialHistoB_b -> SetMinimum(1.e-4);
115  // theDifferentialHistoB_b -> SetMarkerColor(col_b);
116  // theDifferentialHistoB_b -> SetLineColor (col_b);
117  // theDifferentialHistoB_b -> SetMarkerSize(mSize);
118  // theDifferentialHistoB_b -> SetMarkerStyle(mStyle_b);
119  // theDifferentialHistoB_b -> SetStats(false);
120  // theDifferentialHistoB_b -> Draw("pe");
121  // c in magenta
122  theDifferentialHistoB_c ->getTH1F() -> SetMaximum(0.4);
123  theDifferentialHistoB_c ->getTH1F() -> SetMinimum(1.e-4);
124  theDifferentialHistoB_c ->getTH1F() -> SetMarkerColor(col_c);
125  theDifferentialHistoB_c ->getTH1F() -> SetLineColor (col_c);
126  theDifferentialHistoB_c ->getTH1F() -> SetMarkerSize(mSize);
127  theDifferentialHistoB_c ->getTH1F() -> SetMarkerStyle(mStyle_c);
128  theDifferentialHistoB_c ->getTH1F() -> SetStats (false);
129  // theDifferentialHistoB_c -> Draw("peSame");
130  theDifferentialHistoB_c ->getTH1F()-> Draw("pe");
131  if (mcPlots_>2) {
132  // uds in blue
133  theDifferentialHistoB_dus ->getTH1F()-> SetMarkerColor(col_dus);
134  theDifferentialHistoB_dus ->getTH1F()-> SetLineColor (col_dus);
135  theDifferentialHistoB_dus ->getTH1F()-> SetMarkerSize(mSize);
136  theDifferentialHistoB_dus ->getTH1F()-> SetMarkerStyle(mStyle_dus);
137  theDifferentialHistoB_dus ->getTH1F()-> SetStats (false);
138  theDifferentialHistoB_dus ->getTH1F()-> Draw("peSame");
139  // g in green
140  // only uds not to confuse
141  theDifferentialHistoB_g ->getTH1F()-> SetMarkerColor(col_g);
142  theDifferentialHistoB_g ->getTH1F()-> SetLineColor (col_g);
143  theDifferentialHistoB_g ->getTH1F()-> SetMarkerSize(mSize);
144  theDifferentialHistoB_g ->getTH1F()-> SetMarkerStyle(mStyle_g);
145  theDifferentialHistoB_g ->getTH1F()-> SetStats (false);
146  theDifferentialHistoB_g ->getTH1F()-> Draw("peSame");
147  }
148  // NI if wanted
149  if (btppNI) {
150  theDifferentialHistoB_ni ->getTH1F()-> SetMarkerColor(col_ni);
151  theDifferentialHistoB_ni ->getTH1F()-> SetLineColor (col_ni);
152  theDifferentialHistoB_ni ->getTH1F()-> SetMarkerSize(mSize);
153  theDifferentialHistoB_ni ->getTH1F()-> SetMarkerStyle(mStyle_ni);
154  theDifferentialHistoB_ni ->getTH1F()-> SetStats (false);
155  theDifferentialHistoB_ni ->getTH1F()-> Draw("peSame");
156  }
157 }
158 
160 {
161  plot(name, ".eps");
162 }
163 
165 {
166  plot(name, ".ps");
167 }
168 
170 {
171  if (!processed) return;
172  TCanvas tc(commonName.c_str(), commonName.c_str());
173  plot(tc);
174  tc.Print((name + commonName + ext).c_str());
175 }
176 
177 
179  setVariableName(); // also sets noProcessing if not OK
180  if (noProcessing) return;
181  bookHisto(ibook);
182  fillHisto();
183  processed = true;
184 }
185 
186 
188 {
189  if (constVar==constETA) {
190  constVariableName = "eta";
191  diffVariableName = "pt";
192  constVariableValue = make_pair(theBinPlotters[0]->etaPtBin().getEtaMin(), theBinPlotters[0]->etaPtBin().getEtaMax());
193  }
194  if (constVar==constPT ) {
195  constVariableName = "pt";
196  diffVariableName = "eta";
197  constVariableValue = make_pair(theBinPlotters[0]->etaPtBin().getPtMin(), theBinPlotters[0]->etaPtBin().getPtMax());
198  }
199 }
200 
201 
202 
204 
205  // vector with ranges
206  vector<float> variableRanges;
207 
208  for (vector<std::shared_ptr<JetTagPlotter>>::const_iterator iP = theBinPlotters.begin();
209  iP != theBinPlotters.end(); ++iP) {
210  const EtaPtBin & currentBin =(*iP)->etaPtBin();
211  if (diffVariableName == "eta") {
212  // only active bins in the variable on x-axis
213  if (currentBin.getEtaActive()) {
214  variableRanges.push_back(currentBin.getEtaMin());
215  // also max if last one
216  if (iP == --theBinPlotters.end()) variableRanges.push_back(currentBin.getEtaMax());
217  }
218  }
219  if (diffVariableName == "pt") {
220  // only active bins in the variable on x-axis
221  if (currentBin.getPtActive()) {
222  variableRanges.push_back(currentBin.getPtMin());
223  // also max if last one
224  if (iP == --theBinPlotters.end()) variableRanges.push_back(currentBin.getPtMax());
225  }
226  }
227  }
228 
229  // to book histo with variable binning -> put into array
230  int nBins = variableRanges.size() - 1;
231  float * binArray = &variableRanges[0];
232 
233  // part of the name common to all flavours
234  std::stringstream stream("");
235  stream << fixedBEfficiency << "_Const_" << constVariableName << "_" << constVariableValue.first << "-";
236  stream << constVariableValue.second << "_" << "_Vs_" << diffVariableName;
237  commonName += stream.str();
238  std::remove(commonName.begin(), commonName.end(), ' ');
239  std::replace(commonName.begin(), commonName.end(), '.', 'v');
240 
242  HistoProviderDQM prov("Btag",label,ibook);
243 
244  theDifferentialHistoB_b =(prov.book1D("B_" + commonName, "B_" + commonName, nBins, binArray));
245  theDifferentialHistoB_c =(prov.book1D("C_" + commonName, "C_" + commonName, nBins, binArray));
246  theDifferentialHistoB_dusg =(prov.book1D("DUSG_" + commonName, "DUSG_" + commonName, nBins, binArray));
247  theDifferentialHistoB_ni =(prov.book1D("NI_" + commonName, "NI_" + commonName, nBins, binArray));
248  theDifferentialHistoB_pu =(prov.book1D("PU_" + commonName, "PU_" + commonName, nBins, binArray));
249 
250  if (mcPlots_>2) {
251  theDifferentialHistoB_d =(prov.book1D("D_" + commonName, "D_" + commonName, nBins, binArray));
252  theDifferentialHistoB_u =(prov.book1D("U_" + commonName, "U_" + commonName, nBins, binArray));
253  theDifferentialHistoB_s =(prov.book1D("S_" + commonName, "S_" + commonName, nBins, binArray));
254  theDifferentialHistoB_g =(prov.book1D("G_" + commonName, "G_" + commonName, nBins, binArray));
255  theDifferentialHistoB_dus =(prov.book1D("DUS_" + commonName, "DUS_" + commonName, nBins, binArray));
256  }
257 }
258 
259 
261  // loop over bins and find corresponding misid. in the MisIdVs..... histo
262  for (vector<std::shared_ptr<JetTagPlotter>>::const_iterator iP = theBinPlotters.begin();
263  iP != theBinPlotters.end(); ++iP) {
264  const EtaPtBin & currentBin =(*iP)->etaPtBin();
265  EffPurFromHistos & currentEffPurFromHistos =(*iP)->getEffPurFromHistos();
266  //
267  bool isActive = true;
268  double valueXAxis = -999.99;
269  // find right bin based on middle of the interval
270  if (diffVariableName == "eta") {
271  isActive = currentBin.getEtaActive();
272  valueXAxis = 0.5 *(currentBin.getEtaMin() + currentBin.getEtaMax());
273  } else if (diffVariableName == "pt" ) {
274  isActive = currentBin.getPtActive();
275  valueXAxis = 0.5 *(currentBin.getPtMin() + currentBin.getPtMax());
276  } else {
277  throw cms::Exception("Configuration")
278  << "====>>>> BTagDifferentialPlot::fillHisto() : illegal diffVariableName = " << diffVariableName << endl;
279  }
280 
281  // for the moment: ignore inactive bins
282  //(maybe later: if a Bin is inactive -> set value to fill well below left edge of histogram to have it in the underflow)
283 
284  if (!isActive) continue;
285 
286  // to have less lines of code ....
287  vector< pair<TH1F*,TH1F*> > effPurDifferentialPairs;
288 
289  // all flavours(b is a good cross check! must be constant and = fixed b-efficiency)
290  // get histo; find the bin of the fixed b-efficiency in the histo and get misid; fill
291 
292 
293  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_b() , theDifferentialHistoB_b ->getTH1F() ));
294  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_c() , theDifferentialHistoB_c ->getTH1F() ));
295  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_dusg(), theDifferentialHistoB_dusg->getTH1F()));
296  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_ni() , theDifferentialHistoB_ni ->getTH1F() ));
297  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_pu() , theDifferentialHistoB_pu->getTH1F() ));
298  if (mcPlots_>2) {
299  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_d() , theDifferentialHistoB_d ->getTH1F() ));
300  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_u() , theDifferentialHistoB_u ->getTH1F() ));
301  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_s() , theDifferentialHistoB_s ->getTH1F() ));
302  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_g() , theDifferentialHistoB_g ->getTH1F() ));
303  effPurDifferentialPairs.push_back(make_pair(currentEffPurFromHistos.getEffFlavVsBEff_dus(), theDifferentialHistoB_dus->getTH1F() ));
304  }
305 
306  for (vector< pair<TH1F*,TH1F*> >::const_iterator itP = effPurDifferentialPairs.begin();
307  itP != effPurDifferentialPairs.end(); ++itP) {
308  TH1F * effPurHist = itP->first;
309  TH1F * diffHist = itP->second;
310  pair<double, double> mistag = getMistag(fixedBEfficiency, effPurHist);
311  int iBinSet = diffHist->FindBin(valueXAxis);
312  diffHist->SetBinContent(iBinSet, mistag.first);
313  diffHist->SetBinError(iBinSet, mistag.second);
314  }
315  }
316 
317 }
318 
319 pair<double, double>
320 BTagDifferentialPlot::getMistag(const double& fixedBEfficiency, TH1F * effPurHist)
321 {
322  int iBinGet = effPurHist->FindBin(fixedBEfficiency);
323  double effForBEff = effPurHist->GetBinContent(iBinGet);
324  double effForBEffErr = effPurHist->GetBinError (iBinGet);
325 
326  if (effForBEff==0. && effForBEffErr==0.) {
327  // The bin was empty. Could be that it was not filled, as the scan-plot
328  // did not have an entry at the requested value, or that the curve
329  // is above or below.
330  // Fit a plynomial, and evaluate the mistag at the requested value.
331  int fitStatus;
332  //need our own copy for thread safety
333  TF1 myfunc("myfunc","pol4");
334  try {
335  fitStatus = effPurHist->Fit(&myfunc, "q");
336  }catch(cms::Exception& iException) {
337  return pair<double, double>(effForBEff, effForBEffErr);
338  }
339  if (fitStatus != 0) {
340  edm::LogWarning("BTagDifferentialPlot")<<"Fit failed to hisogram " << effPurHist->GetTitle() << ", perhaps because too few entries = " << effPurHist->GetEntries() <<". This bin will be missing in plots at fixed b efficiency.";
341  // } else {
342  // edm::LogInfo("BTagDifferentialPlot")<<"Fit OK to hisogram " << effPurHist->GetTitle() << " entries = " << effPurHist->GetEntries();
343  return pair<double, double>(effForBEff, effForBEffErr);
344  }
345  effForBEff = myfunc.Eval(fixedBEfficiency);
346  effPurHist->RecursiveRemove(&myfunc);
347  //Error: first non-empty bin on the right and take its error
348  for (int i = iBinGet+1; i< effPurHist->GetNbinsX(); ++i) {
349  if (effPurHist->GetBinContent(i)!=0) {
350  effForBEffErr = effPurHist->GetBinError(i);
351  break;
352  }
353  }
354  }
355 
356  return pair<double, double>(effForBEff, effForBEffErr);
357 }
358 
MonitorElement * theDifferentialHistoB_b
MonitorElement * theDifferentialHistoB_g
TH1F * getEffFlavVsBEff_g()
BTagDifferentialPlot(double bEff, const ConstVarType &constVariable, const std::string &tagName, unsigned int mc)
TH1F * getTH1F() const
MonitorElement * theDifferentialHistoB_s
TH1F * getEffFlavVsBEff_dusg()
TH1F * getEffFlavVsBEff_c()
#define nullptr
MonitorElement * theDifferentialHistoB_pu
MonitorElement * theDifferentialHistoB_dusg
def replace(string, replacements)
double getEtaMax() const
Definition: EtaPtBin.h:37
double getEtaMin() const
Definition: EtaPtBin.h:36
std::vector< std::shared_ptr< JetTagPlotter > > theBinPlotters
void process(DQMStore::IBooker &ibook)
double getPtMin() const
Definition: EtaPtBin.h:40
void psPlot(const std::string &name)
bool getEtaActive() const
Get rapidity/pt ranges and check whether rapidity/pt cuts are active.
Definition: EtaPtBin.h:35
MonitorElement * theDifferentialHistoB_dus
char const * label
virtual MonitorElement * book1D(const std::string &name, const std::string &title, const int &nchX, const double &lowX, const double &highX)
bool getPtActive() const
Definition: EtaPtBin.h:39
TH1F * getEffFlavVsBEff_u()
void bookHisto(DQMStore::IBooker &ibook)
std::pair< double, double > getMistag(const double &fixedBEfficiency, TH1F *effPurHist)
std::pair< double, double > constVariableValue
double getPtMax() const
Definition: EtaPtBin.h:41
TH1F * getEffFlavVsBEff_dus()
TH1F * getEffFlavVsBEff_d()
MonitorElement * theDifferentialHistoB_u
def remove(d, key, TELL=False)
Definition: MatrixUtil.py:212
TH1F * getEffFlavVsBEff_ni()
MonitorElement * theDifferentialHistoB_d
TH1F * getEffFlavVsBEff_pu()
TH1F * getEffFlavVsBEff_b()
Definition: memstream.h:15
void plot(TCanvas &theCanvas)
TH1F * getEffFlavVsBEff_s()
void epsPlot(const std::string &name)
MonitorElement * theDifferentialHistoB_c
MonitorElement * theDifferentialHistoB_ni