GenVisTauProducer.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "DataFormats/Common/interface/View.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "DataFormats/JetReco/interface/GenJet.h"
#include "DataFormats/JetReco/interface/GenJetCollection.h"
#include "PhysicsTools/JetMCUtils/interface/JetMCTag.h"
#include "DataFormats/TauReco/interface/PFTau.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "CommonTools/Utils/interface/PtComparator.h"

#include <vector>
#include <iostream>

class GenVisTauProducer : public edm::global::EDProducer<> 
{
 public:
  GenVisTauProducer(const edm::ParameterSet& params) 
    : src_(consumes<reco::GenJetCollection>(params.getParameter<edm::InputTag>("src")))
    , srcGenParticles_(consumes<reco::GenParticleCollection>(params.getParameter<edm::InputTag>("srcGenParticles")))
    , pTComparator_()
  {
    produces<reco::GenParticleCollection>();
  }
  
  ~GenVisTauProducer() override {}

  void produce(edm::StreamID id, edm::Event& evt, const edm::EventSetup& es) const override 
  {
    edm::Handle<reco::GenJetCollection> genTauJets;
    evt.getByToken(src_, genTauJets);

    edm::Handle<reco::GenParticleCollection> genParticles;
    evt.getByToken(srcGenParticles_, genParticles);
    size_t numGenParticles = genParticles->size();

    auto genVisTaus = std::make_unique<reco::GenParticleCollection>();

    for (const auto & genTauJet : *genTauJets) {
      std::string decayMode_string = JetMCTagUtils::genTauDecayMode(genTauJet);
      // CV: store hadronic tau decays only
      if ( decayMode_string == "electron" || decayMode_string == "muon" ) continue;
      int decayMode = reco::PFTau::kNull;
      if      ( decayMode_string == "oneProng0Pi0"   ) decayMode = reco::PFTau::kOneProng0PiZero;
      else if ( decayMode_string == "oneProng1Pi0"   ) decayMode = reco::PFTau::kOneProng1PiZero;
      else if ( decayMode_string == "oneProng2Pi0"   ) decayMode = reco::PFTau::kOneProng2PiZero;
      else if ( decayMode_string == "threeProng0Pi0" ) decayMode = reco::PFTau::kThreeProng0PiZero;
      else if ( decayMode_string == "threeProng1Pi0" ) decayMode = reco::PFTau::kThreeProng1PiZero;
      else                                             decayMode = reco::PFTau::kRareDecayMode;

      int pdgId = ( genTauJet.charge() > 0 ) ? -15 : +15;

      // CV: store decayMode in status flag of GenParticle object
      reco::GenParticle genVisTau(genTauJet.charge(), genTauJet.p4(), genTauJet.vertex(), pdgId, decayMode, true);

      // CV: find tau lepton "mother" particle
      for ( size_t idxGenParticle = 0; idxGenParticle < numGenParticles; ++idxGenParticle ) {
    const reco::GenParticle & genTau = (*genParticles)[idxGenParticle];
    if ( abs(genTau.pdgId()) == 15 && genTau.status() == 2 ) {
      reco::Candidate::LorentzVector daughterVisP4;
          for (const reco::GenParticleRef & daughter : genTau.daughterRefVector()) {
        int abs_pdgId = abs(daughter->pdgId());
        // CV: skip neutrinos
        if ( abs_pdgId == 12 || abs_pdgId == 14 || abs_pdgId == 16 ) continue;
        daughterVisP4 += daughter->p4();
      }
      double dR2 = deltaR2(daughterVisP4, genVisTau);
      if ( dR2 < 1.e-4 ) {    
        genVisTau.addMother(reco::GenParticleRef(genParticles, idxGenParticle));
        break;
      }
    }
      }

      genVisTaus->push_back(genVisTau);
    }

    std::sort(genVisTaus->begin(), genVisTaus->end(), pTComparator_);
    evt.put(std::move(genVisTaus));
  }

  static void fillDescriptions(edm::ConfigurationDescriptions & descriptions) {
      edm::ParameterSetDescription desc;
      desc.add<edm::InputTag>("src")->setComment("collection of visible gen taus (as reco::GenJetCollection)");
      desc.add<edm::InputTag>("srcGenParticles")->setComment("collections of gen particles");
      descriptions.add("genVisTaus", desc);
  }


 private:
  const edm::EDGetTokenT<reco::GenJetCollection> src_;
  const edm::EDGetTokenT<reco::GenParticleCollection> srcGenParticles_;
  const GreaterByPt<reco::GenParticle> pTComparator_;
};

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