TtFullHadKinFitProducer.cc

Go to the documentation of this file.

#include "TopQuarkAnalysis/TopKinFitter/plugins/TtFullHadKinFitProducer.h"

static const unsigned int nPartons = 6;

TtFullHadKinFitProducer::TtFullHadKinFitProducer(const edm::ParameterSet& cfg)
    : jetsToken_(consumes<std::vector<pat::Jet> >(cfg.getParameter<edm::InputTag>("jets"))),
      matchToken_(mayConsume<std::vector<std::vector<int> > >(cfg.getParameter<edm::InputTag>("match"))),
      useOnlyMatch_(cfg.getParameter<bool>("useOnlyMatch")),
      bTagAlgo_(cfg.getParameter<std::string>("bTagAlgo")),
      minBTagValueBJet_(cfg.getParameter<double>("minBTagValueBJet")),
      maxBTagValueNonBJet_(cfg.getParameter<double>("maxBTagValueNonBJet")),
      useBTagging_(cfg.getParameter<bool>("useBTagging")),
      bTags_(cfg.getParameter<unsigned int>("bTags")),
      jetCorrectionLevel_(cfg.getParameter<std::string>("jetCorrectionLevel")),
      maxNJets_(cfg.getParameter<int>("maxNJets")),
      maxNComb_(cfg.getParameter<int>("maxNComb")),
      maxNrIter_(cfg.getParameter<unsigned int>("maxNrIter")),
      maxDeltaS_(cfg.getParameter<double>("maxDeltaS")),
      maxF_(cfg.getParameter<double>("maxF")),
      jetParam_(cfg.getParameter<unsigned>("jetParametrisation")),
      constraints_(cfg.getParameter<std::vector<unsigned> >("constraints")),
      mW_(cfg.getParameter<double>("mW")),
      mTop_(cfg.getParameter<double>("mTop")),
      jetEnergyResolutionScaleFactors_(cfg.getParameter<std::vector<double> >("jetEnergyResolutionScaleFactors")),
      jetEnergyResolutionEtaBinning_(cfg.getParameter<std::vector<double> >("jetEnergyResolutionEtaBinning")) {
  if (cfg.exists("udscResolutions") && cfg.exists("bResolutions")) {
    udscResolutions_ = cfg.getParameter<std::vector<edm::ParameterSet> >("udscResolutions");
    bResolutions_ = cfg.getParameter<std::vector<edm::ParameterSet> >("bResolutions");
  } else if (cfg.exists("udscResolutions") || cfg.exists("bResolutions")) {
    if (cfg.exists("udscResolutions"))
      throw cms::Exception("Configuration")
          << "Parameter 'bResolutions' is needed if parameter 'udscResolutions' is defined!\n";
    else
      throw cms::Exception("Configuration")
          << "Parameter 'udscResolutions' is needed if parameter 'bResolutions' is defined!\n";
  }

  // define kinematic fit interface
  kinFitter = new TtFullHadKinFitter::KinFit(useBTagging_,
                                             bTags_,
                                             bTagAlgo_,
                                             minBTagValueBJet_,
                                             maxBTagValueNonBJet_,
                                             udscResolutions_,
                                             bResolutions_,
                                             jetEnergyResolutionScaleFactors_,
                                             jetEnergyResolutionEtaBinning_,
                                             jetCorrectionLevel_,
                                             maxNJets_,
                                             maxNComb_,
                                             maxNrIter_,
                                             maxDeltaS_,
                                             maxF_,
                                             jetParam_,
                                             constraints_,
                                             mW_,
                                             mTop_);

  // produces the following collections
  produces<std::vector<pat::Particle> >("PartonsB");
  produces<std::vector<pat::Particle> >("PartonsBBar");
  produces<std::vector<pat::Particle> >("PartonsLightQ");
  produces<std::vector<pat::Particle> >("PartonsLightQBar");
  produces<std::vector<pat::Particle> >("PartonsLightP");
  produces<std::vector<pat::Particle> >("PartonsLightPBar");

  produces<std::vector<std::vector<int> > >();
  produces<std::vector<double> >("Chi2");
  produces<std::vector<double> >("Prob");
  produces<std::vector<int> >("Status");
}

TtFullHadKinFitProducer::~TtFullHadKinFitProducer() { delete kinFitter; }

void TtFullHadKinFitProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
  // get jet collection
  edm::Handle<std::vector<pat::Jet> > jets;
  event.getByToken(jetsToken_, jets);

  // get match in case that useOnlyMatch_ is true
  std::vector<int> match;
  bool invalidMatch = false;
  if (useOnlyMatch_) {
    kinFitter->setUseOnlyMatch(true);
    // in case that only a ceratin match should be used, get match here
    edm::Handle<std::vector<std::vector<int> > > matches;
    event.getByToken(matchToken_, matches);
    match = *(matches->begin());
    // check if match is valid
    if (match.size() != nPartons) {
      invalidMatch = true;
    } else {
      for (unsigned int idx = 0; idx < match.size(); ++idx) {
        if (match[idx] < 0 || match[idx] >= (int)jets->size()) {
          invalidMatch = true;
          break;
        }
      }
    }
    kinFitter->setMatch(match);
  }

  kinFitter->setMatchInvalidity(invalidMatch);

  std::list<TtFullHadKinFitter::KinFitResult> fitResults = kinFitter->fit(*jets);

  // pointer for output collections
  std::unique_ptr<std::vector<pat::Particle> > pPartonsB(new std::vector<pat::Particle>);
  std::unique_ptr<std::vector<pat::Particle> > pPartonsBBar(new std::vector<pat::Particle>);
  std::unique_ptr<std::vector<pat::Particle> > pPartonsLightQ(new std::vector<pat::Particle>);
  std::unique_ptr<std::vector<pat::Particle> > pPartonsLightQBar(new std::vector<pat::Particle>);
  std::unique_ptr<std::vector<pat::Particle> > pPartonsLightP(new std::vector<pat::Particle>);
  std::unique_ptr<std::vector<pat::Particle> > pPartonsLightPBar(new std::vector<pat::Particle>);
  // pointer for meta information
  std::unique_ptr<std::vector<std::vector<int> > > pCombi(new std::vector<std::vector<int> >);
  std::unique_ptr<std::vector<double> > pChi2(new std::vector<double>);
  std::unique_ptr<std::vector<double> > pProb(new std::vector<double>);
  std::unique_ptr<std::vector<int> > pStatus(new std::vector<int>);

  unsigned int iComb = 0;
  for (std::list<TtFullHadKinFitter::KinFitResult>::const_iterator res = fitResults.begin(); res != fitResults.end();
       ++res) {
    if (maxNComb_ >= 1 && iComb == (unsigned int)maxNComb_) {
      break;
    }
    ++iComb;

    pPartonsB->push_back(res->B);
    pPartonsBBar->push_back(res->BBar);
    pPartonsLightQ->push_back(res->LightQ);
    pPartonsLightQBar->push_back(res->LightQBar);
    pPartonsLightP->push_back(res->LightP);
    pPartonsLightPBar->push_back(res->LightPBar);

    pCombi->push_back(res->JetCombi);
    pChi2->push_back(res->Chi2);
    pProb->push_back(res->Prob);
    pStatus->push_back(res->Status);
  }

  event.put(std::move(pCombi));
  event.put(std::move(pPartonsB), "PartonsB");
  event.put(std::move(pPartonsBBar), "PartonsBBar");
  event.put(std::move(pPartonsLightQ), "PartonsLightQ");
  event.put(std::move(pPartonsLightQBar), "PartonsLightQBar");
  event.put(std::move(pPartonsLightP), "PartonsLightP");
  event.put(std::move(pPartonsLightPBar), "PartonsLightPBar");
  event.put(std::move(pChi2), "Chi2");
  event.put(std::move(pProb), "Prob");
  event.put(std::move(pStatus), "Status");
}

#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(TtFullHadKinFitProducer);