PATTauDiscriminantCutMultiplexer.cc

Go to the documentation of this file.

/*
 * \class PATTauDiscriminantCutMultiplexer
 * 
 * Takes two PFTauDiscriminators.
 *
 * The "key" discriminantor is rounded to the nearest integer.
 *
 * A set of cuts for different keys on the "toMultiplex" discriminantor is
 * provided in the config file.
 *
 * Both the key and toMultiplex discriminators should map to the same PFTau
 * collection.
 * 
 * Adopted from RecoTauTag/RecoTau/plugins/RecoTauDiscriminantCutMultiplexer.cc
 * 
 * \author Alexander Nehrkorn, RWTH Aachen
 */

#include <boost/foreach.hpp>
#include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "CommonTools/Utils/interface/StringObjectFunction.h"
#include "DataFormats/PatCandidates/interface/Tau.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"

#include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
#include <FWCore/ParameterSet/interface/ParameterSetDescription.h>

#include "CondFormats/PhysicsToolsObjects/interface/PhysicsTGraphPayload.h"
#include "CondFormats/DataRecord/interface/PhysicsTGraphPayloadRcd.h"
#include "CondFormats/PhysicsToolsObjects/interface/PhysicsTFormulaPayload.h"
#include "CondFormats/DataRecord/interface/PhysicsTFormulaPayloadRcd.h"

#include "TMath.h"
#include "TGraph.h"
#include "TFormula.h"
#include "TFile.h"

class PATTauDiscriminantCutMultiplexer : public PATTauDiscriminationProducerBase {
public:
  explicit PATTauDiscriminantCutMultiplexer(const edm::ParameterSet& pset);

  ~PATTauDiscriminantCutMultiplexer() override;
  double discriminate(const pat::TauRef&) const override;
  void beginEvent(const edm::Event& event, const edm::EventSetup& eventSetup) override;

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
  std::string moduleLabel_;

  bool loadMVAfromDB_;
  edm::FileInPath inputFileName_;

  struct DiscriminantCutEntry {
    DiscriminantCutEntry() : cutVariable_(), cutFunction_(), mode_(kUndefined) {}
    ~DiscriminantCutEntry() {}
    double cutValue_;
    std::string cutName_;
    std::unique_ptr<StringObjectFunction<pat::Tau>> cutVariable_;
    std::unique_ptr<const TGraph> cutFunction_;
    enum { kUndefined, kFixedCut, kVariableCut };
    int mode_;
  };
  typedef std::map<int, std::unique_ptr<DiscriminantCutEntry>> DiscriminantCutMap;
  DiscriminantCutMap cuts_;

  std::string mvaOutputNormalizationName_;
  std::unique_ptr<const TFormula> mvaOutput_normalization_;

  bool isInitialized_;

  edm::InputTag toMultiplex_;
  edm::InputTag key_;
  edm::Handle<pat::PATTauDiscriminator> toMultiplexHandle_;
  edm::Handle<pat::PATTauDiscriminator> keyHandle_;
  edm::EDGetTokenT<pat::PATTauDiscriminator> toMultiplex_token;
  edm::EDGetTokenT<pat::PATTauDiscriminator> key_token;

  int verbosity_;
};

namespace {
  std::unique_ptr<TFile> openInputFile(const edm::FileInPath& inputFileName) {
    if (inputFileName.location() == edm::FileInPath::Unknown) {
      throw cms::Exception("PATTauDiscriminantCutMultiplexer::loadObjectFromFile")
          << " Failed to find File = " << inputFileName << " !!\n";
    }
    return std::unique_ptr<TFile>{new TFile(inputFileName.fullPath().data())};
  }

  template <typename T>
  std::unique_ptr<const T> loadObjectFromFile(TFile& inputFile, const std::string& objectName) {
    const T* object = dynamic_cast<T*>(inputFile.Get(objectName.data()));
    if (!object)
      throw cms::Exception("PATTauDiscriminantCutMultiplexer::loadObjectFromFile")
          << " Failed to load Object = " << objectName.data() << " from file = " << inputFile.GetName() << " !!\n";
    //Need to use TObject::Clone since the type T might be a base class
    return std::unique_ptr<const T>{static_cast<T*>(object->Clone())};
  }

  std::unique_ptr<const TGraph> loadTGraphFromDB(const edm::EventSetup& es,
                                                 const std::string& graphName,
                                                 const int& verbosity_ = 0) {
    if (verbosity_) {
      std::cout << "<loadTGraphFromDB>:" << std::endl;
      std::cout << " graphName = " << graphName << std::endl;
    }
    edm::ESHandle<PhysicsTGraphPayload> graphPayload;
    es.get<PhysicsTGraphPayloadRcd>().get(graphName, graphPayload);
    return std::unique_ptr<const TGraph>{new TGraph(*graphPayload.product())};
  }

  std::unique_ptr<TFormula> loadTFormulaFromDB(const edm::EventSetup& es,
                                               const std::string& formulaName,
                                               const TString& newName,
                                               const int& verbosity_ = 0) {
    if (verbosity_) {
      std::cout << "<loadTFormulaFromDB>:" << std::endl;
      std::cout << " formulaName = " << formulaName << std::endl;
    }
    edm::ESHandle<PhysicsTFormulaPayload> formulaPayload;
    es.get<PhysicsTFormulaPayloadRcd>().get(formulaName, formulaPayload);

    if (formulaPayload->formulas().size() == 1 && formulaPayload->limits().size() == 1) {
      return std::unique_ptr<TFormula>{new TFormula(newName, formulaPayload->formulas().at(0).data())};
    } else {
      throw cms::Exception("PATTauDiscriminantCutMultiplexer::loadTFormulaFromDB")
          << "Failed to load TFormula = " << formulaName << " from Database !!\n";
    }
    return std::unique_ptr<TFormula>{};
  }
}  // namespace

PATTauDiscriminantCutMultiplexer::PATTauDiscriminantCutMultiplexer(const edm::ParameterSet& cfg)
    : PATTauDiscriminationProducerBase(cfg),
      moduleLabel_(cfg.getParameter<std::string>("@module_label")),
      mvaOutput_normalization_(),
      isInitialized_(false) {
  toMultiplex_ = cfg.getParameter<edm::InputTag>("toMultiplex");
  toMultiplex_token = consumes<pat::PATTauDiscriminator>(toMultiplex_);
  key_ = cfg.getParameter<edm::InputTag>("key");
  key_token = consumes<pat::PATTauDiscriminator>(key_);

  verbosity_ = cfg.getParameter<int>("verbosity");
  loadMVAfromDB_ = cfg.getParameter<bool>("loadMVAfromDB");
  if (!loadMVAfromDB_) {
    inputFileName_ = cfg.getParameter<edm::FileInPath>("inputFileName");
  }
  if (verbosity_)
    std::cout << moduleLabel_ << " loadMVA = " << loadMVAfromDB_ << std::endl;
  mvaOutputNormalizationName_ = cfg.getParameter<std::string>(
      "mvaOutput_normalization");  // default value is "" which should just overwrite existing value with same empty string

  // Setup our cut map
  typedef std::vector<edm::ParameterSet> VPSet;
  VPSet mapping = cfg.getParameter<VPSet>("mapping");
  for (VPSet::const_iterator mappingEntry = mapping.begin(); mappingEntry != mapping.end(); ++mappingEntry) {
    unsigned category = mappingEntry->getParameter<uint32_t>("category");
    std::unique_ptr<DiscriminantCutEntry> cut{new DiscriminantCutEntry()};
    if (mappingEntry->existsAs<double>("cut")) {
      cut->cutValue_ = mappingEntry->getParameter<double>("cut");
      cut->mode_ = DiscriminantCutEntry::kFixedCut;
    } else if (mappingEntry->existsAs<std::string>("cut")) {
      cut->cutName_ = mappingEntry->getParameter<std::string>("cut");
      std::string cutVariable_string = mappingEntry->getParameter<std::string>("variable");
      cut->cutVariable_.reset(new StringObjectFunction<pat::Tau>(cutVariable_string));
      cut->mode_ = DiscriminantCutEntry::kVariableCut;
    } else {
      throw cms::Exception("PATTauDiscriminantCutMultiplexer") << " Undefined Configuration Parameter 'cut' !!\n";
    }
    cuts_[category] = std::move(cut);
  }
  verbosity_ = cfg.getParameter<int>("verbosity");
  if (verbosity_)
    std::cout << "constructed " << moduleLabel_ << std::endl;
}

PATTauDiscriminantCutMultiplexer::~PATTauDiscriminantCutMultiplexer() {}

void PATTauDiscriminantCutMultiplexer::beginEvent(const edm::Event& evt, const edm::EventSetup& es) {
  if (verbosity_)
    std::cout << " begin! " << moduleLabel_ << " " << isInitialized_ << std::endl;
  if (!isInitialized_) {
    //Only open the file once and we can close it when this routine is done
    // since all objects gotten from the file will have been copied
    std::unique_ptr<TFile> inputFile;
    if (!mvaOutputNormalizationName_.empty()) {
      if (!loadMVAfromDB_) {
        inputFile = openInputFile(inputFileName_);
        mvaOutput_normalization_ = loadObjectFromFile<TFormula>(*inputFile, mvaOutputNormalizationName_);
      } else {
        auto temp = loadTFormulaFromDB(
            es, mvaOutputNormalizationName_, Form("%s_mvaOutput_normalization", moduleLabel_.data()), verbosity_);
        mvaOutput_normalization_ = std::move(temp);
      }
    }
    for (DiscriminantCutMap::iterator cut = cuts_.begin(); cut != cuts_.end(); ++cut) {
      if (cut->second->mode_ == DiscriminantCutEntry::kVariableCut) {
        if (!loadMVAfromDB_) {
          if (not inputFile) {
            inputFile = openInputFile(inputFileName_);
          }
          if (verbosity_)
            std::cout << "Loading from file" << inputFileName_ << std::endl;
          cut->second->cutFunction_ = loadObjectFromFile<TGraph>(*inputFile, cut->second->cutName_);
        } else {
          if (verbosity_)
            std::cout << "Loading from DB" << std::endl;
          cut->second->cutFunction_ = loadTGraphFromDB(es, cut->second->cutName_, verbosity_);
        }
      }
    }
    isInitialized_ = true;
  }

  evt.getByToken(toMultiplex_token, toMultiplexHandle_);
  evt.getByToken(key_token, keyHandle_);
}

double PATTauDiscriminantCutMultiplexer::discriminate(const pat::TauRef& tau) const {
  if (verbosity_) {
    std::cout << "<PATTauDiscriminantCutMultiplexer::discriminate>:" << std::endl;
    std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
  }

  double disc_result = (*toMultiplexHandle_)[tau];
  if (verbosity_) {
    std::cout << "disc_result = " << disc_result << std::endl;
  }
  if (mvaOutput_normalization_) {
    disc_result = mvaOutput_normalization_->Eval(disc_result);
    //if ( disc_result > 1. ) disc_result = 1.;
    //if ( disc_result < 0. ) disc_result = 0.;
    if (verbosity_) {
      std::cout << "disc_result (normalized) = " << disc_result << std::endl;
    }
  }
  double key_result = (*keyHandle_)[tau];
  DiscriminantCutMap::const_iterator cutIter = cuts_.find(TMath::Nint(key_result));

  // Return null if it doesn't exist
  if (cutIter == cuts_.end()) {
    return prediscriminantFailValue_;
  }
  // See if the discriminator passes our cuts
  bool passesCuts = false;
  if (cutIter->second->mode_ == DiscriminantCutEntry::kFixedCut) {
    passesCuts = (disc_result > cutIter->second->cutValue_);
    if (verbosity_) {
      std::cout << "cutValue (fixed) = " << cutIter->second->cutValue_ << " --> passesCuts = " << passesCuts
                << std::endl;
    }
  } else if (cutIter->second->mode_ == DiscriminantCutEntry::kVariableCut) {
    double cutVariable = (*cutIter->second->cutVariable_)(*tau);
    double xMin, xMax, dummy;
    cutIter->second->cutFunction_->GetPoint(0, xMin, dummy);
    cutIter->second->cutFunction_->GetPoint(cutIter->second->cutFunction_->GetN() - 1, xMax, dummy);
    const double epsilon = 1.e-3;
    if (cutVariable < (xMin + epsilon))
      cutVariable = xMin + epsilon;
    else if (cutVariable > (xMax - epsilon))
      cutVariable = xMax - epsilon;
    double cutValue = cutIter->second->cutFunction_->Eval(cutVariable);
    passesCuts = (disc_result > cutValue);
    if (verbosity_) {
      std::cout << "cutValue (@" << cutVariable << ") = " << cutValue << " --> passesCuts = " << passesCuts
                << std::endl;
    }
  } else
    assert(0);

  return passesCuts;
}

void PATTauDiscriminantCutMultiplexer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  // patTauDiscriminantCutMultiplexer
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("toMultiplex", edm::InputTag("fixme"));
  desc.add<int>("verbosity", 0);
  {
    edm::ParameterSetDescription desc_mapping;
    desc_mapping.add<unsigned int>("category");
    desc_mapping.setAllowAnything();  //option cut can be double or (string plus additional option variable)

    std::vector<edm::ParameterSet> vpsd_mapping;
    edm::ParameterSet psd1;
    psd1.addParameter<unsigned int>("category", 0);
    psd1.addParameter<double>("cut", 0.5);
    vpsd_mapping.push_back(psd1);
    edm::ParameterSet psd2;
    psd2.addParameter<unsigned int>("category", 1);
    psd2.addParameter<double>("cut", 0.2);
    vpsd_mapping.push_back(psd2);
    desc.addVPSet("mapping", desc_mapping, vpsd_mapping);
  }
  desc.add<edm::FileInPath>("inputFileName", edm::FileInPath("RecoTauTag/RecoTau/data/emptyMVAinputFile"));
  desc.add<bool>("loadMVAfromDB", true);
  fillProducerDescriptions(desc);  // inherited from the base
  desc.add<std::string>("mvaOutput_normalization", "");
  desc.add<edm::InputTag>("key", edm::InputTag("fixme"));
  descriptions.add("patTauDiscriminantCutMultiplexer", desc);
}

DEFINE_FWK_MODULE(PATTauDiscriminantCutMultiplexer);