PartonSelector.cc

Go to the documentation of this file.

//
// Select the partons status 2 and 3 for MC Jet Flavour
// Author: Attilio
// Date: 10.10.2007
//

//=======================================================================

// user include files
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/Ref.h"
//#include "DataFormats/Candidate/interface/CandidateFwd.h"
//#include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"

#include <memory>
#include <string>
#include <iostream>
#include <vector>

using namespace std;
using namespace reco;
using namespace edm;

class PartonSelector : public edm::global::EDProducer<> {
public:
  PartonSelector(const edm::ParameterSet&);
  ~PartonSelector() override;

private:
  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
  bool withLeptons;         // Optionally specify leptons
  bool withTop;             // Optionally include top quarks in the list
  bool acceptNoDaughters;   // Parton with zero daugthers are not considered by default, make it configurable
  unsigned int skipFirstN;  // Default skips first 6 particles, make it configurable
  edm::EDGetTokenT<reco::GenParticleCollection> tokenGenParticles_;  // input collection
};
//=========================================================================

PartonSelector::PartonSelector(const edm::ParameterSet& iConfig) {
  produces<reco::GenParticleRefVector>();
  withLeptons = iConfig.getParameter<bool>("withLeptons");
  tokenGenParticles_ = consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("src"));
  if (iConfig.exists("acceptNoDaughters")) {
    acceptNoDaughters = iConfig.getParameter<bool>("acceptNoDaughters");
  } else {
    acceptNoDaughters = false;
  }
  if (iConfig.exists("skipFirstN")) {
    skipFirstN = iConfig.getParameter<unsigned int>("skipFirstN");
  } else {
    skipFirstN = 6;
  }
  if (iConfig.exists("withTop")) {
    withTop = iConfig.getParameter<bool>("withTop");
  } else {
    withTop = false;
  }
}

//=========================================================================

PartonSelector::~PartonSelector() {}

// ------------ method called to produce the data  ------------

void PartonSelector::produce(StreamID, Event& iEvent, const EventSetup& iEs) const {
  //edm::Handle <reco::CandidateView> particles;
  edm::Handle<reco::GenParticleCollection> particles;
  iEvent.getByToken(tokenGenParticles_, particles);
  edm::LogVerbatim("PartonSelector") << "=== GenParticle size:" << particles->size();
  int nPart = 0;

  auto thePartons = std::make_unique<GenParticleRefVector>();

  for (size_t m = 0; m < particles->size(); m++) {
    // Don't take into account first 6 particles in generator list
    if (m < skipFirstN)
      continue;

    const GenParticle& aParticle = (*particles)[m];

    bool isAParton = false;
    bool isALepton = false;
    int flavour = abs(aParticle.pdgId());
    if (flavour == 1 || flavour == 2 || flavour == 3 || flavour == 4 || flavour == 5 || (flavour == 6 && withTop) ||
        flavour == 21)
      isAParton = true;
    if (flavour == 11 || flavour == 12 || flavour == 13 || flavour == 14 || flavour == 15 || flavour == 16)
      isALepton = true;

    //Add Partons status 3
    if (aParticle.status() == 3 && isAParton) {
      thePartons->push_back(GenParticleRef(particles, m));
      nPart++;
    }

    //Add Partons status 2
    int nparton_daughters = 0;
    if ((aParticle.numberOfDaughters() > 0 || acceptNoDaughters) && isAParton) {
      for (unsigned int i = 0; i < aParticle.numberOfDaughters(); i++) {
        int daughterFlavour = abs(aParticle.daughter(i)->pdgId());
        if ((daughterFlavour == 1 || daughterFlavour == 2 || daughterFlavour == 3 || daughterFlavour == 4 ||
             daughterFlavour == 5 || daughterFlavour == 6 || daughterFlavour == 21)) {
          nparton_daughters++;
        }
      }
      if (nparton_daughters == 0) {
        nPart++;
        thePartons->push_back(GenParticleRef(particles, m));
      }
    }

    //Add Leptons
    // Here you have to decide what to do with taus....
    // Now all leptons, including e and mu from leptonic tau decays, are added
    if (withLeptons && aParticle.status() == 3 && isALepton) {
      thePartons->push_back(GenParticleRef(particles, m));
      nPart++;
    }
  }

  edm::LogVerbatim("PartonSelector") << "=== GenParticle selected:" << nPart;
  iEvent.put(std::move(thePartons));
}

DEFINE_FWK_MODULE(PartonSelector);