001 #include "JetTagTools/NewLikelihoodTool.h"
002
003 #include "JetTagTools/HistoHelperRoot.h"
004 #include "JetTagCalibration/CalibrationBroker.h"
005
006 #include <cmath>
007 #include <string>
008 #include <utility>
009 #include <vector>
010
011 #include "TH1.h"
012 #include "TH2.h"
013 #include "TH3.h"
014
015 namespace Analysis {
016
017 NewLikelihoodTool::NewLikelihoodTool(const std::string& t, const std::string& n, const IInterface* p) :
018 AthAlgTool(t,n,p),
019 m_taggerName("undefined"),
020 m_hypotheses(std::vector<std::string>()),
021 m_calibrationTool("BTagCalibrationBroker"),
022 m_normalizedProb(true),
023 m_interpolate(false),
024 m_smoothNTimes(1),
025 m_vetoSmoothingOf(std::vector<std::string>()),
026 m_lhVariableValues(std::vector<Slice> ()),
027 m_likelihoodVector(std::vector<double>()),
028 m_histograms(std::vector<std::string>()) {
029
030 declareInterface<NewLikelihoodTool>(this);
031 declareProperty("taggerName", m_taggerName);
032 declareProperty("hypotheses", m_hypotheses);
033 declareProperty("calibrationTool", m_calibrationTool);
034 declareProperty("normalizedProb", m_normalizedProb);
035 declareProperty("interpolate", m_interpolate);
036 declareProperty("smoothNTimes", m_smoothNTimes);
037 declareProperty("vetoSmoothingOf", m_vetoSmoothingOf);
038 m_hypotheses.push_back("B");
039 m_hypotheses.push_back("U");
040 m_vetoSmoothingOf.push_back("/N2T");
041 m_vetoSmoothingOf.push_back("/Sip3D");
042 }
043
044 NewLikelihoodTool::~NewLikelihoodTool() {
045 }
046
047 StatusCode NewLikelihoodTool::initialize() {
048
049 if ( m_calibrationTool.retrieve().isFailure() ) {
050 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_calibrationTool);
051 return StatusCode::FAILURE;
052 } else {
053 ATH_MSG_DEBUG("#BTAG# Retrieved tool " << m_calibrationTool);
054 }
055 return StatusCode::SUCCESS;
056 }
057
058 StatusCode NewLikelihoodTool::finalize() {
059 return StatusCode::SUCCESS;
060 }
061
062 void NewLikelihoodTool::setLhVariableValue(std::vector<Slice>& value) {
063 m_lhVariableValues=value;
064 }
065
066 void NewLikelihoodTool::setLhVariableValue(Slice& value) {
067 std::vector<Slice> tmpVector;
068 tmpVector.push_back(value);
069 setLhVariableValue(tmpVector);
070 }
071
072 void NewLikelihoodTool::defineHypotheses(const std::vector<std::string>& hyp) {
073 m_hypotheses = hyp;
074 }
075
076 void NewLikelihoodTool::defineHistogram(const std::string& hname) {
077 m_histograms.push_back(hname);
078 std::vector<std::string> gradeList = this->gradeList(hname);
079 m_calibrationTool->registerHistogram(m_taggerName, hname);
080 }
081
082 void NewLikelihoodTool::printStatus() const {
083 msg(MSG::INFO) << "#BTAG# - hypotheses : ";
084 for(unsigned int ih=0;ih<m_hypotheses.size();ih++) msg(MSG::INFO) << m_hypotheses[ih] << ", ";
085 msg(MSG::INFO) << endreq;
086 msg(MSG::INFO) << "#BTAG# - histograms : " << endreq;
087 for(unsigned int ih=0;ih<m_histograms.size();ih++) msg(MSG::INFO) << m_histograms[ih] << endreq;
088 msg(MSG::INFO) << "#BTAG# - status of underlying calibrations: " << endreq;
089 m_calibrationTool->printStatus();
090 }
091
092 void NewLikelihoodTool::clear() {
093 m_lhVariableValues.clear();
094 m_likelihoodVector.clear();
095 }
096
097 std::vector<std::string> NewLikelihoodTool::gradeList(const std::string& histoName) {
098 ATH_MSG_VERBOSE("#BTAG# gradeList() called for " << histoName);
099 const std::string delimSlash("/");
100
101 unsigned int nSlash = 0;
102 std::string::size_type slashPos;
103 for(unsigned int i = 0;
104 (slashPos = histoName.find(delimSlash, i)) != std::string::npos;
105 i = slashPos + 1) nSlash++;
106
107 slashPos = histoName.find_first_of(delimSlash);
108 std::string newName = histoName.substr(slashPos+1);
109 slashPos = newName.find_first_of(delimSlash);
110 std::string grades = newName.substr(0,slashPos);
111 std::string hhname = newName.substr(slashPos+1);
112 if(nSlash<2) grades = "";
113 ATH_MSG_VERBOSE("#BTAG# -> grades: " << grades);
114 ATH_MSG_VERBOSE("#BTAG# -> hhname: " << hhname);
115
116 std::vector<std::string> gradeList;
117 if(grades!="") {
118 const std::string delimUds("_");
119 std::string::size_type sPos, sEnd, sLen;
120 sPos = grades.find_first_not_of(delimUds);
121 while ( sPos != std::string::npos ) {
122 sEnd = grades.find_first_of(delimUds, sPos);
123 if(sEnd==std::string::npos) sEnd = grades.length();
124 sLen = sEnd - sPos;
125 std::string word = grades.substr(sPos,sLen);
126 ATH_MSG_DEBUG("#BTAG# --> grade = " << word);
127 gradeList.push_back(word);
128 sPos = grades.find_first_not_of(delimUds, sEnd);
129 }
130 }
131
132 gradeList.push_back(hhname);
133 return gradeList;
134 }
135
136 TH1* NewLikelihoodTool::prepareHistogram(const std::string& hypo, const std::string& hname) {
137 TH1* histoSum(0);
138 bool updated = false;
139 std::string channelName = m_calibrationTool->channelName(hname);
140 std::string histoName = m_calibrationTool->histoName(hname);
141 std::string actualName = hypo + "/" + histoName;
142 std::string longName = channelName + "#" + actualName;
143 ATH_MSG_VERBOSE("#BTAG# preparing histogram " << longName);
144 ATH_MSG_VERBOSE("#BTAG# -> channel is " << channelName);
145 ATH_MSG_VERBOSE("#BTAG# -> histo is " << histoName);
146
147 std::vector<std::string> gradeList = this->gradeList(histoName);
148
149 if(1==gradeList.size()) {
150 ATH_MSG_DEBUG("#BTAG# Histo "<<actualName<<" has no grade: direct retrieval");
151 std::pair<TH1*, bool> rth = m_calibrationTool->retrieveHistogram(m_taggerName,
152 channelName,
153 actualName);
154 histoSum = rth.first;
155 updated = rth.second;
156 if(updated) {
157 ATH_MSG_VERBOSE("#BTAG# Smoothing histogram " << longName << " ...");
158 this->smoothAndNormalizeHistogram(histoSum, longName);
159 m_calibrationTool->updateHistogramStatus(m_taggerName, channelName, actualName, false);
160 }
161 }
162
163 if(2==gradeList.size()) {
164 ATH_MSG_DEBUG("#BTAG# Histo "<<actualName<<" has only one grade: direct retrieval");
165 std::pair<TH1*, bool> rth = m_calibrationTool->retrieveHistogram(m_taggerName,
166 channelName,
167 actualName);
168 histoSum = rth.first;
169 updated = rth.second;
170 if(updated) {
171 ATH_MSG_VERBOSE("#BTAG# Smoothing histogram " << longName << " ...");
172 this->smoothAndNormalizeHistogram(histoSum, longName);
173 m_calibrationTool->updateHistogramStatus(m_taggerName, channelName, actualName, false);
174 }
175 }
176
177 if(gradeList.size()>2) {
178 TH1* histo(0);
179 bool AllUpdated = true;
180 ATH_MSG_DEBUG("Histo " << actualName << " has " << (gradeList.size()-1) << " grades:");
181 for(unsigned int i=0;i<(gradeList.size()-1);i++) {
182 actualName = hypo+"/"+gradeList[i]+gradeList[gradeList.size()-1];
183 ATH_MSG_VERBOSE("#BTAG# -> retrieving histo for grade " << i << " " << gradeList[i] << ": ");
184 std::pair<TH1*, bool> rthh = m_calibrationTool->retrieveHistogram(m_taggerName,
185 channelName,
186 actualName);
187 histo = rthh.first;
188 updated = rthh.second;
189 if(histo) ATH_MSG_VERBOSE("#BTAG# histo " << actualName
190 << " has " << histo->GetEntries() << " entries. Updated: "
191 << updated );
192 if(updated) {
193 m_calibrationTool->updateHistogramStatus(m_taggerName, channelName, actualName, false);
194 }
195 AllUpdated = AllUpdated && updated;
196 if(0==i) {
197 histoSum = histo;
198 } else {
199 if(histo&&histoSum) histoSum->Add(histo,1.);
200 }
201 }
202 ATH_MSG_VERBOSE("#BTAG# Smoothing histogram " << longName << " ...");
203 this->smoothAndNormalizeHistogram(histoSum, longName);
204 }
205 return histoSum;
206 }
207
208 void NewLikelihoodTool::smoothAndNormalizeHistogram(TH1* h, const std::string& hname = "") {
209 MsgStream mlog(msgSvc(),name());
210 if(h) {
211 double norm = h->Integral();
212 if(norm) {
213
214 bool veto = false;
215 for(unsigned int iv=0; iv<m_vetoSmoothingOf.size(); iv++) {
216 if(hname.find(m_vetoSmoothingOf[iv])!=std::string::npos) {
217 veto = true;
218 ATH_MSG_VERBOSE("#BTAG# Smoothing of " << hname << " is vetoed !");
219 break;
220 }
221 }
222 if(1==h->GetDimension() && m_smoothNTimes) {
223 if(!veto) {
224 if(norm>10000)h->Smooth(m_smoothNTimes);
225 else h->Smooth((int)(m_smoothNTimes*100./sqrt(norm)));
226 }
227 }
228 if(2==h->GetDimension()) {
229 int m2d=3;
230 if(hname.find("#B/")==std::string::npos) m2d=5;
231 if(!veto) {
232 HistoHelperRoot::smoothASH2D(dynamic_cast<TH2*>(h), m2d, m2d, msgLvl(MSG::DEBUG));
233 }
234 }
235 if(3==h->GetDimension()) {
236 int m3d1=3;
237 int m3d3=2;
238 if(!veto) {
239 HistoHelperRoot::smoothASH3D(dynamic_cast<TH3*>(h), m3d1, m3d1, m3d3, msgLvl(MSG::DEBUG));
240 }
241 }
242
243 norm = h->Integral();
244 h->Scale(1./norm);
245 } else {
246 ATH_MSG_WARNING("#BTAG# Histo "<<h<<" is empty!");
247 }
248 }
249 }
250
251 std::vector<double> NewLikelihoodTool::calculateLikelihood() {
252 m_likelihoodVector.clear();
253 ATH_MSG_VERBOSE("#BTAG# calculate called for" << m_hypotheses.size() << "hypotheses.");
254 std::vector<double> probDensityPerEventClassAllVariables;
255 probDensityPerEventClassAllVariables.resize(m_hypotheses.size());
256 for (unsigned int i = 0 ; i < m_hypotheses.size(); ++i) {
257 probDensityPerEventClassAllVariables[i]=1.;
258 }
259 ATH_MSG_VERBOSE("#BTAG# -- lhVarVal size= " << m_lhVariableValues.size());
260
261 for (unsigned int iel = 0; iel<m_lhVariableValues.size(); iel++) {
262 msg(MSG::VERBOSE) << "#BTAG# -- element " << iel << " "
263 << m_lhVariableValues[iel].name << endreq;
264 int ncompo = m_lhVariableValues[iel].composites.size();
265 msg(MSG::VERBOSE) << "#BTAG# -- element " << iel << " "
266 << m_lhVariableValues[iel].name
267 << " has " << ncompo << " composites." << endreq;
268
269
270 for (int icompo = 0;icompo<ncompo;icompo++) {
271 double sum(0.);
272 std::vector<double> tmpVector;
273 std::string histName = m_lhVariableValues[iel].composites[icompo].name;
274 int idim = m_lhVariableValues[iel].composites[icompo].atoms.size();
275 msg(MSG::VERBOSE) << "#BTAG# -- composite " << icompo << " histo= "
276 << histName << " dim= " << idim << endreq;
277 for (unsigned int ihyp = 0 ; ihyp < m_hypotheses.size(); ++ihyp) {
278 TH1* tmpHisto = this->prepareHistogram(m_hypotheses[ihyp],histName);
279 if(tmpHisto) {
280 if(1==idim) {
281 double valuex = m_lhVariableValues[iel].composites[icompo].atoms[0].value;
282 int binx = (tmpHisto->GetXaxis())->FindBin(valuex);
283 if(valuex >= tmpHisto->GetXaxis()->GetXmax()) binx = tmpHisto->GetXaxis()->GetNbins();
284 if(valuex <= tmpHisto->GetXaxis()->GetXmin()) binx = 1;
285 double tmp = tmpHisto->GetBinContent(binx);
286 msg(MSG::VERBOSE) << "#BTAG# For hypothesis= " << m_hypotheses[ihyp]
287 << " (1D) actual value= " << valuex
288 << " --> bin= " << binx << " f = " << tmp;
289 if(m_interpolate) {
290 tmp = HistoHelperRoot::Interpol1d(valuex, dynamic_cast<TH1*>(tmpHisto));
291 msg(MSG::VERBOSE) << " interpolated f = " << tmp;
292 }
293 if(m_normalizedProb) {
294 msg(MSG::VERBOSE) << " (normalized)" << endreq;
295 } else {
296 double binw = tmpHisto->GetBinWidth(binx);
297 msg(MSG::VERBOSE) << " binw= " << binw;
298 if(0==binw) {
299 msg(MSG::ERROR) << "bin width is 0" << endreq;
300 } else {
301 tmp /= binw;
302 }
303 msg(MSG::VERBOSE) << " normalized f = " << tmp << endreq;
304 }
305 tmpVector.push_back(tmp);
306 sum += tmp;
307 }
308 if(2==idim) {
309 double valuex = m_lhVariableValues[iel].composites[icompo].atoms[0].value;
310 double valuey = m_lhVariableValues[iel].composites[icompo].atoms[1].value;
311 int binx = (tmpHisto->GetXaxis())->FindBin(valuex);
312 int biny = (tmpHisto->GetYaxis())->FindBin(valuey);
313 if(valuex >= tmpHisto->GetXaxis()->GetXmax()) binx = tmpHisto->GetXaxis()->GetNbins();
314 if(valuex <= tmpHisto->GetXaxis()->GetXmin()) binx = 1;
315 if(valuey >= tmpHisto->GetYaxis()->GetXmax()) biny = tmpHisto->GetYaxis()->GetNbins();
316 if(valuey <= tmpHisto->GetYaxis()->GetXmin()) biny = 1;
317 double tmp = tmpHisto->GetBinContent(binx, biny);
318 msg(MSG::VERBOSE) << "#BTAG# For hypothesis= " << m_hypotheses[ihyp]
319 << " (2D) actual value= " << valuex << " " << valuey
320 << " --> bin= " << binx << " " << biny << " f = " << tmp;
321 if(m_interpolate) {
322 tmp = HistoHelperRoot::Interpol2d(valuex, valuey, dynamic_cast<TH2*>(tmpHisto));
323 msg(MSG::VERBOSE) << " interpolated f = " << tmp;
324 }
325 if(m_normalizedProb) {
326 msg(MSG::VERBOSE) << " (normalized)" << endreq;
327 } else {
328 double binw = tmpHisto->GetXaxis()->GetBinWidth(binx)
329 * tmpHisto->GetYaxis()->GetBinWidth(biny);
330 msg(MSG::VERBOSE) << " binw= " << binw;
331 if(0==binw) {
332 msg(MSG::ERROR) << "bin width is 0" << endreq;
333 } else {
334 tmp /= binw;
335 }
336 msg(MSG::VERBOSE) << " normalized f = " << tmp << endreq;
337 }
338 tmpVector.push_back(tmp);
339 sum += tmp;
340 }
341 if(3==idim) {
342 double valuex = m_lhVariableValues[iel].composites[icompo].atoms[0].value;
343 double valuey = m_lhVariableValues[iel].composites[icompo].atoms[1].value;
344 double valuez = m_lhVariableValues[iel].composites[icompo].atoms[2].value;
345 int binx = (tmpHisto->GetXaxis())->FindBin(valuex);
346 int biny = (tmpHisto->GetYaxis())->FindBin(valuey);
347 int binz = (tmpHisto->GetZaxis())->FindBin(valuez);
348 if(valuex >= tmpHisto->GetXaxis()->GetXmax()) binx = tmpHisto->GetXaxis()->GetNbins();
349 if(valuex <= tmpHisto->GetXaxis()->GetXmin()) binx = 1;
350 if(valuey >= tmpHisto->GetYaxis()->GetXmax()) biny = tmpHisto->GetYaxis()->GetNbins();
351 if(valuey <= tmpHisto->GetYaxis()->GetXmin()) biny = 1;
352 if(valuez >= tmpHisto->GetZaxis()->GetXmax()) binz = tmpHisto->GetZaxis()->GetNbins();
353 if(valuez <= tmpHisto->GetZaxis()->GetXmin()) binz = 1;
354 double tmp = tmpHisto->GetBinContent(binx, biny, binz);
355 msg(MSG::VERBOSE) << "#BTAG# For hypothesis= " << m_hypotheses[ihyp]
356 << " (3D) actual value= " << valuex << " " << valuey << " " << valuez
357 << " --> bin= " << binx << " " << biny << " " << binz << " f = " << tmp;
358 if(m_interpolate) {
359 tmp = HistoHelperRoot::Interpol3d(valuex, valuey, valuez, dynamic_cast<TH3*>(tmpHisto));
360 msg(MSG::VERBOSE) << " interpolated f = " << tmp;
361 }
362 if(m_normalizedProb) {
363 msg(MSG::VERBOSE) << " (normalized)" << endreq;
364 } else {
365 double binw = tmpHisto->GetXaxis()->GetBinWidth(binx)
366 * tmpHisto->GetYaxis()->GetBinWidth(biny)
367 * tmpHisto->GetZaxis()->GetBinWidth(binz);
368 msg(MSG::VERBOSE) << " binw= " << binw;
369 if(0==binw) {
370 msg(MSG::ERROR) << "bin width is 0" << endreq;
371 } else {
372 tmp /= binw;
373 }
374 msg(MSG::VERBOSE) << " normalized f = " << tmp << endreq;
375 }
376 tmpVector.push_back(tmp);
377 sum += tmp;
378 }
379 if(idim>3 || idim<1 ) msg(MSG::WARNING) << "#BTAG# " << idim
380 << " is not a correct dimensionality for pdfs !"
381 << endreq;
382 }
383 }
384 unsigned int classCount(0);
385 for(std::vector<double>::iterator itr3 = tmpVector.begin();
386 itr3 != tmpVector.end(); ++itr3) {
387 if(sum != 0.) {
388 msg(MSG::VERBOSE) << "#BTAG# pb + pu = " << sum << endreq;
389 double p = (*itr3);
390 if(m_normalizedProb) p /= sum;
391 probDensityPerEventClassAllVariables[classCount] *= p;
392 } else {
393 msg(MSG::WARNING) << "#BTAG# Empty bins for all hypothesis... "
394 << "The discriminating variables are not taken into "
395 << "account in this jet." << endreq;
396 msg(MSG::WARNING) << "#BTAG# Please check your inputs" << endreq;
397 }
398 msg(MSG::VERBOSE) << "#BTAG# probDensity= "
399 << probDensityPerEventClassAllVariables[classCount]
400 << " ic= " << classCount << endreq;
401 classCount++;
402 }
403 msg(MSG::VERBOSE) << "#BTAG# Final probDensity= "
404 << probDensityPerEventClassAllVariables << endreq;
405 }
406 }
407 msg(MSG::VERBOSE) << "#BTAG# Ending ..." << endreq;
408 m_likelihoodVector=probDensityPerEventClassAllVariables;
409 return m_likelihoodVector;
410 }
411
412 std::vector<double> NewLikelihoodTool::tagLikelihood() {
413 return m_likelihoodVector;
414 }
415
416 double NewLikelihoodTool::getEff(const std::string& hypo, const std::string& hname, const std::string& suffix) {
417 double eff(0.);
418 TH1* numH = this->prepareHistogram(hypo, hname+"Eff"+suffix);
419 TH1* denH = this->prepareHistogram(hypo, hname+"Norm"+suffix);
420 if(numH && denH) {
421 double nnum = numH->GetEntries();
422 double nden = denH->GetEntries();
423 if(nden!=0) {
424 eff = nnum / nden;
425 msg(MSG::VERBOSE) << "#BTAG# Efficiency for " << hypo << "#" << hname << " "
426 << suffix << ": " << nnum << "/" << nden << "= " << eff << endreq;
427 } else {
428 msg(MSG::WARNING) << "#BTAG# Problem with Efficiency for " << hypo << "#" << hname << " "
429 << suffix << ": " << nnum << "/" << nden << "= " << eff << endreq;
430 }
431 } else {
432 msg(MSG::WARNING) << "#BTAG# Unknown histogram for Efficiency for " << hypo
433 << "#" << hname << " " << suffix << endreq;
434 }
435 return eff;
436 }
437
438
439 }
| Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems |
|
This page was automatically generated by the
LXR engine.
|
|