File indexing completed on 2023-08-19 01:35:02
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
0018 #include "FWCore/Utilities/interface/InputTag.h"
0019 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0020 #include "DataFormats/TauReco/interface/PFTau.h"
0021 #include "FWCore/ParameterSet/interface/FileInPath.h"
0022
0023 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0024 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0025
0026 #include <memory>
0027
0028 #include "CondFormats/PhysicsToolsObjects/interface/PhysicsTGraphPayload.h"
0029 #include "CondFormats/DataRecord/interface/PhysicsTGraphPayloadRcd.h"
0030 #include "CondFormats/PhysicsToolsObjects/interface/PhysicsTFormulaPayload.h"
0031 #include "CondFormats/DataRecord/interface/PhysicsTFormulaPayloadRcd.h"
0032
0033 #include "TGraph.h"
0034 #include "TFormula.h"
0035 #include "TFile.h"
0036
0037 template <class TauType, class TauTypeRef, class ParentClass>
0038 class TauDiscriminantCutMultiplexerT : public ParentClass {
0039 public:
0040 explicit TauDiscriminantCutMultiplexerT(const edm::ParameterSet& pset);
0041
0042 ~TauDiscriminantCutMultiplexerT() override;
0043 reco::SingleTauDiscriminatorContainer discriminate(const TauTypeRef&) const override;
0044 void beginEvent(const edm::Event& event, const edm::EventSetup& eventSetup) override;
0045
0046 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047
0048 private:
0049 std::string moduleLabel_;
0050
0051 bool loadMVAfromDB_;
0052 edm::FileInPath inputFileName_;
0053
0054 struct DiscriminantCutEntry {
0055 DiscriminantCutEntry() : cutVariable_(), cutFunction_(), mode_(kUndefined) {}
0056 ~DiscriminantCutEntry() {}
0057 double cutValue_;
0058 std::string cutName_;
0059 edm::ESGetToken<PhysicsTGraphPayload, PhysicsTGraphPayloadRcd> cutToken_;
0060 std::unique_ptr<StringObjectFunction<TauType>> cutVariable_;
0061 std::unique_ptr<const TGraph> cutFunction_;
0062 enum { kUndefined, kFixedCut, kVariableCut };
0063 int mode_;
0064 };
0065 typedef std::map<int, std::vector<std::unique_ptr<DiscriminantCutEntry>>> DiscriminantCutMap;
0066 DiscriminantCutMap cuts_;
0067 uint n_raws_;
0068 int raw_discriminator_idx_ = -1;
0069 int raw_category_idx_ = -1;
0070
0071 std::string mvaOutputNormalizationName_;
0072 edm::ESGetToken<PhysicsTFormulaPayload, PhysicsTFormulaPayloadRcd> formulaToken_;
0073 std::unique_ptr<const TFormula> mvaOutput_normalization_;
0074
0075 bool isInitialized_;
0076
0077 edm::InputTag toMultiplex_;
0078 edm::Handle<reco::TauDiscriminatorContainer> toMultiplexHandle_;
0079 edm::EDGetTokenT<reco::TauDiscriminatorContainer> toMultiplex_token;
0080
0081 int verbosity_;
0082 };
0083
0084 namespace {
0085 std::unique_ptr<TFile> openInputFile(const edm::FileInPath& inputFileName) {
0086 if (inputFileName.location() == edm::FileInPath::Unknown) {
0087 throw cms::Exception("TauDiscriminantCutMultiplexerT::loadObjectFromFile")
0088 << " Failed to find File = " << inputFileName << " !!\n";
0089 }
0090 return std::make_unique<TFile>(inputFileName.fullPath().data());
0091 }
0092
0093 template <typename T>
0094 std::unique_ptr<const T> loadObjectFromFile(TFile& inputFile, const std::string& objectName) {
0095 const T* object = dynamic_cast<T*>(inputFile.Get(objectName.data()));
0096 if (!object)
0097 throw cms::Exception("TauDiscriminantCutMultiplexerT::loadObjectFromFile")
0098 << " Failed to load Object = " << objectName.data() << " from file = " << inputFile.GetName() << " !!\n";
0099
0100 return std::unique_ptr<const T>{static_cast<T*>(object->Clone())};
0101 }
0102
0103 std::unique_ptr<const TGraph> loadTGraphFromDB(
0104 const edm::EventSetup& es,
0105 const std::string& graphName,
0106 const edm::ESGetToken<PhysicsTGraphPayload, PhysicsTGraphPayloadRcd>& graphToken,
0107 const int& verbosity_ = 0) {
0108 if (verbosity_) {
0109 std::cout << "<loadTGraphFromDB>:" << std::endl;
0110 std::cout << " graphName = " << graphName << std::endl;
0111 }
0112 return std::make_unique<TGraph>(es.getData(graphToken));
0113 }
0114
0115 std::unique_ptr<TFormula> loadTFormulaFromDB(
0116 const edm::EventSetup& es,
0117 const std::string& formulaName,
0118 const edm::ESGetToken<PhysicsTFormulaPayload, PhysicsTFormulaPayloadRcd>& formulaToken_,
0119 const TString& newName,
0120 const int& verbosity_ = 0) {
0121 if (verbosity_) {
0122 std::cout << "<loadTFormulaFromDB>:" << std::endl;
0123 std::cout << " formulaName = " << formulaName << std::endl;
0124 }
0125 auto const& formulaPayload = es.getData(formulaToken_);
0126
0127 if (formulaPayload.formulas().size() == 1 && formulaPayload.limits().size() == 1) {
0128 return std::make_unique<TFormula>(newName, formulaPayload.formulas().at(0).data());
0129 } else {
0130 throw cms::Exception("TauDiscriminantCutMultiplexerT::loadTFormulaFromDB")
0131 << "Failed to load TFormula = " << formulaName << " from Database !!\n";
0132 }
0133 return std::unique_ptr<TFormula>{};
0134 }
0135 }
0136
0137 template <class TauType, class TauTypeRef, class ParentClass>
0138 TauDiscriminantCutMultiplexerT<TauType, TauTypeRef, ParentClass>::TauDiscriminantCutMultiplexerT(
0139 const edm::ParameterSet& cfg)
0140 : ParentClass(cfg),
0141 moduleLabel_(cfg.getParameter<std::string>("@module_label")),
0142 mvaOutput_normalization_(),
0143 isInitialized_(false) {
0144 toMultiplex_ = cfg.getParameter<edm::InputTag>("toMultiplex");
0145 toMultiplex_token = this->template consumes<reco::TauDiscriminatorContainer>(toMultiplex_);
0146
0147 verbosity_ = cfg.getParameter<int>("verbosity");
0148
0149 mvaOutputNormalizationName_ = cfg.getParameter<std::string>("mvaOutput_normalization");
0150
0151 loadMVAfromDB_ = cfg.getParameter<bool>("loadMVAfromDB");
0152 if (!loadMVAfromDB_) {
0153 inputFileName_ = cfg.getParameter<edm::FileInPath>("inputFileName");
0154 } else if (not mvaOutputNormalizationName_.empty()) {
0155 formulaToken_ = this->esConsumes(edm::ESInputTag{"", mvaOutputNormalizationName_});
0156 }
0157 if (verbosity_)
0158 std::cout << moduleLabel_ << " loadMVA = " << loadMVAfromDB_ << std::endl;
0159
0160
0161 typedef std::vector<edm::ParameterSet> VPSet;
0162 typedef std::vector<std::string> VString;
0163 typedef std::vector<double> VDouble;
0164 VString rawValueConfig = cfg.getParameter<VString>("rawValues");
0165 for (uint i = 0; i < rawValueConfig.size(); i++) {
0166 if (rawValueConfig[i] == "discriminator")
0167 raw_discriminator_idx_ = i;
0168 else if (rawValueConfig[i] == "category")
0169 raw_category_idx_ = i;
0170 else
0171 throw cms::Exception("TauDiscriminantCutMultiplexerT")
0172 << " Configuration Parameter 'rawValues' containes unknown values. Must be 'discriminator' or 'category'!!\n";
0173 }
0174 n_raws_ = rawValueConfig.size();
0175 VPSet mapping = cfg.getParameter<VPSet>("mapping");
0176 for (auto const& mappingEntry : mapping) {
0177 unsigned category = mappingEntry.getParameter<uint32_t>("category");
0178 std::vector<std::unique_ptr<DiscriminantCutEntry>> cutWPs;
0179 std::string categoryname = mappingEntry.getParameter<std::string>("cut");
0180 bool localWPs = false;
0181 bool wpsAsDouble = false;
0182 if (mappingEntry.exists("workingPoints")) {
0183 localWPs = true;
0184 if (mappingEntry.existsAs<VDouble>("workingPoints")) {
0185 wpsAsDouble = true;
0186 } else if (mappingEntry.existsAs<VString>("workingPoints")) {
0187 wpsAsDouble = false;
0188 } else {
0189 throw cms::Exception("TauDiscriminantCutMultiplexerT")
0190 << " Configuration Parameter 'workingPoints' must be filled with cms.string or cms.double!!\n";
0191 }
0192 } else if (cfg.exists("workingPoints")) {
0193 localWPs = false;
0194 if (cfg.existsAs<VDouble>("workingPoints")) {
0195 wpsAsDouble = true;
0196 } else if (cfg.existsAs<VString>("workingPoints")) {
0197 wpsAsDouble = false;
0198 } else {
0199 throw cms::Exception("TauDiscriminantCutMultiplexerT")
0200 << " Configuration Parameter 'workingPoints' must be filled with cms.string or cms.double!!\n";
0201 }
0202 } else {
0203 throw cms::Exception("TauDiscriminantCutMultiplexerT")
0204 << " Undefined Configuration Parameter 'workingPoints' !!\n";
0205 }
0206 if (wpsAsDouble) {
0207 VDouble workingPoints;
0208 if (localWPs)
0209 workingPoints = mappingEntry.getParameter<VDouble>("workingPoints");
0210 else
0211 workingPoints = cfg.getParameter<VDouble>("workingPoints");
0212 for (auto const& wp : workingPoints) {
0213 auto cut = std::make_unique<DiscriminantCutEntry>();
0214 cut->cutValue_ = wp;
0215 cut->mode_ = DiscriminantCutEntry::kFixedCut;
0216 cutWPs.push_back(std::move(cut));
0217 }
0218 } else {
0219 VString workingPoints;
0220 if (localWPs)
0221 workingPoints = mappingEntry.getParameter<VString>("workingPoints");
0222 else
0223 workingPoints = cfg.getParameter<VString>("workingPoints");
0224 for (auto const& wp : workingPoints) {
0225 auto cut = std::make_unique<DiscriminantCutEntry>();
0226 cut->cutName_ = categoryname + wp;
0227 if (loadMVAfromDB_) {
0228 cut->cutToken_ = this->esConsumes(edm::ESInputTag{"", cut->cutName_});
0229 }
0230 std::string cutVariable_string = mappingEntry.getParameter<std::string>("variable");
0231 cut->cutVariable_.reset(new StringObjectFunction<TauType>(cutVariable_string));
0232 cut->mode_ = DiscriminantCutEntry::kVariableCut;
0233 cutWPs.push_back(std::move(cut));
0234 }
0235 }
0236 cuts_[category] = std::move(cutWPs);
0237 }
0238
0239 verbosity_ = cfg.getParameter<int>("verbosity");
0240 if (verbosity_)
0241 std::cout << "constructed " << moduleLabel_ << std::endl;
0242 }
0243
0244 template <class TauType, class TauTypeRef, class ParentClass>
0245 TauDiscriminantCutMultiplexerT<TauType, TauTypeRef, ParentClass>::~TauDiscriminantCutMultiplexerT() {}
0246
0247 template <class TauType, class TauTypeRef, class ParentClass>
0248 void TauDiscriminantCutMultiplexerT<TauType, TauTypeRef, ParentClass>::beginEvent(const edm::Event& evt,
0249 const edm::EventSetup& es) {
0250 if (verbosity_)
0251 std::cout << " begin! " << moduleLabel_ << " " << isInitialized_ << std::endl;
0252 if (!isInitialized_) {
0253
0254
0255 std::unique_ptr<TFile> inputFile;
0256 if (!mvaOutputNormalizationName_.empty()) {
0257 if (!loadMVAfromDB_) {
0258 inputFile = openInputFile(inputFileName_);
0259 mvaOutput_normalization_ = loadObjectFromFile<TFormula>(*inputFile, mvaOutputNormalizationName_);
0260 } else {
0261 auto temp = loadTFormulaFromDB(es,
0262 mvaOutputNormalizationName_,
0263 formulaToken_,
0264 Form("%s_mvaOutput_normalization", moduleLabel_.data()),
0265 verbosity_);
0266 mvaOutput_normalization_ = std::move(temp);
0267 }
0268 }
0269 for (auto const& cutWPs : cuts_) {
0270 for (auto const& cut : cutWPs.second) {
0271 if (cut->mode_ == DiscriminantCutEntry::kVariableCut) {
0272 if (!loadMVAfromDB_) {
0273 if (not inputFile) {
0274 inputFile = openInputFile(inputFileName_);
0275 }
0276 if (verbosity_)
0277 std::cout << "Loading from file" << inputFileName_ << std::endl;
0278 cut->cutFunction_ = loadObjectFromFile<TGraph>(*inputFile, cut->cutName_);
0279 } else {
0280 if (verbosity_)
0281 std::cout << "Loading from DB" << std::endl;
0282 cut->cutFunction_ = loadTGraphFromDB(es, cut->cutName_, cut->cutToken_, verbosity_);
0283 }
0284 }
0285 }
0286 }
0287 isInitialized_ = true;
0288 }
0289
0290 evt.getByToken(toMultiplex_token, toMultiplexHandle_);
0291 }
0292
0293 template <class TauType, class TauTypeRef, class ParentClass>
0294 reco::SingleTauDiscriminatorContainer TauDiscriminantCutMultiplexerT<TauType, TauTypeRef, ParentClass>::discriminate(
0295 const TauTypeRef& tau) const {
0296 if (verbosity_) {
0297 std::cout << "<TauDiscriminantCutMultiplexerT::discriminate>:" << std::endl;
0298 std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
0299 }
0300
0301 reco::SingleTauDiscriminatorContainer result;
0302 result.rawValues.resize(n_raws_);
0303 double disc_result = (*toMultiplexHandle_)[tau].rawValues.at(0);
0304 if (verbosity_) {
0305 std::cout << "disc_result = " << disc_result << std::endl;
0306 }
0307 if (raw_discriminator_idx_ >= 0)
0308 result.rawValues[raw_discriminator_idx_] = disc_result;
0309 if (mvaOutput_normalization_) {
0310 disc_result = mvaOutput_normalization_->Eval(disc_result);
0311 if (verbosity_) {
0312 std::cout << "disc_result (normalized) = " << disc_result << std::endl;
0313 }
0314 }
0315 double key_result = 0.0;
0316 if ((*toMultiplexHandle_)[tau].rawValues.size() == 2) {
0317 key_result = (*toMultiplexHandle_)[tau].rawValues.at(1);
0318 if (raw_category_idx_ >= 0)
0319 result.rawValues[raw_category_idx_] = key_result;
0320 }
0321 typename DiscriminantCutMap::const_iterator cutWPsIter = cuts_.find(std::round(key_result));
0322
0323
0324 if (cutWPsIter == cuts_.end()) {
0325 return result;
0326 }
0327
0328 for (auto const& cutIter : cutWPsIter->second) {
0329 bool passesCuts = false;
0330 if (cutIter->mode_ == DiscriminantCutEntry::kFixedCut) {
0331 passesCuts = (disc_result > cutIter->cutValue_);
0332 if (verbosity_) {
0333 std::cout << "cutValue (fixed) = " << cutIter->cutValue_ << " --> passesCuts = " << passesCuts << std::endl;
0334 }
0335 } else if (cutIter->mode_ == DiscriminantCutEntry::kVariableCut) {
0336 double cutVariable = (*cutIter->cutVariable_)(*tau);
0337 double xMin, xMax, dummy;
0338 cutIter->cutFunction_->GetPoint(0, xMin, dummy);
0339 cutIter->cutFunction_->GetPoint(cutIter->cutFunction_->GetN() - 1, xMax, dummy);
0340 const double epsilon = 1.e-3;
0341 if (cutVariable < (xMin + epsilon))
0342 cutVariable = xMin + epsilon;
0343 else if (cutVariable > (xMax - epsilon))
0344 cutVariable = xMax - epsilon;
0345 double cutValue = cutIter->cutFunction_->Eval(cutVariable);
0346 passesCuts = (disc_result > cutValue);
0347 if (verbosity_) {
0348 std::cout << "cutValue (@" << cutVariable << ") = " << cutValue << " --> passesCuts = " << passesCuts
0349 << std::endl;
0350 }
0351 }
0352 result.workingPoints.push_back(passesCuts);
0353 }
0354 return result;
0355 }
0356
0357
0358 template <class TauType>
0359 std::string getDefaultConfigString() {
0360
0361
0362 throw cms::Exception("TauDiscriminantCutMultiplexerT")
0363 << "Unsupported TauType used. You must use either PFTau or PATTau.";
0364 }
0365
0366 template <>
0367 std::string getDefaultConfigString<reco::PFTau>() {
0368 return "recoTauDiscriminantCutMultiplexerDefault";
0369 }
0370 template <>
0371 std::string getDefaultConfigString<pat::Tau>() {
0372 return "patTauDiscriminantCutMultiplexerDefault";
0373 }
0374
0375 template <class TauType, class TauTypeRef, class ParentClass>
0376 void TauDiscriminantCutMultiplexerT<TauType, TauTypeRef, ParentClass>::fillDescriptions(
0377 edm::ConfigurationDescriptions& descriptions) {
0378
0379 edm::ParameterSetDescription desc;
0380 desc.add<edm::InputTag>("toMultiplex", edm::InputTag("fixme"));
0381 desc.add<int>("verbosity", 0);
0382
0383 {
0384 edm::ParameterSet pset_mapping;
0385 pset_mapping.addParameter<unsigned int>("category", 0);
0386 pset_mapping.addParameter<std::string>("cut", "fixme");
0387 edm::ParameterSetDescription desc_mapping;
0388 desc_mapping.add<unsigned int>("category", 0);
0389 desc_mapping.add<std::string>("cut");
0390
0391
0392 desc_mapping.addOptional<std::string>("variable")
0393 ->setComment("the parameter is required when \"workingPoints\" are string");
0394 std::vector<edm::ParameterSet> vpsd_mapping;
0395 vpsd_mapping.push_back(pset_mapping);
0396 desc.addVPSet("mapping", desc_mapping, vpsd_mapping);
0397 }
0398
0399 std::vector<std::string> defaultRaws{"discriminator"};
0400 desc.add<std::vector<std::string>>("rawValues", defaultRaws);
0401 std::vector<double> defaultWP{0.0};
0402 desc.addNode(edm::ParameterDescription<std::vector<double>>("workingPoints", defaultWP, true) xor
0403 edm::ParameterDescription<std::vector<std::string>>("workingPoints", true));
0404 desc.add<edm::FileInPath>("inputFileName", edm::FileInPath("RecoTauTag/RecoTau/data/emptyMVAinputFile"));
0405 desc.add<bool>("loadMVAfromDB", true);
0406 ParentClass::fillProducerDescriptions(desc);
0407 desc.add<std::string>("mvaOutput_normalization", "");
0408 descriptions.add(getDefaultConfigString<TauType>(), desc);
0409 }
0410
0411
0412 typedef TauDiscriminantCutMultiplexerT<reco::PFTau, reco::PFTauRef, PFTauDiscriminationContainerProducerBase>
0413 RecoTauDiscriminantCutMultiplexer;
0414 typedef TauDiscriminantCutMultiplexerT<pat::Tau, pat::TauRef, PATTauDiscriminationContainerProducerBase>
0415 PATTauDiscriminantCutMultiplexer;
0416
0417 DEFINE_FWK_MODULE(RecoTauDiscriminantCutMultiplexer);
0418 DEFINE_FWK_MODULE(PATTauDiscriminantCutMultiplexer);