CMS 3D CMS Logo

JetTagPlotter.cc
Go to the documentation of this file.
5 
6 #include <iostream>
7 
8 using namespace std;
9 using namespace RecoBTag;
10 
11 
13  const edm::ParameterSet& pSet, unsigned int mc,
14  bool wf, DQMStore::IBooker & ibook, bool doCTagPlots/*=false*/,
15  bool doDifferentialPlots/*=false*/, double discrCut/*=-999.*/) :
16  BaseBTagPlotter(tagName, etaPtBin),
17  discrStart_(pSet.getParameter<double>("discriminatorStart")),
18  discrEnd_(pSet.getParameter<double>("discriminatorEnd")),
19  nBinEffPur_(pSet.getParameter<int>("nBinEffPur")),
20  startEffPur_(pSet.getParameter<double>("startEffPur")),
21  endEffPur_(pSet.getParameter<double>("endEffPur")),
22  mcPlots_(mc), willFinalize_(wf), doCTagPlots_(doCTagPlots),
23  doDifferentialPlots_(doDifferentialPlots),
24  cutValue_(discrCut) {
25 
26  // to have a shorter name .....
27  const std::string & es = theExtensionString;
28  const std::string jetTagDir(es.substr(1));
29 
30  if (willFinalize_) return;
31 
32  //added to count the number of jets by event : 0=DATA or NI, 1to5=quarks u,d,s,c,b , 6=gluon
33  int nFl = 1;
34  if (mcPlots_) nFl = 8;
35  nJets_.resize(nFl, 0);
36 
37  if (mcPlots_) {
38  // jet flavour
39  dJetFlav_ = std::make_unique<FlavourHistograms<int>>
40  ("jetFlavour" + es, "Jet Flavour", 22, -0.5, 21.5,
41  false, false, false, "b", jetTagDir, mcPlots_, ibook);
42  }
43 
44  // jet multiplicity
45  jetMultiplicity_ = std::make_unique<FlavourHistograms<int>>
46  ("jetMultiplicity" + es, "Jet Multiplicity", 11, -0.5, 10.5,
47  false, true, true, "b", jetTagDir, mcPlots_, ibook);
48 
49  // Discriminator: again with reasonable binning
50  dDiscriminator_ = std::make_unique<FlavourHistograms<double>>
51  ("discr" + es, "Discriminator", 102, discrStart_, discrEnd_,
52  false, true, true, "b", jetTagDir, mcPlots_, ibook);
53  dDiscriminator_->settitle("Discriminant");
54  // reconstructed jet momentum
55  dJetRecMomentum_ = std::make_unique<FlavourHistograms<double>>
56  ("jetMomentum" + es, "jet momentum", 350, 0.0, 350.0,
57  false, false, true, "b", jetTagDir, mcPlots_, ibook);
58 
59  // reconstructed jet transverse momentum
60  dJetRecPt_ = std::make_unique<FlavourHistograms<double>>
61  ("jetPt" + es, "jet pt", 350, 0.0, 350.0,
62  false, false, true, "b", jetTagDir, mcPlots_, ibook);
63 
64  // reconstructed jet eta
65  dJetRecPseudoRapidity_ = std::make_unique<FlavourHistograms<double>>
66  ("jetEta" + es, "jet eta", 20, -etaPtBin.getEtaMax(), etaPtBin.getEtaMax(),
67  false, false, true, "b", jetTagDir, mcPlots_, ibook);
68 
69  // reconstructed jet phi
70  dJetRecPhi_ = std::make_unique<FlavourHistograms<double>>
71  ("jetPhi" + es, "jet phi", 20, -M_PI, M_PI,
72  false, false, true, "b", jetTagDir, mcPlots_, ibook);
73 
75  // jet Phi larger than requested discrimnator cut
76  dJetPhiDiscrCut_ = std::make_unique<FlavourHistograms<double>>("jetPhi_diffEff" + es, "Efficiency vs. jet Phi for discriminator above cut",
77  20, -M_PI, M_PI, false, false, true, "b", jetTagDir, mcPlots_, ibook);
78 
79  // jet Eta larger than requested discrimnator cut
80  dJetPseudoRapidityDiscrCut_ = std::make_unique<FlavourHistograms<double>>("jetEta_diffEff" + es, "Efficiency vs. jet eta for discriminator above cut",
81  20, -etaPtBin.getEtaMax(), etaPtBin.getEtaMax(), false, false, true, "b", jetTagDir, mcPlots_, ibook);
82  }
83 }
84 
85 
87 
89 {
90  if (!willFinalize_) {
91  dJetFlav_->epsPlot(name);
92  jetMultiplicity_->epsPlot(name);
93  dDiscriminator_->epsPlot(name);
94  dJetRecMomentum_->epsPlot(name);
95  dJetRecPt_->epsPlot(name);
96  dJetRecPseudoRapidity_->epsPlot(name);
97  dJetRecPhi_->epsPlot(name);
98  }
99  else {
100  effPurFromHistos_->epsPlot(name);
101  }
102 }
103 
105 {
106  std::string cName = "JetTagPlots"+ theExtensionString;
107  setTDRStyle()->cd();
108  TCanvas canvas(cName.c_str(), cName.c_str(), 600, 900);
109  canvas.UseCurrentStyle();
110 
111  canvas.Divide(2,3);
112  canvas.Print((name + cName + ".ps[").c_str());
113  if (!willFinalize_) {
114  canvas.cd(1);
115  dJetFlav_->plot();
116  canvas.cd(2);
117  canvas.cd(3);
118  dDiscriminator_->plot();
119  canvas.cd(4);
120  dJetRecMomentum_->plot();
121  canvas.cd(5);
122  dJetRecPt_->plot();
123  canvas.cd(6);
124  dJetRecPseudoRapidity_->plot();
125  canvas.Print((name + cName + ".ps").c_str());
126  canvas.Clear();
127  canvas.Divide(2,3);
128 
129  jetMultiplicity_->plot();
130  canvas.Print((name + cName + ".ps").c_str());
131  canvas.Clear();
132 
133  canvas.cd(1);
134  dJetRecPhi_->plot();
135  canvas.cd(2);
136  canvas.cd(3);
137  canvas.cd(4);
138  } else {
139  canvas.cd(5);
140  effPurFromHistos_->discriminatorNoCutEffic().plot();
141  canvas.cd(6);
142  effPurFromHistos_->discriminatorCutEfficScan().plot();
143  canvas.Print((name + cName + ".ps").c_str());
144  canvas.Clear();
145  canvas.Divide(2,3);
146  canvas.cd(1);
147  effPurFromHistos_->plot();
148  }
149  canvas.Print((name + cName + ".ps").c_str());
150  canvas.Print((name + cName + ".ps]").c_str());
151 }
152 
153 void JetTagPlotter::analyzeTag() //here jetFlavour not needed
154 {
155  //to use on data
156  jetMultiplicity_->fill(-1, nJets_[0]);
157  nJets_[0] = 0; //reset to 0 before the next event
158 }
159 
161 {
162  if (mcPlots_) {
163  //to use with MC
164  int totNJets = 0;
165  int udsNJets = 0;
166  int udsgNJets = 0;
167  for (int i = 0; i < 8; i++) {
168  totNJets += nJets_[i];
169  if(i > 0 && i < 4) udsNJets += nJets_[i];
170  if((i > 0 && i < 4) || i == 6) udsgNJets += nJets_[i];
171  if(i <= 5 && i >= 1) jetMultiplicity_->fill(i, nJets_[i], w);
172  else if (i==6) jetMultiplicity_->fill(21, nJets_[i], w);
173  else if (i==7) jetMultiplicity_->fill(20, nJets_[i], w);
174  else jetMultiplicity_->fill(0, nJets_[i], w);
175  nJets_[i] = 0; //reset to 0 before the next event
176  }
177  jetMultiplicity_->fill(-1, totNJets, w); //total number of jets in the event
178  jetMultiplicity_->fill(123, udsNJets, w);
179  jetMultiplicity_->fill(12321, udsgNJets, w);
180  } else {
181  int totNJets = 0;
182  for (int i = 0; i < 8; i++) {
183  totNJets += nJets_[i];
184  nJets_[i] = 0;
185  }
186  jetMultiplicity_->fill(-1, totNJets, w);
187  }
188 }
189 
191  double jec,
192  float discriminator,
193  int jetFlavour,
194  float w/*=1*/)
195 {
196  if (mcPlots_) {
197  dJetFlav_->fill(jetFlavour, jetFlavour, w);
198  if (abs(jetFlavour) > 0 && abs(jetFlavour) < 6) nJets_[abs(jetFlavour)] += 1; //quarks 1 to 5
199  else if (abs(jetFlavour) == 21) nJets_[6] += 1; //gluons
200  else if (jetFlavour == 20) nJets_[7] += 1; //PU
201  else nJets_[0] += 1; //NI
202  } else {
203  nJets_[0] += 1;
204  }
205  if (edm::isNotFinite(discriminator)) dDiscriminator_->fill(jetFlavour, -999.0, w);
206  else dDiscriminator_->fill(jetFlavour, discriminator, w);
207  dJetRecMomentum_->fill(jetFlavour, jet.p() * jec, w);
208  dJetRecPt_->fill(jetFlavour, jet.pt() * jec, w);
209  dJetRecPseudoRapidity_->fill(jetFlavour, jet.eta(), w);
210  dJetRecPhi_->fill(jetFlavour, jet.phi(), w);
211  if (doDifferentialPlots_) {
212  if (edm::isFinite(discriminator) && discriminator > cutValue_) {
213  dJetPhiDiscrCut_->fill(jetFlavour, jet.phi(), w);
214  dJetPseudoRapidityDiscrCut_->fill(jetFlavour, jet.eta(), w);
215  }
216  }
217 }
218 
219 
221  double jec,
222  int jetFlavour,
223  float w/*=1*/)
224 {
225  if (mcPlots_) {
226  dJetFlav_->fill(jetFlavour, jetFlavour, w);
227  if (abs(jetFlavour) > 0 && abs(jetFlavour) < 6) nJets_[abs(jetFlavour)] += 1; //quarks 1 to 5
228  else if (abs(jetFlavour) == 21) nJets_[6] += 1; //gluons
229  else if (jetFlavour == 20) nJets_[7] += 1; //PU
230  else nJets_[0] += 1; //NI
231  } else {
232  nJets_[0] += 1;
233  }
234  const auto& discriminator = jetTag.second;
235  if (edm::isNotFinite(discriminator)) dDiscriminator_->fill(jetFlavour, -999.0, w);
236  else dDiscriminator_->fill(jetFlavour, discriminator, w);
237  dJetRecMomentum_->fill(jetFlavour, jetTag.first->p() * jec, w);
238  dJetRecPt_->fill(jetFlavour, jetTag.first->pt() * jec, w);
239  dJetRecPseudoRapidity_->fill(jetFlavour, jetTag.first->eta(), w);
240  dJetRecPhi_->fill(jetFlavour, jetTag.first->phi(), w);
241  if (doDifferentialPlots_) {
243  dJetPhiDiscrCut_->fill(jetFlavour, jetTag.first->phi(), w);
244  dJetPseudoRapidityDiscrCut_->fill(jetFlavour, jetTag.first->eta(), w);
245  }
246  }
247 }
248 
250 {
251  //
252  // final processing:
253  // produce the misid. vs. eff histograms
254  //
255  const std::string & es = theExtensionString;
256  const std::string jetTagDir(es.substr(1));
257  dDiscriminator_ = std::make_unique<FlavourHistograms<double>>("discr" + es, "Discriminator", 102, discrStart_, discrEnd_, "b", jetTagDir, mcPlots_, igetter_);
258 
259  effPurFromHistos_ = std::make_unique<EffPurFromHistos>(*dDiscriminator_, jetTagDir, mcPlots_, ibook_, nBinEffPur_, startEffPur_, endEffPur_);
260  effPurFromHistos_->doCTagPlots(doCTagPlots_);
261  effPurFromHistos_->compute(ibook_);
262 
263  // Produce the differentiel efficiency vs. kinematical variables
264  if (doDifferentialPlots_) {
265  dJetRecPhi_ = std::make_unique<FlavourHistograms<double>>("jetPhi" + es, "jet phi", 20, -M_PI, M_PI, "b", jetTagDir, mcPlots_, igetter_);
266  dJetPhiDiscrCut_ = std::make_unique<FlavourHistograms<double>>("jetPhi_diffEff" + es, "Efficiency vs. jet Phi for discriminator above cut", 20, -M_PI, M_PI, "b", jetTagDir, mcPlots_, igetter_);
267  dJetPhiDiscrCut_->divide(*dJetRecPhi_);
268  dJetPhiDiscrCut_->setEfficiencyFlag();
269 
270  dJetRecPseudoRapidity_ = std::make_unique<FlavourHistograms<double>>("jetEta" + es, "jet eta", 20, -etaPtBin_.getEtaMax(), etaPtBin_.getEtaMax(), "b", jetTagDir, mcPlots_, igetter_);
271  dJetPseudoRapidityDiscrCut_ = std::make_unique<FlavourHistograms<double>>("jetEta_diffEff" + es, "Efficiency vs. jet eta for discriminator above cut", 20, -etaPtBin_.getEtaMax(), etaPtBin_.getEtaMax(), "b", jetTagDir, mcPlots_, igetter_);
272  dJetPseudoRapidityDiscrCut_->divide(*dJetRecPseudoRapidity_);
273  dJetPseudoRapidityDiscrCut_->setEfficiencyFlag();
274  }
275 }
276 
std::unique_ptr< FlavourHistograms< double > > dJetPseudoRapidityDiscrCut_
Definition: JetTagPlotter.h:92
double eta() const final
momentum pseudorapidity
std::unique_ptr< FlavourHistograms< double > > dJetPhiDiscrCut_
Definition: JetTagPlotter.h:89
const double w
Definition: UKUtility.cc:23
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
std::unique_ptr< FlavourHistograms< int > > dJetFlav_
Definition: JetTagPlotter.h:71
const std::string theExtensionString
void finalize(DQMStore::IBooker &ibook_, DQMStore::IGetter &igetter_) override
Base class for all types of Jets.
Definition: Jet.h:20
double endEffPur_
Definition: JetTagPlotter.h:52
JetFloatAssociation::value_type JetTag
Definition: JetTag.h:17
void psPlot(const std::string &name) override
bool doDifferentialPlots_
Definition: JetTagPlotter.h:60
std::unique_ptr< FlavourHistograms< int > > jetMultiplicity_
Definition: JetTagPlotter.h:66
double getEtaMax() const
Definition: EtaPtBin.h:37
double pt() const final
transverse momentum
double discrStart_
Definition: JetTagPlotter.h:48
JetTagPlotter(const std::string &tagName, const EtaPtBin &etaPtBin, const edm::ParameterSet &pSet, unsigned int mc, bool willFinalize, DQMStore::IBooker &ibook, bool doCTagPlots=false, bool doDifferentialPlots=false, double discrCut=-999.)
unsigned int mcPlots_
Definition: JetTagPlotter.h:54
constexpr bool isFinite(T x)
~JetTagPlotter() override
def setTDRStyle()
Definition: plotscripts.py:89
std::unique_ptr< FlavourHistograms< double > > dJetRecPseudoRapidity_
Definition: JetTagPlotter.h:83
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::unique_ptr< FlavourHistograms< double > > dJetRecMomentum_
Definition: JetTagPlotter.h:77
double startEffPur_
Definition: JetTagPlotter.h:51
double cutValue_
Definition: JetTagPlotter.h:61
#define M_PI
const EtaPtBin etaPtBin_
double p() const final
magnitude of momentum vector
Definition: Tools.h:23
std::unique_ptr< FlavourHistograms< double > > dDiscriminator_
Definition: JetTagPlotter.h:74
def canvas(sub, attr)
Definition: svgfig.py:482
void epsPlot(const std::string &name) override
double discrEnd_
Definition: JetTagPlotter.h:49
std::unique_ptr< FlavourHistograms< double > > dJetRecPhi_
Definition: JetTagPlotter.h:86
double phi() const final
momentum azimuthal angle
std::unique_ptr< FlavourHistograms< double > > dJetRecPt_
Definition: JetTagPlotter.h:80
std::vector< int > nJets_
Definition: JetTagPlotter.h:63
std::unique_ptr< EffPurFromHistos > effPurFromHistos_
Definition: JetTagPlotter.h:69