CMS 3D CMS Logo

PFRecoTauDiscriminationByMVAIsolationRun2.cc
Go to the documentation of this file.
1 
11 
18 
20 
23 
32 
35 
36 #include <TFile.h>
37 
38 #include <iostream>
39 
40 using namespace reco;
41 
42 namespace {
43  const GBRForest* loadMVAfromFile(const edm::FileInPath& inputFileName,
44  const std::string& mvaName,
45  std::vector<TFile*>& inputFilesToDelete) {
46  if (inputFileName.location() == edm::FileInPath::Unknown)
47  throw cms::Exception("PFRecoTauDiscriminationByIsolationMVA2::loadMVA")
48  << " Failed to find File = " << inputFileName << " !!\n";
49  TFile* inputFile = new TFile(inputFileName.fullPath().data());
50 
51  //const GBRForest* mva = dynamic_cast<GBRForest*>(inputFile->Get(mvaName.data())); // CV: dynamic_cast<GBRForest*> fails for some reason ?!
52  const GBRForest* mva = (GBRForest*)inputFile->Get(mvaName.data());
53  if (!mva)
54  throw cms::Exception("PFRecoTauDiscriminationByIsolationMVA2::loadMVA")
55  << " Failed to load MVA = " << mvaName.data() << " from file = " << inputFileName.fullPath().data()
56  << " !!\n";
57 
58  inputFilesToDelete.push_back(inputFile);
59 
60  return mva;
61  }
62 } // namespace
63 
64 namespace reco {
65  namespace tau {
66 
68  public:
71  moduleLabel_(cfg.getParameter<std::string>("@module_label")),
72  mvaReader_(nullptr),
73  mvaInput_(nullptr) {
74  mvaName_ = cfg.getParameter<std::string>("mvaName");
75  loadMVAfromDB_ = cfg.getParameter<bool>("loadMVAfromDB");
76  if (!loadMVAfromDB_) {
77  inputFileName_ = cfg.getParameter<edm::FileInPath>("inputFileName");
78  } else {
80  }
81  std::string mvaOpt_string = cfg.getParameter<std::string>("mvaOpt");
82  if (mvaOpt_string == "oldDMwoLT")
84  else if (mvaOpt_string == "oldDMwLT")
86  else if (mvaOpt_string == "newDMwoLT")
88  else if (mvaOpt_string == "newDMwLT")
90  else if (mvaOpt_string == "DBoldDMwLT")
92  else if (mvaOpt_string == "DBnewDMwLT")
94  else if (mvaOpt_string == "PWoldDMwLT")
96  else if (mvaOpt_string == "PWnewDMwLT")
98  else if (mvaOpt_string == "DBoldDMwLTwGJ")
100  else if (mvaOpt_string == "DBnewDMwLTwGJ")
102  else
103  throw cms::Exception("PFRecoTauDiscriminationByMVAIsolationRun2")
104  << " Invalid Configuration Parameter 'mvaOpt' = " << mvaOpt_string << " !!\n";
105 
106  if (mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT)
107  mvaInput_ = new float[6];
108  else if (mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT)
109  mvaInput_ = new float[12];
112  mvaInput_ = new float[23];
113  else
114  assert(0);
115 
117  consumes<PFTauTIPAssociationByRef>(cfg.getParameter<edm::InputTag>("srcTauTransverseImpactParameters"));
118 
120  consumes<reco::TauDiscriminatorContainer>(cfg.getParameter<edm::InputTag>("srcBasicTauDiscriminators"));
121  input_id_name_suffix_ = cfg.getParameter<std::string>("inputIDNameSuffix");
122 
123  verbosity_ = cfg.getParameter<int>("verbosity");
124  }
125 
126  void beginEvent(const edm::Event&, const edm::EventSetup&) override;
127 
129 
131  if (!loadMVAfromDB_)
132  delete mvaReader_;
133  delete[] mvaInput_;
134  for (std::vector<TFile*>::iterator it = inputFilesToDelete_.begin(); it != inputFilesToDelete_.end(); ++it) {
135  delete (*it);
136  }
137  }
138 
139  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
140 
141  private:
143 
149  int mvaOpt_;
150  float* mvaInput_;
151 
156 
166 
168 
169  std::vector<TFile*> inputFilesToDelete_;
170 
172  };
173 
175  if (!mvaReader_) {
176  if (loadMVAfromDB_) {
178  } else {
180  }
181  }
182 
184 
186 
187  evt.getByToken(Tau_token, taus_);
188 
189  // load indices from input provenance config if process history changed, in particular for the first event
190  // skip missing IDs and leave treatment to produce/discriminate function
191  if (evt.processHistoryID() != phID_ && basicTauDiscriminators_.isValid()) {
192  phID_ = evt.processHistoryID();
193  const edm::Provenance* prov = basicTauDiscriminators_.provenance();
194  const std::vector<edm::ParameterSet> psetsFromProvenance =
195  edm::parameterSet(prov->stable(), evt.processHistory())
196  .getParameter<std::vector<edm::ParameterSet>>("IDdefinitions");
197  for (uint i = 0; i < psetsFromProvenance.size(); i++) {
198  if (psetsFromProvenance[i].getParameter<std::string>("IDname") == "ChargedIsoPtSum" + input_id_name_suffix_)
200  else if (psetsFromProvenance[i].getParameter<std::string>("IDname") ==
201  "NeutralIsoPtSum" + input_id_name_suffix_)
203  else if (psetsFromProvenance[i].getParameter<std::string>("IDname") == "PUcorrPtSum" + input_id_name_suffix_)
205  else if (psetsFromProvenance[i].getParameter<std::string>("IDname") ==
206  "PhotonPtSumOutsideSignalCone" + input_id_name_suffix_)
208  else if (psetsFromProvenance[i].getParameter<std::string>("IDname") ==
209  "TauFootprintCorrection" + input_id_name_suffix_)
211  }
212  }
213  }
214 
216  const PFTauRef& tau) const {
218  result.rawValues = {-1.};
219 
220  // CV: computation of MVA value requires presence of leading charged hadron
221  if (tau->leadChargedHadrCand().isNull())
222  return 0.;
223 
224  int tauDecayMode = tau->decayMode();
225 
227  mvaOpt_ == kDBoldDMwLTwGJ) &&
228  (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10)) ||
230  mvaOpt_ == kDBnewDMwLTwGJ) &&
231  (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 ||
232  tauDecayMode == 10 || tauDecayMode == 11))) {
233  auto const rawValues = (*basicTauDiscriminators_)[tau].rawValues;
234  float chargedIsoPtSum = rawValues.at(chargedIsoPtSum_index_);
235  float neutralIsoPtSum = rawValues.at(neutralIsoPtSum_index_);
236  float puCorrPtSum = rawValues.at(pucorrPtSum_index_);
237  float photonPtSumOutsideSignalCone = rawValues.at(photonPtSumOutsideSignalCone_index_);
238  float footprintCorrection = rawValues.at(footprintCorrection_index_);
239 
240  const reco::PFTauTransverseImpactParameter& tauLifetimeInfo = *(*tauLifetimeInfos)[tau];
241 
242  float decayDistX = tauLifetimeInfo.flightLength().x();
243  float decayDistY = tauLifetimeInfo.flightLength().y();
244  float decayDistZ = tauLifetimeInfo.flightLength().z();
245  float decayDistMag = std::sqrt(decayDistX * decayDistX + decayDistY * decayDistY + decayDistZ * decayDistZ);
246 
248  float ptWeightedDetaStrip = reco::tau::pt_weighted_deta_strip(*tau, tauDecayMode);
249  float ptWeightedDphiStrip = reco::tau::pt_weighted_dphi_strip(*tau, tauDecayMode);
250  float ptWeightedDrSignal = reco::tau::pt_weighted_dr_signal(*tau, tauDecayMode);
251  float ptWeightedDrIsolation = reco::tau::pt_weighted_dr_iso(*tau, tauDecayMode);
252  float leadingTrackChi2 = reco::tau::lead_track_chi2(*tau);
253  float eRatio = reco::tau::eratio(*tau);
254 
255  // Difference between measured and maximally allowed Gottfried-Jackson angle
256  float gjAngleDiff = -999;
257  if (tauDecayMode == 10) {
258  double mTau = 1.77682;
259  double mAOne = tau->p4().M();
260  double pAOneMag = tau->p();
261  double argumentThetaGJmax = (std::pow(mTau, 2) - std::pow(mAOne, 2)) / (2 * mTau * pAOneMag);
262  double argumentThetaGJmeasured =
263  (tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ) /
264  (pAOneMag * decayDistMag);
265  if (std::abs(argumentThetaGJmax) <= 1. && std::abs(argumentThetaGJmeasured) <= 1.) {
266  double thetaGJmax = std::asin(argumentThetaGJmax);
267  double thetaGJmeasured = std::acos(argumentThetaGJmeasured);
268  gjAngleDiff = thetaGJmeasured - thetaGJmax;
269  }
270  }
271 
272  if (mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT) {
273  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
274  mvaInput_[1] = std::abs((float)tau->eta());
275  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
276  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum - 0.125f * puCorrPtSum));
277  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
278  mvaInput_[5] = tauDecayMode;
279  } else if (mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT) {
280  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
281  mvaInput_[1] = std::abs((float)tau->eta());
282  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
283  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum - 0.125f * puCorrPtSum));
284  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
285  mvaInput_[5] = tauDecayMode;
286  mvaInput_[6] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
287  mvaInput_[7] = std::sqrt(std::abs(std::min(1.f, (float)tauLifetimeInfo.dxy())));
288  mvaInput_[8] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
289  mvaInput_[9] = (tauLifetimeInfo.hasSecondaryVertex()) ? 1. : 0.;
290  mvaInput_[10] = std::sqrt(decayDistMag);
291  mvaInput_[11] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
292  } else if (mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kDBnewDMwLT) {
293  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
294  mvaInput_[1] = std::abs((float)tau->eta());
295  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
296  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
297  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
298  mvaInput_[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
299  mvaInput_[6] = tauDecayMode;
300  mvaInput_[7] = std::min(30.f, nPhoton);
301  mvaInput_[8] = std::min(0.5f, ptWeightedDetaStrip);
302  mvaInput_[9] = std::min(0.5f, ptWeightedDphiStrip);
303  mvaInput_[10] = std::min(0.5f, ptWeightedDrSignal);
304  mvaInput_[11] = std::min(0.5f, ptWeightedDrIsolation);
305  mvaInput_[12] = std::min(100.f, leadingTrackChi2);
306  mvaInput_[13] = std::min(1.f, eRatio);
307  mvaInput_[14] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
308  mvaInput_[15] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.dxy())));
309  mvaInput_[16] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
310  mvaInput_[17] = std::copysign(+1.f, (float)tauLifetimeInfo.ip3d());
311  mvaInput_[18] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.ip3d())));
312  mvaInput_[19] = std::min(10.f, std::abs((float)tauLifetimeInfo.ip3d_Sig()));
313  mvaInput_[20] = (tauLifetimeInfo.hasSecondaryVertex()) ? 1. : 0.;
314  mvaInput_[21] = std::sqrt(decayDistMag);
315  mvaInput_[22] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
316  } else if (mvaOpt_ == kPWoldDMwLT || mvaOpt_ == kPWnewDMwLT) {
317  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
318  mvaInput_[1] = std::abs((float)tau->eta());
319  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
320  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
321  mvaInput_[4] = std::log(std::max(1.e-2f, footprintCorrection));
322  mvaInput_[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
323  mvaInput_[6] = tauDecayMode;
324  mvaInput_[7] = std::min(30.f, nPhoton);
325  mvaInput_[8] = std::min(0.5f, ptWeightedDetaStrip);
326  mvaInput_[9] = std::min(0.5f, ptWeightedDphiStrip);
327  mvaInput_[10] = std::min(0.5f, ptWeightedDrSignal);
328  mvaInput_[11] = std::min(0.5f, ptWeightedDrIsolation);
329  mvaInput_[12] = std::min(100.f, leadingTrackChi2);
330  mvaInput_[13] = std::min(1.f, eRatio);
331  mvaInput_[14] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
332  mvaInput_[15] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.dxy())));
333  mvaInput_[16] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
334  mvaInput_[17] = std::copysign(+1.f, (float)tauLifetimeInfo.ip3d());
335  mvaInput_[18] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.ip3d())));
336  mvaInput_[19] = std::min(10.f, std::abs((float)tauLifetimeInfo.ip3d_Sig()));
337  mvaInput_[20] = (tauLifetimeInfo.hasSecondaryVertex()) ? 1. : 0.;
338  mvaInput_[21] = std::sqrt(decayDistMag);
339  mvaInput_[22] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
340  } else if (mvaOpt_ == kDBoldDMwLTwGJ || mvaOpt_ == kDBnewDMwLTwGJ) {
341  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
342  mvaInput_[1] = std::abs((float)tau->eta());
343  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
344  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
345  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
346  mvaInput_[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
347  mvaInput_[6] = tauDecayMode;
348  mvaInput_[7] = std::min(30.f, nPhoton);
349  mvaInput_[8] = std::min(0.5f, ptWeightedDetaStrip);
350  mvaInput_[9] = std::min(0.5f, ptWeightedDphiStrip);
351  mvaInput_[10] = std::min(0.5f, ptWeightedDrSignal);
352  mvaInput_[11] = std::min(0.5f, ptWeightedDrIsolation);
353  mvaInput_[12] = std::min(1.f, eRatio);
354  mvaInput_[13] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
355  mvaInput_[14] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.dxy())));
356  mvaInput_[15] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
357  mvaInput_[16] = std::copysign(+1.f, (float)tauLifetimeInfo.ip3d());
358  mvaInput_[17] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.ip3d())));
359  mvaInput_[18] = std::min(10.f, std::abs((float)tauLifetimeInfo.ip3d_Sig()));
360  mvaInput_[19] = (tauLifetimeInfo.hasSecondaryVertex()) ? 1. : 0.;
361  mvaInput_[20] = std::sqrt(decayDistMag);
362  mvaInput_[21] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
363  mvaInput_[22] = std::max(-1.f, gjAngleDiff);
364  }
365 
366  double mvaValue = mvaReader_->GetClassifier(mvaInput_);
367  if (verbosity_) {
368  edm::LogPrint("PFTauDiscByMVAIsol2") << "<PFRecoTauDiscriminationByMVAIsolationRun2::discriminate>:";
369  edm::LogPrint("PFTauDiscByMVAIsol2") << " tau: Pt = " << tau->pt() << ", eta = " << tau->eta();
370  edm::LogPrint("PFTauDiscByMVAIsol2") << " isolation: charged = " << chargedIsoPtSum
371  << ", neutral = " << neutralIsoPtSum << ", PUcorr = " << puCorrPtSum;
372  edm::LogPrint("PFTauDiscByMVAIsol2") << " decay mode = " << tauDecayMode;
373  edm::LogPrint("PFTauDiscByMVAIsol2") << " impact parameter: distance = " << tauLifetimeInfo.dxy()
374  << ", significance = " << tauLifetimeInfo.dxy_Sig();
375  edm::LogPrint("PFTauDiscByMVAIsol2")
376  << " has decay vertex = " << tauLifetimeInfo.hasSecondaryVertex() << ":"
377  << " distance = " << decayDistMag << ", significance = " << tauLifetimeInfo.flightLengthSig();
378  edm::LogPrint("PFTauDiscByMVAIsol2") << "--> mvaValue = " << mvaValue;
379  }
380  result.rawValues.at(0) = mvaValue;
381  }
382  return result;
383  }
384 
386  // pfRecoTauDiscriminationByMVAIsolationRun2
388 
389  desc.add<std::string>("mvaName", "tauIdMVAnewDMwLT");
390  desc.add<bool>("loadMVAfromDB", true);
391  desc.addOptional<edm::FileInPath>("inputFileName");
392  desc.add<std::string>("mvaOpt", "newDMwLT");
393 
394  desc.add<edm::InputTag>("srcTauTransverseImpactParameters", edm::InputTag(""));
395  desc.add<edm::InputTag>("srcBasicTauDiscriminators", edm::InputTag("hpsPFTauBasicDiscriminators"));
396  desc.add<std::string>("inputIDNameSuffix", "");
397 
398  desc.add<int>("verbosity", 0);
399 
400  fillProducerDescriptions(desc); // inherited from the base
401 
402  descriptions.add("pfRecoTauDiscriminationByMVAIsolationRun2", desc);
403  }
404 
406 
407  } // namespace tau
408 } // namespace reco
edm::EDGetTokenT< reco::TauDiscriminatorContainer > basicTauDiscriminators_token
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
unsigned int n_photons_total(const reco::PFTau &tau)
return total number of pf photon candidates with pT>500 MeV, which are associated to signal ...
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
float pt_weighted_dr_signal(const reco::PFTau &tau, int dm)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:526
assert(be >=bs)
ParameterSet const & parameterSet(StableProvenance const &provenance, ProcessHistory const &history)
Definition: Provenance.cc:11
float pt_weighted_deta_strip(const reco::PFTau &tau, int dm)
T sqrt(T t)
Definition: SSEVec.h:23
static void fillProducerDescriptions(edm::ParameterSetDescription &desc)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
float pt_weighted_dr_iso(const reco::PFTau &tau, int dm)
double f[11][100]
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Log< level::Warning, true > LogPrint
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
edm::AssociationVector< reco::PFTauRefProd, std::vector< reco::PFTauTransverseImpactParameterRef > > PFTauTIPAssociationByRef
edm::EDGetTokenT< PFTauTIPAssociationByRef > TauTransverseImpactParameters_token
double GetClassifier(const float *vector) const
Definition: GBRForest.h:33
float pt_weighted_dphi_strip(const reco::PFTau &tau, int dm)
float eratio(const reco::PFTau &tau)
return ratio of energy in ECAL over sum of energy in ECAL and HCAL
void add(std::string const &label, ParameterSetDescription const &psetDescription)
edm::EDGetTokenT< TauCollection > Tau_token
void beginEvent(const edm::Event &, const edm::EventSetup &) override
fixed size matrix
float lead_track_chi2(const reco::PFTau &tau)
return chi2 of the leading track ==> deprecated? <==
reco::SingleTauDiscriminatorContainer discriminate(const PFTauRef &) const override
ProcessHistoryID const & processHistoryID() const
Definition: Event.cc:116
StableProvenance const & stable() const
Definition: Provenance.h:42
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
ProcessHistory const & processHistory() const override
Definition: Event.cc:250