CMS 3D CMS Logo

PFRecoTauDiscriminationByMVAIsolationRun2.cc
Go to the documentation of this file.
1 
11 
17 
19 
22 
30 
34 
35 #include <TFile.h>
36 
37 #include <iostream>
38 
39 using namespace reco;
40 
41 namespace {
42  const GBRForest* loadMVAfromFile(const edm::FileInPath& inputFileName,
43  const std::string& mvaName,
44  std::vector<TFile*>& inputFilesToDelete) {
45  if (inputFileName.location() == edm::FileInPath::Unknown)
46  throw cms::Exception("PFRecoTauDiscriminationByIsolationMVA2::loadMVA")
47  << " Failed to find File = " << inputFileName << " !!\n";
48  TFile* inputFile = new TFile(inputFileName.fullPath().data());
49 
50  //const GBRForest* mva = dynamic_cast<GBRForest*>(inputFile->Get(mvaName.data())); // CV: dynamic_cast<GBRForest*> fails for some reason ?!
51  const GBRForest* mva = (GBRForest*)inputFile->Get(mvaName.data());
52  if (!mva)
53  throw cms::Exception("PFRecoTauDiscriminationByIsolationMVA2::loadMVA")
54  << " Failed to load MVA = " << mvaName.data() << " from file = " << inputFileName.fullPath().data()
55  << " !!\n";
56 
57  inputFilesToDelete.push_back(inputFile);
58 
59  return mva;
60  }
61 
62  const GBRForest* loadMVAfromDB(const edm::EventSetup& es, const std::string& mvaName) {
64  es.get<GBRWrapperRcd>().get(mvaName, mva);
65  return mva.product();
66  }
67 } // namespace
68 
69 namespace reco {
70  namespace tau {
71 
73  public:
76  moduleLabel_(cfg.getParameter<std::string>("@module_label")),
80  mvaName_ = cfg.getParameter<std::string>("mvaName");
81  loadMVAfromDB_ = cfg.getParameter<bool>("loadMVAfromDB");
82  if (!loadMVAfromDB_) {
83  inputFileName_ = cfg.getParameter<edm::FileInPath>("inputFileName");
84  }
85  std::string mvaOpt_string = cfg.getParameter<std::string>("mvaOpt");
86  if (mvaOpt_string == "oldDMwoLT")
88  else if (mvaOpt_string == "oldDMwLT")
90  else if (mvaOpt_string == "newDMwoLT")
92  else if (mvaOpt_string == "newDMwLT")
94  else if (mvaOpt_string == "DBoldDMwLT")
96  else if (mvaOpt_string == "DBnewDMwLT")
98  else if (mvaOpt_string == "PWoldDMwLT")
100  else if (mvaOpt_string == "PWnewDMwLT")
102  else if (mvaOpt_string == "DBoldDMwLTwGJ")
104  else if (mvaOpt_string == "DBnewDMwLTwGJ")
106  else
107  throw cms::Exception("PFRecoTauDiscriminationByMVAIsolationRun2")
108  << " Invalid Configuration Parameter 'mvaOpt' = " << mvaOpt_string << " !!\n";
109 
110  if (mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT)
111  mvaInput_ = new float[6];
112  else if (mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT)
113  mvaInput_ = new float[12];
116  mvaInput_ = new float[23];
117  else
118  assert(0);
119 
121  consumes<PFTauTIPAssociationByRef>(cfg.getParameter<edm::InputTag>("srcTauTransverseImpactParameters"));
122 
124  consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcChargedIsoPtSum"));
126  consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcNeutralIsoPtSum"));
127  PUcorrPtSum_token = consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcPUcorrPtSum"));
129  consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcPhotonPtSumOutsideSignalCone"));
131  consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcFootprintCorrection"));
132 
133  verbosity_ = cfg.getParameter<int>("verbosity");
134 
135  produces<PFTauDiscriminator>("category");
136  }
137 
138  void beginEvent(const edm::Event&, const edm::EventSetup&) override;
139 
140  double discriminate(const PFTauRef&) const override;
141 
142  void endEvent(edm::Event&) override;
143 
145  if (!loadMVAfromDB_)
146  delete mvaReader_;
147  delete[] mvaInput_;
148  for (std::vector<TFile*>::iterator it = inputFilesToDelete_.begin(); it != inputFilesToDelete_.end(); ++it) {
149  delete (*it);
150  }
151  }
152 
153  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
154 
155  private:
157 
162  int mvaOpt_;
163  float* mvaInput_;
164 
169 
180 
182  std::unique_ptr<PFTauDiscriminator> category_output_;
183 
184  std::vector<TFile*> inputFilesToDelete_;
185 
187  };
188 
190  if (!mvaReader_) {
191  if (loadMVAfromDB_) {
193  } else {
195  }
196  }
197 
199 
205 
206  evt.getByToken(Tau_token, taus_);
208  }
209 
211  // CV: define dummy category index in order to use RecoTauDiscriminantCutMultiplexer module to appy WP cuts
212  double category = 0.;
213  category_output_->setValue(tauIndex_, category);
214 
215  // CV: computation of MVA value requires presence of leading charged hadron
216  if (tau->leadChargedHadrCand().isNull())
217  return 0.;
218 
219  int tauDecayMode = tau->decayMode();
220 
222  mvaOpt_ == kDBoldDMwLTwGJ) &&
223  (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10)) ||
225  mvaOpt_ == kDBnewDMwLTwGJ) &&
226  (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 ||
227  tauDecayMode == 10 || tauDecayMode == 11))) {
228  float chargedIsoPtSum = (*chargedIsoPtSums_)[tau];
229  float neutralIsoPtSum = (*neutralIsoPtSums_)[tau];
230  float puCorrPtSum = (*puCorrPtSums_)[tau];
231  float photonPtSumOutsideSignalCone = (*photonPtSumOutsideSignalCone_)[tau];
232  float footprintCorrection = (*footprintCorrection_)[tau];
233 
234  const reco::PFTauTransverseImpactParameter& tauLifetimeInfo = *(*tauLifetimeInfos)[tau];
235 
236  float decayDistX = tauLifetimeInfo.flightLength().x();
237  float decayDistY = tauLifetimeInfo.flightLength().y();
238  float decayDistZ = tauLifetimeInfo.flightLength().z();
239  float decayDistMag = std::sqrt(decayDistX * decayDistX + decayDistY * decayDistY + decayDistZ * decayDistZ);
240 
241  float nPhoton = (float)reco::tau::n_photons_total(*tau);
242  float ptWeightedDetaStrip = reco::tau::pt_weighted_deta_strip(*tau, tauDecayMode);
243  float ptWeightedDphiStrip = reco::tau::pt_weighted_dphi_strip(*tau, tauDecayMode);
244  float ptWeightedDrSignal = reco::tau::pt_weighted_dr_signal(*tau, tauDecayMode);
245  float ptWeightedDrIsolation = reco::tau::pt_weighted_dr_iso(*tau, tauDecayMode);
246  float leadingTrackChi2 = reco::tau::lead_track_chi2(*tau);
247  float eRatio = reco::tau::eratio(*tau);
248 
249  // Difference between measured and maximally allowed Gottfried-Jackson angle
250  float gjAngleDiff = -999;
251  if (tauDecayMode == 10) {
252  double mTau = 1.77682;
253  double mAOne = tau->p4().M();
254  double pAOneMag = tau->p();
255  double argumentThetaGJmax = (std::pow(mTau, 2) - std::pow(mAOne, 2)) / (2 * mTau * pAOneMag);
256  double argumentThetaGJmeasured =
257  (tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ) /
258  (pAOneMag * decayDistMag);
259  if (std::abs(argumentThetaGJmax) <= 1. && std::abs(argumentThetaGJmeasured) <= 1.) {
260  double thetaGJmax = std::asin(argumentThetaGJmax);
261  double thetaGJmeasured = std::acos(argumentThetaGJmeasured);
262  gjAngleDiff = thetaGJmeasured - thetaGJmax;
263  }
264  }
265 
266  if (mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT) {
267  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
268  mvaInput_[1] = std::abs((float)tau->eta());
269  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
270  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum - 0.125f * puCorrPtSum));
271  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
272  mvaInput_[5] = tauDecayMode;
273  } else if (mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT) {
274  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
275  mvaInput_[1] = std::abs((float)tau->eta());
276  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
277  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum - 0.125f * puCorrPtSum));
278  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
279  mvaInput_[5] = tauDecayMode;
280  mvaInput_[6] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
281  mvaInput_[7] = std::sqrt(std::abs(std::min(1.f, (float)tauLifetimeInfo.dxy())));
282  mvaInput_[8] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
283  mvaInput_[9] = (tauLifetimeInfo.hasSecondaryVertex()) ? 1. : 0.;
284  mvaInput_[10] = std::sqrt(decayDistMag);
285  mvaInput_[11] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
286  } else if (mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kDBnewDMwLT) {
287  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
288  mvaInput_[1] = std::abs((float)tau->eta());
289  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
290  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
291  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
292  mvaInput_[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
293  mvaInput_[6] = tauDecayMode;
294  mvaInput_[7] = std::min(30.f, nPhoton);
295  mvaInput_[8] = std::min(0.5f, ptWeightedDetaStrip);
296  mvaInput_[9] = std::min(0.5f, ptWeightedDphiStrip);
297  mvaInput_[10] = std::min(0.5f, ptWeightedDrSignal);
298  mvaInput_[11] = std::min(0.5f, ptWeightedDrIsolation);
299  mvaInput_[12] = std::min(100.f, leadingTrackChi2);
300  mvaInput_[13] = std::min(1.f, eRatio);
301  mvaInput_[14] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
302  mvaInput_[15] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.dxy())));
303  mvaInput_[16] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
304  mvaInput_[17] = std::copysign(+1.f, (float)tauLifetimeInfo.ip3d());
305  mvaInput_[18] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.ip3d())));
306  mvaInput_[19] = std::min(10.f, std::abs((float)tauLifetimeInfo.ip3d_Sig()));
307  mvaInput_[20] = (tauLifetimeInfo.hasSecondaryVertex()) ? 1. : 0.;
308  mvaInput_[21] = std::sqrt(decayDistMag);
309  mvaInput_[22] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
310  } else if (mvaOpt_ == kPWoldDMwLT || mvaOpt_ == kPWnewDMwLT) {
311  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
312  mvaInput_[1] = std::abs((float)tau->eta());
313  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
314  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
315  mvaInput_[4] = std::log(std::max(1.e-2f, footprintCorrection));
316  mvaInput_[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
317  mvaInput_[6] = tauDecayMode;
318  mvaInput_[7] = std::min(30.f, nPhoton);
319  mvaInput_[8] = std::min(0.5f, ptWeightedDetaStrip);
320  mvaInput_[9] = std::min(0.5f, ptWeightedDphiStrip);
321  mvaInput_[10] = std::min(0.5f, ptWeightedDrSignal);
322  mvaInput_[11] = std::min(0.5f, ptWeightedDrIsolation);
323  mvaInput_[12] = std::min(100.f, leadingTrackChi2);
324  mvaInput_[13] = std::min(1.f, eRatio);
325  mvaInput_[14] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
326  mvaInput_[15] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.dxy())));
327  mvaInput_[16] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
328  mvaInput_[17] = std::copysign(+1.f, (float)tauLifetimeInfo.ip3d());
329  mvaInput_[18] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.ip3d())));
330  mvaInput_[19] = std::min(10.f, std::abs((float)tauLifetimeInfo.ip3d_Sig()));
331  mvaInput_[20] = (tauLifetimeInfo.hasSecondaryVertex()) ? 1. : 0.;
332  mvaInput_[21] = std::sqrt(decayDistMag);
333  mvaInput_[22] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
334  } else if (mvaOpt_ == kDBoldDMwLTwGJ || mvaOpt_ == kDBnewDMwLTwGJ) {
335  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
336  mvaInput_[1] = std::abs((float)tau->eta());
337  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
338  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
339  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
340  mvaInput_[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
341  mvaInput_[6] = tauDecayMode;
342  mvaInput_[7] = std::min(30.f, nPhoton);
343  mvaInput_[8] = std::min(0.5f, ptWeightedDetaStrip);
344  mvaInput_[9] = std::min(0.5f, ptWeightedDphiStrip);
345  mvaInput_[10] = std::min(0.5f, ptWeightedDrSignal);
346  mvaInput_[11] = std::min(0.5f, ptWeightedDrIsolation);
347  mvaInput_[12] = std::min(1.f, eRatio);
348  mvaInput_[13] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
349  mvaInput_[14] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.dxy())));
350  mvaInput_[15] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
351  mvaInput_[16] = std::copysign(+1.f, (float)tauLifetimeInfo.ip3d());
352  mvaInput_[17] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.ip3d())));
353  mvaInput_[18] = std::min(10.f, std::abs((float)tauLifetimeInfo.ip3d_Sig()));
354  mvaInput_[19] = (tauLifetimeInfo.hasSecondaryVertex()) ? 1. : 0.;
355  mvaInput_[20] = std::sqrt(decayDistMag);
356  mvaInput_[21] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
357  mvaInput_[22] = std::max(-1.f, gjAngleDiff);
358  }
359 
360  double mvaValue = mvaReader_->GetClassifier(mvaInput_);
361  if (verbosity_) {
362  edm::LogPrint("PFTauDiscByMVAIsol2") << "<PFRecoTauDiscriminationByMVAIsolationRun2::discriminate>:";
363  edm::LogPrint("PFTauDiscByMVAIsol2") << " tau: Pt = " << tau->pt() << ", eta = " << tau->eta();
364  edm::LogPrint("PFTauDiscByMVAIsol2") << " isolation: charged = " << chargedIsoPtSum
365  << ", neutral = " << neutralIsoPtSum << ", PUcorr = " << puCorrPtSum;
366  edm::LogPrint("PFTauDiscByMVAIsol2") << " decay mode = " << tauDecayMode;
367  edm::LogPrint("PFTauDiscByMVAIsol2") << " impact parameter: distance = " << tauLifetimeInfo.dxy()
368  << ", significance = " << tauLifetimeInfo.dxy_Sig();
369  edm::LogPrint("PFTauDiscByMVAIsol2")
370  << " has decay vertex = " << tauLifetimeInfo.hasSecondaryVertex() << ":"
371  << " distance = " << decayDistMag << ", significance = " << tauLifetimeInfo.flightLengthSig();
372  edm::LogPrint("PFTauDiscByMVAIsol2") << "--> mvaValue = " << mvaValue;
373  }
374  return mvaValue;
375  } else {
376  return -1.;
377  }
378  }
379 
381  // add all category indices to event
382  evt.put(std::move(category_output_), "category");
383  }
384 
386  // pfRecoTauDiscriminationByMVAIsolationRun2
388 
389  desc.add<std::string>("mvaName");
390  desc.add<bool>("loadMVAfromDB");
391  desc.addOptional<edm::FileInPath>("inputFileName");
392  desc.add<std::string>("mvaOpt");
393 
394  desc.add<edm::InputTag>("srcTauTransverseImpactParameters");
395  desc.add<edm::InputTag>("srcChargedIsoPtSum");
396  desc.add<edm::InputTag>("srcNeutralIsoPtSum");
397  desc.add<edm::InputTag>("srcPUcorrPtSum");
398 
399  desc.add<edm::InputTag>("srcPhotonPtSumOutsideSignalCone");
400  desc.add<edm::InputTag>("srcFootprintCorrection");
401 
402  desc.add<int>("verbosity", 0);
403 
404  fillProducerDescriptions(desc); // inherited from the base
405 
406  descriptions.add("pfRecoTauDiscriminationByMVAIsolationRun2", desc);
407  }
408 
410 
411  } // namespace tau
412 } // namespace reco
T getParameter(std::string const &) const
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 ...
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:131
edm::EDGetTokenT< reco::PFTauDiscriminator > PhotonPtSumOutsideSignalCone_token
ParameterDescriptionBase * addOptional(U const &iLabel, T const &value)
edm::RefProd< TauCollection > TauRefProd
float pt_weighted_dr_signal(const reco::PFTau &tau, int dm)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:525
#define nullptr
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
static void fillProducerDescriptions(edm::ParameterSetDescription &desc)
float pt_weighted_deta_strip(const reco::PFTau &tau, int dm)
T sqrt(T t)
Definition: SSEVec.h:19
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]
T min(T a, T b)
Definition: MathUtil.h:58
LocationCode location() const
Where was the file found?
Definition: FileInPath.cc:161
ParameterDescriptionBase * add(U const &iLabel, T const &value)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
bool isNull() const
Checks for null.
Definition: Ref.h:235
edm::AssociationVector< reco::PFTauRefProd, std::vector< reco::PFTauTransverseImpactParameterRef > > PFTauTIPAssociationByRef
edm::EDGetTokenT< TauCollection > Tau_token
edm::EDGetTokenT< PFTauTIPAssociationByRef > TauTransverseImpactParameters_token
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)
void beginEvent(const edm::Event &, const edm::EventSetup &) override
fixed size matrix
T get() const
Definition: EventSetup.h:73
std::string fullPath() const
Definition: FileInPath.cc:163
float lead_track_chi2(const reco::PFTau &tau)
return chi2 of the leading track ==> deprecated? <==
double GetClassifier(const float *vector) const
Definition: GBRForest.h:34
T const * product() const
Definition: ESHandle.h:86
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:30
def move(src, dest)
Definition: eostools.py:511