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 {
43  const GBRForest* loadMVAfromFile(const edm::FileInPath& inputFileName, const std::string& mvaName, std::vector<TFile*>& inputFilesToDelete)
44  {
45  if ( inputFileName.location() == edm::FileInPath::Unknown ) throw cms::Exception("PFRecoTauDiscriminationByIsolationMVA2::loadMVA")
46  << " Failed to find File = " << inputFileName << " !!\n";
47  TFile* inputFile = new TFile(inputFileName.fullPath().data());
48 
49  //const GBRForest* mva = dynamic_cast<GBRForest*>(inputFile->Get(mvaName.data())); // CV: dynamic_cast<GBRForest*> fails for some reason ?!
50  const GBRForest* mva = (GBRForest*)inputFile->Get(mvaName.data());
51  if ( !mva )
52  throw cms::Exception("PFRecoTauDiscriminationByIsolationMVA2::loadMVA")
53  << " Failed to load MVA = " << mvaName.data() << " from file = " << inputFileName.fullPath().data() << " !!\n";
54 
55  inputFilesToDelete.push_back(inputFile);
56 
57  return mva;
58  }
59 
60  const GBRForest* loadMVAfromDB(const edm::EventSetup& es, const std::string& mvaName)
61  {
63  es.get<GBRWrapperRcd>().get(mvaName, mva);
64  return mva.product();
65  }
66 }
67 
68 namespace reco { namespace tau {
69 
71 {
72  public:
75  moduleLabel_(cfg.getParameter<std::string>("@module_label")),
79  {
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" ) mvaOpt_ = kOldDMwoLT;
87  else if ( mvaOpt_string == "oldDMwLT" ) mvaOpt_ = kOldDMwLT;
88  else if ( mvaOpt_string == "newDMwoLT" ) mvaOpt_ = kNewDMwoLT;
89  else if ( mvaOpt_string == "newDMwLT" ) mvaOpt_ = kNewDMwLT;
90  else if ( mvaOpt_string == "DBoldDMwLT" ) mvaOpt_ = kDBoldDMwLT;
91  else if ( mvaOpt_string == "DBnewDMwLT" ) mvaOpt_ = kDBnewDMwLT;
92  else if ( mvaOpt_string == "PWoldDMwLT" ) mvaOpt_ = kPWoldDMwLT;
93  else if ( mvaOpt_string == "PWnewDMwLT" ) mvaOpt_ = kPWnewDMwLT;
94  else if ( mvaOpt_string == "DBoldDMwLTwGJ" ) mvaOpt_ = kDBoldDMwLTwGJ;
95  else if ( mvaOpt_string == "DBnewDMwLTwGJ" ) mvaOpt_ = kDBnewDMwLTwGJ;
96  else throw cms::Exception("PFRecoTauDiscriminationByMVAIsolationRun2")
97  << " Invalid Configuration Parameter 'mvaOpt' = " << mvaOpt_string << " !!\n";
98 
99  if ( mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT ) mvaInput_ = new float[6];
100  else if ( mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT ) mvaInput_ = new float[12];
101  else if ( mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kDBnewDMwLT ||
103  mvaOpt_ == kDBoldDMwLTwGJ || mvaOpt_ == kDBnewDMwLTwGJ) mvaInput_ = new float[23];
104  else assert(0);
105 
106  TauTransverseImpactParameters_token = consumes<PFTauTIPAssociationByRef>(cfg.getParameter<edm::InputTag>("srcTauTransverseImpactParameters"));
107 
108  ChargedIsoPtSum_token = consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcChargedIsoPtSum"));
109  NeutralIsoPtSum_token = consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcNeutralIsoPtSum"));
110  PUcorrPtSum_token = consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcPUcorrPtSum"));
111  PhotonPtSumOutsideSignalCone_token = consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcPhotonPtSumOutsideSignalCone"));
112  FootprintCorrection_token = consumes<reco::PFTauDiscriminator>(cfg.getParameter<edm::InputTag>("srcFootprintCorrection"));
113 
114  verbosity_ = cfg.getParameter<int>("verbosity");
115 
116  produces<PFTauDiscriminator>("category");
117  }
118 
119  void beginEvent(const edm::Event&, const edm::EventSetup&) override;
120 
121  double discriminate(const PFTauRef&) const override;
122 
123  void endEvent(edm::Event&) override;
124 
126  {
127  if(!loadMVAfromDB_) delete mvaReader_;
128  delete[] mvaInput_;
129  for ( std::vector<TFile*>::iterator it = inputFilesToDelete_.begin();
130  it != inputFilesToDelete_.end(); ++it ) {
131  delete (*it);
132  }
133  }
134 
135  static void fillDescriptions(edm::ConfigurationDescriptions & descriptions);
136  private:
137 
139 
144  int mvaOpt_;
145  float* mvaInput_;
146 
150 
161 
163  std::unique_ptr<PFTauDiscriminator> category_output_;
164 
165  std::vector<TFile*> inputFilesToDelete_;
166 
168 };
169 
171 {
172  if ( !mvaReader_ ) {
173  if ( loadMVAfromDB_ ) {
175  } else {
177  }
178  }
179 
181 
187 
188  evt.getByToken(Tau_token, taus_);
190 }
191 
193 {
194  // CV: define dummy category index in order to use RecoTauDiscriminantCutMultiplexer module to appy WP cuts
195  double category = 0.;
196  category_output_->setValue(tauIndex_, category);
197 
198  // CV: computation of MVA value requires presence of leading charged hadron
199  if ( tau->leadChargedHadrCand().isNull() ) return 0.;
200 
201  int tauDecayMode = tau->decayMode();
202 
204  && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10))
205  ||
207  && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 || tauDecayMode == 10 || tauDecayMode == 11))
208  ) {
209 
210  float chargedIsoPtSum = (*chargedIsoPtSums_)[tau];
211  float neutralIsoPtSum = (*neutralIsoPtSums_)[tau];
212  float puCorrPtSum = (*puCorrPtSums_)[tau];
213  float photonPtSumOutsideSignalCone = (*photonPtSumOutsideSignalCone_)[tau];
214  float footprintCorrection = (*footprintCorrection_)[tau];
215 
216  const reco::PFTauTransverseImpactParameter& tauLifetimeInfo = *(*tauLifetimeInfos)[tau];
217 
218  float decayDistX = tauLifetimeInfo.flightLength().x();
219  float decayDistY = tauLifetimeInfo.flightLength().y();
220  float decayDistZ = tauLifetimeInfo.flightLength().z();
221  float decayDistMag = std::sqrt(decayDistX*decayDistX + decayDistY*decayDistY + decayDistZ*decayDistZ);
222 
223  float nPhoton = (float)reco::tau::n_photons_total(*tau);
224  float ptWeightedDetaStrip = reco::tau::pt_weighted_deta_strip(*tau, tauDecayMode);
225  float ptWeightedDphiStrip = reco::tau::pt_weighted_dphi_strip(*tau, tauDecayMode);
226  float ptWeightedDrSignal = reco::tau::pt_weighted_dr_signal(*tau, tauDecayMode);
227  float ptWeightedDrIsolation = reco::tau::pt_weighted_dr_iso(*tau, tauDecayMode);
228  float leadingTrackChi2 = reco::tau::lead_track_chi2(*tau);
229  float eRatio = reco::tau::eratio(*tau);
230 
231  // Difference between measured and maximally allowed Gottfried-Jackson angle
232  float gjAngleDiff = -999;
233  if ( tauDecayMode == 10 ) {
234  double mTau = 1.77682;
235  double mAOne = tau->p4().M();
236  double pAOneMag = tau->p();
237  double argumentThetaGJmax = (std::pow(mTau,2) - std::pow(mAOne,2) ) / ( 2 * mTau * pAOneMag );
238  double argumentThetaGJmeasured = ( tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ ) / ( pAOneMag * decayDistMag );
239  if ( std::abs(argumentThetaGJmax) <= 1. && std::abs(argumentThetaGJmeasured) <= 1. ) {
240  double thetaGJmax = std::asin( argumentThetaGJmax );
241  double thetaGJmeasured = std::acos( argumentThetaGJmeasured );
242  gjAngleDiff = thetaGJmeasured - thetaGJmax;
243  }
244  }
245 
246  if ( mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT ) {
247  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
248  mvaInput_[1] = std::abs((float)tau->eta());
249  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
250  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum - 0.125f*puCorrPtSum));
251  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
252  mvaInput_[5] = tauDecayMode;
253  } else if ( mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT ) {
254  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
255  mvaInput_[1] = std::abs((float)tau->eta());
256  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
257  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum - 0.125f*puCorrPtSum));
258  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
259  mvaInput_[5] = tauDecayMode;
260  mvaInput_[6] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
261  mvaInput_[7] = std::sqrt(std::abs(std::min(1.f, (float)tauLifetimeInfo.dxy())));
262  mvaInput_[8] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
263  mvaInput_[9] = ( tauLifetimeInfo.hasSecondaryVertex() ) ? 1. : 0.;
264  mvaInput_[10] = std::sqrt(decayDistMag);
265  mvaInput_[11] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
266  } else if ( mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kDBnewDMwLT ) {
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));
271  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
272  mvaInput_[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
273  mvaInput_[6] = tauDecayMode;
274  mvaInput_[7] = std::min(30.f, nPhoton);
275  mvaInput_[8] = std::min(0.5f, ptWeightedDetaStrip);
276  mvaInput_[9] = std::min(0.5f, ptWeightedDphiStrip);
277  mvaInput_[10] = std::min(0.5f, ptWeightedDrSignal);
278  mvaInput_[11] = std::min(0.5f, ptWeightedDrIsolation);
279  mvaInput_[12] = std::min(100.f, leadingTrackChi2);
280  mvaInput_[13] = std::min(1.f, eRatio);
281  mvaInput_[14] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
282  mvaInput_[15] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.dxy())));
283  mvaInput_[16] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
284  mvaInput_[17] = std::copysign(+1.f, (float)tauLifetimeInfo.ip3d());
285  mvaInput_[18] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.ip3d())));
286  mvaInput_[19] = std::min(10.f, std::abs((float)tauLifetimeInfo.ip3d_Sig()));
287  mvaInput_[20] = ( tauLifetimeInfo.hasSecondaryVertex() ) ? 1. : 0.;
288  mvaInput_[21] = std::sqrt(decayDistMag);
289  mvaInput_[22] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
290  } else if ( mvaOpt_ == kPWoldDMwLT || mvaOpt_ == kPWnewDMwLT ) {
291  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
292  mvaInput_[1] = std::abs((float)tau->eta());
293  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
294  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
295  mvaInput_[4] = std::log(std::max(1.e-2f, footprintCorrection));
296  mvaInput_[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
297  mvaInput_[6] = tauDecayMode;
298  mvaInput_[7] = std::min(30.f, nPhoton);
299  mvaInput_[8] = std::min(0.5f, ptWeightedDetaStrip);
300  mvaInput_[9] = std::min(0.5f, ptWeightedDphiStrip);
301  mvaInput_[10] = std::min(0.5f, ptWeightedDrSignal);
302  mvaInput_[11] = std::min(0.5f, ptWeightedDrIsolation);
303  mvaInput_[12] = std::min(100.f, leadingTrackChi2);
304  mvaInput_[13] = std::min(1.f, eRatio);
305  mvaInput_[14] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
306  mvaInput_[15] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.dxy())));
307  mvaInput_[16] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
308  mvaInput_[17] = std::copysign(+1.f, (float)tauLifetimeInfo.ip3d());
309  mvaInput_[18] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.ip3d())));
310  mvaInput_[19] = std::min(10.f, std::abs((float)tauLifetimeInfo.ip3d_Sig()));
311  mvaInput_[20] = ( tauLifetimeInfo.hasSecondaryVertex() ) ? 1. : 0.;
312  mvaInput_[21] = std::sqrt(decayDistMag);
313  mvaInput_[22] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
314  } else if ( mvaOpt_ == kDBoldDMwLTwGJ || mvaOpt_ == kDBnewDMwLTwGJ ) {
315  mvaInput_[0] = std::log(std::max(1.f, (float)tau->pt()));
316  mvaInput_[1] = std::abs((float)tau->eta());
317  mvaInput_[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
318  mvaInput_[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
319  mvaInput_[4] = std::log(std::max(1.e-2f, puCorrPtSum));
320  mvaInput_[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
321  mvaInput_[6] = tauDecayMode;
322  mvaInput_[7] = std::min(30.f, nPhoton);
323  mvaInput_[8] = std::min(0.5f, ptWeightedDetaStrip);
324  mvaInput_[9] = std::min(0.5f, ptWeightedDphiStrip);
325  mvaInput_[10] = std::min(0.5f, ptWeightedDrSignal);
326  mvaInput_[11] = std::min(0.5f, ptWeightedDrIsolation);
327  mvaInput_[12] = std::min(1.f, eRatio);
328  mvaInput_[13] = std::copysign(+1.f, (float)tauLifetimeInfo.dxy());
329  mvaInput_[14] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.dxy())));
330  mvaInput_[15] = std::min(10.f, std::abs((float)tauLifetimeInfo.dxy_Sig()));
331  mvaInput_[16] = std::copysign(+1.f, (float)tauLifetimeInfo.ip3d());
332  mvaInput_[17] = std::sqrt(std::min(1.f, std::abs((float)tauLifetimeInfo.ip3d())));
333  mvaInput_[18] = std::min(10.f, std::abs((float)tauLifetimeInfo.ip3d_Sig()));
334  mvaInput_[19] = ( tauLifetimeInfo.hasSecondaryVertex() ) ? 1. : 0.;
335  mvaInput_[20] = std::sqrt(decayDistMag);
336  mvaInput_[21] = std::min(10.f, (float)tauLifetimeInfo.flightLengthSig());
337  mvaInput_[22] = std::max(-1.f, gjAngleDiff);
338  }
339 
340  double mvaValue = mvaReader_->GetClassifier(mvaInput_);
341  if ( verbosity_ ) {
342  edm::LogPrint("PFTauDiscByMVAIsol2") << "<PFRecoTauDiscriminationByMVAIsolationRun2::discriminate>:";
343  edm::LogPrint("PFTauDiscByMVAIsol2") << " tau: Pt = " << tau->pt() << ", eta = " << tau->eta();
344  edm::LogPrint("PFTauDiscByMVAIsol2") << " isolation: charged = " << chargedIsoPtSum << ", neutral = " << neutralIsoPtSum << ", PUcorr = " << puCorrPtSum;
345  edm::LogPrint("PFTauDiscByMVAIsol2") << " decay mode = " << tauDecayMode;
346  edm::LogPrint("PFTauDiscByMVAIsol2") << " impact parameter: distance = " << tauLifetimeInfo.dxy() << ", significance = " << tauLifetimeInfo.dxy_Sig();
347  edm::LogPrint("PFTauDiscByMVAIsol2") << " has decay vertex = " << tauLifetimeInfo.hasSecondaryVertex() << ":"
348  << " distance = " << decayDistMag << ", significance = " << tauLifetimeInfo.flightLengthSig();
349  edm::LogPrint("PFTauDiscByMVAIsol2") << "--> mvaValue = " << mvaValue;
350  }
351  return mvaValue;
352  } else {
353  return -1.;
354  }
355 }
356 
358 {
359  // add all category indices to event
360  evt.put(std::move(category_output_), "category");
361 }
362 
363 
364 void
366  // pfRecoTauDiscriminationByMVAIsolationRun2
368 
369  desc.add<std::string>("mvaName");
370  desc.add<bool>("loadMVAfromDB");
371  desc.addOptional<edm::FileInPath>("inputFileName");
372  desc.add<std::string>("mvaOpt");
373 
374  desc.add<edm::InputTag>("srcTauTransverseImpactParameters");
375  desc.add<edm::InputTag>("srcChargedIsoPtSum");
376  desc.add<edm::InputTag>("srcNeutralIsoPtSum");
377  desc.add<edm::InputTag>("srcPUcorrPtSum");
378 
379  desc.add<edm::InputTag>("srcPhotonPtSumOutsideSignalCone");
380  desc.add<edm::InputTag>("srcFootprintCorrection");
381 
382  desc.add<int>("verbosity", 0);
383 
384  fillProducerDescriptions(desc); // inherited from the base
385 
386  descriptions.add("pfRecoTauDiscriminationByMVAIsolationRun2", desc);
387 }
388 
389 
391 
392 }} //namespace
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:125
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:517
#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:18
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:248
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:71
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:36
T const * product() const
Definition: ESHandle.h:86
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
def move(src, dest)
Definition: eostools.py:511