L1TkEmParticleProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
//
// dummy producer for a TkEm
//
 
// system include files
#include <memory>
 
// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
 
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
 
#include "DataFormats/Common/interface/Handle.h"
 
#include "DataFormats/L1TCorrelator/interface/TkEm.h"
#include "DataFormats/L1TCorrelator/interface/TkEmFwd.h"
#include "DataFormats/L1TCorrelator/interface/TkPrimaryVertex.h"
 
#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/Math/interface/deltaR.h"
 
// for L1Tracks:
#include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
 
#include <string>
#include <cmath>
 
static constexpr float EtaECal = 1.479;
static constexpr float REcal = 129.;
static constexpr float ZEcal = 315.4;
using namespace l1t;
//
// class declaration
//
 
class L1TkEmParticleProducer : public edm::global::EDProducer<> {
public:
  typedef TTTrack<Ref_Phase2TrackerDigi_> L1TTTrackType;
  typedef std::vector<L1TTTrackType> L1TTTrackCollectionType;
 
  explicit L1TkEmParticleProducer(const edm::ParameterSet&);
  ~L1TkEmParticleProducer() override;
 
  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
 
  float CorrectedEta(float eta, float zv) const;
 
private:
  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
 
  // ----------member data ---------------------------
 
  const edm::EDGetTokenT<EGammaBxCollection> egToken_;
  const edm::EDGetTokenT<L1TTTrackCollectionType> trackToken_;
  const edm::EDGetTokenT<TkPrimaryVertexCollection> vertexToken_;
 
  const std::string label_;
 
  const float etMin_;  // min ET in GeV of L1EG objects
  const float zMax_;   // |z_track| < zMax_ in cm
  const float chi2Max_;
  const float pTMinTra_;
  const float dR2Min_;
  const float dR2Max_;
  const bool primaryVtxConstrain_;  // use the primary vertex (default = false)
  const float deltaZMax_;           // | z_track - z_primaryvtx | < deltaZMax_ in cm.
                                    // Used only when primaryVtxConstrain_ = True.
  const float isoCut_;
  const bool relativeIsolation_;
};
 
//
// constructors and destructor
//
L1TkEmParticleProducer::L1TkEmParticleProducer(const edm::ParameterSet& iConfig)
    : egToken_(consumes<EGammaBxCollection>(iConfig.getParameter<edm::InputTag>("L1EGammaInputTag"))),
      trackToken_(consumes<L1TTTrackCollectionType>(iConfig.getParameter<edm::InputTag>("L1TrackInputTag"))),
      vertexToken_(consumes<TkPrimaryVertexCollection>(iConfig.getParameter<edm::InputTag>("L1VertexInputTag"))),
      label_(iConfig.getParameter<std::string>("label")),
      etMin_((float)iConfig.getParameter<double>("ETmin")),
      zMax_((float)iConfig.getParameter<double>("ZMAX")),
      chi2Max_((float)iConfig.getParameter<double>("CHI2MAX")),
      pTMinTra_((float)iConfig.getParameter<double>("PTMINTRA")),
      dR2Min_((float)iConfig.getParameter<double>("DRmin") * (float)iConfig.getParameter<double>("DRmin")),
      dR2Max_((float)iConfig.getParameter<double>("DRmax") * (float)iConfig.getParameter<double>("DRmax")),
      primaryVtxConstrain_(iConfig.getParameter<bool>("PrimaryVtxConstrain")),
      deltaZMax_((float)iConfig.getParameter<double>("DeltaZMax")),
      isoCut_((float)iConfig.getParameter<double>("IsoCut")),
      relativeIsolation_(iConfig.getParameter<bool>("RelativeIsolation")) {
  produces<TkEmCollection>(label_);
}
 
L1TkEmParticleProducer::~L1TkEmParticleProducer() {}
 
// ------------ method called to produce the data  ------------
void L1TkEmParticleProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
  using namespace edm;
 
  auto result = std::make_unique<TkEmCollection>();
 
  // the L1Tracks
  edm::Handle<L1TTTrackCollectionType> L1TTTrackHandle;
  iEvent.getByToken(trackToken_, L1TTTrackHandle);
  if (!L1TTTrackHandle.isValid()) {
    LogError("L1TkEmParticleProducer") << "\nWarning: L1TTTrackCollectionType not found in the event. Exit."
                                       << std::endl;
    return;
  }
 
  // the primary vertex (used only if primaryVtxConstrain_ = true)
  float zvtxL1tk = -999;
  bool primaryVtxConstrain = primaryVtxConstrain_;
  edm::Handle<TkPrimaryVertexCollection> L1VertexHandle;
  iEvent.getByToken(vertexToken_, L1VertexHandle);
  if (!L1VertexHandle.isValid()) {
    LogWarning("L1TkEmParticleProducer")
        << "Warning: TkPrimaryVertexCollection not found in the event. Won't use any PrimaryVertex constraint."
        << std::endl;
    primaryVtxConstrain = false;
  } else {
    std::vector<TkPrimaryVertex>::const_iterator vtxIter = L1VertexHandle->begin();
    // by convention, the first vertex in the collection is the one that should
    // be used by default
    zvtxL1tk = vtxIter->zvertex();
  }
 
  // Now loop over the L1EGamma objects
 
  edm::Handle<EGammaBxCollection> eGammaHandle;
  iEvent.getByToken(egToken_, eGammaHandle);
  if (!eGammaHandle.isValid()) {
    LogError("L1TkEmParticleProducer") << "\nWarning: L1EmCollection not found in the event. Exit." << std::endl;
    return;
  }
  EGammaBxCollection eGammaCollection = (*eGammaHandle.product());
  EGammaBxCollection::const_iterator egIter;
 
  int ieg = 0;
  for (egIter = eGammaCollection.begin(0); egIter != eGammaCollection.end(0); ++egIter)  // considering BX = only
  {
    edm::Ref<EGammaBxCollection> EGammaRef(eGammaHandle, ieg);
    ieg++;
 
    float et = egIter->et();
    if (et < etMin_)
      continue;
 
    float eta = egIter->eta();
    // The eta of the L1EG object is seen from (0,0,0).
    // if primaryVtxConstrain_ = true, and for the PV constrained iso, use the zvtxL1tk to correct the eta(L1EG)
    // that is used in the calculation of DeltaR.
    float etaPV = CorrectedEta(eta, zvtxL1tk);
 
    float phi = egIter->phi();
 
    // calculate the isolation of the L1EG object with
    // respect to L1Tracks.
 
    float sumPt = 0;
    float sumPtPV = 0;
 
    for (const auto& track : *L1TTTrackHandle) {
      float Pt = track.momentum().perp();
      float Eta = track.momentum().eta();
      float Phi = track.momentum().phi();
      float z = track.POCA().z();
      if (fabs(z) > zMax_)
        continue;
      if (Pt < pTMinTra_)
        continue;
      float chi2 = track.chi2();
      if (chi2 > chi2Max_)
        continue;
 
      float dr2 = reco::deltaR2(Eta, Phi, eta, phi);
      if (dr2 < dR2Max_ && dr2 >= dR2Min_) {
        sumPt += Pt;
      }
 
      if (zvtxL1tk > -999 && std::abs(z - zvtxL1tk) >= deltaZMax_)
        continue;  // Now, PV constrained trackSum:
 
      dr2 = reco::deltaR2(Eta, Phi, etaPV, phi);  // recompute using the corrected eta
      if (dr2 < dR2Max_ && dr2 >= dR2Min_) {
        sumPtPV += Pt;
      }
    }  // end loop over tracks
 
    float trkisol = -999;
    float trkisolPV = -999;
    if (relativeIsolation_) {
      if (et > 0) {
        trkisol = sumPt / et;
        trkisolPV = sumPtPV / et;
      }       // relative isolation
    } else {  // absolute isolation
      trkisol = sumPt;
      trkisolPV = sumPtPV;
    }
 
    float isolation = primaryVtxConstrain ? trkisolPV : trkisol;
 
    const math::XYZTLorentzVector P4 = egIter->p4();
    TkEm trkEm(P4, EGammaRef, trkisol, trkisolPV);
 
    if (isoCut_ <= 0) {
      // write the L1TkEm particle to the collection,
      // irrespective of its relative isolation
      result->push_back(trkEm);
    } else {
      // the object is written to the collection only
      // if it passes the isolation cut
      if (isolation <= isoCut_)
        result->push_back(trkEm);
    }
 
  }  // end loop over EGamma objects
 
  iEvent.put(std::move(result), label_);
}
 
// --------------------------------------------------------------------------------------
 
float L1TkEmParticleProducer::CorrectedEta(float eta, float zv) const {
  // Correct the eta of the L1EG object once we know the zvertex
 
  if (zv == 0.)
    return eta;
 
  if (eta == 0) {
    float thetaprime = atan(-REcal / zv);
    if (thetaprime < 0)
      thetaprime = thetaprime + M_PI;
    float etaprime = -log(tan(0.5 * thetaprime));
    return etaprime;
  }
 
  bool IsBarrel = (std::abs(eta) < EtaECal);
 
  float tanhalftheta = exp(-eta);
  float tantheta = 2. * tanhalftheta / (1. - tanhalftheta * tanhalftheta);
 
  float delta;
  if (IsBarrel)
    delta = REcal / tantheta;
  else
    delta = eta > 0 ? ZEcal : -ZEcal;
 
  float etaprime;
  if (delta == zv) {
    etaprime = 0.;
  } else {
    float tanthetaprime = delta * tantheta / (delta - zv);
    float thetaprime = atan(tanthetaprime);
    if (thetaprime < 0)
      thetaprime = thetaprime + M_PI;
    etaprime = -log(tan(0.5 * thetaprime));
  }
 
  return etaprime;
}
 
// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void L1TkEmParticleProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  //The following says we do not know what parameters are allowed so do no validation
  // Please change this to state exactly what you do use, even if it is no parameters
  edm::ParameterSetDescription desc;
  desc.add<string>("label", "EG");
  desc.add<edm::InputTag>("L1EGammaInputTag", edm::InputTag("simCaloStage2Digis"));
  desc.add<edm::InputTag>("L1TrackInputTag", edm::InputTag("TTTracksFromTrackletEmulation", "Level1TTTracks"));
  desc.add<edm::InputTag>("L1VertexInputTag", edm::InputTag("L1TkPrimaryVertex"));
  desc.add<double>("ETmin", -1.);
  desc.add<bool>("RelativeIsolation", true);
  desc.add<double>("IsoCut", 0.23);
  desc.add<double>("ZMAX", 25.);
  desc.add<double>("CHI2MAX", 100.);
  desc.add<double>("PTMINTRA", 2.);
  desc.add<double>("DRmin", 0.07);
  desc.add<double>("DRmax", 0.30);
  desc.add<bool>("PrimaryVtxConstrain", false);
  desc.add<double>("DeltaZMax", 0.6);
  descriptions.add("l1TkEmParticleProducer", desc);
}
 
//define this as a plug-in
DEFINE_FWK_MODULE(L1TkEmParticleProducer);