EcalDetailedTimeRecHitProducer.cc

Go to the documentation of this file.

#include "RecoLocalCalo/EcalRecProducers/plugins/EcalDetailedTimeRecHitProducer.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalRecHitSimpleAlgo.h"
#include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
#include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
#include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
#include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
#include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
#include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h"

#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <iostream>
#include <cmath>
#include <vector>

#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"

#include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include "DataFormats/EcalDigi/interface/EcalTimeDigi.h"

#include "CLHEP/Units/GlobalPhysicalConstants.h"
#include "CLHEP/Units/GlobalSystemOfUnits.h"

EcalDetailedTimeRecHitProducer::EcalDetailedTimeRecHitProducer(const edm::ParameterSet& ps) : m_geometry(nullptr) {
  EBRecHitCollection_ = consumes<EBRecHitCollection>(ps.getParameter<edm::InputTag>("EBRecHitCollection"));
  EERecHitCollection_ = consumes<EERecHitCollection>(ps.getParameter<edm::InputTag>("EERecHitCollection"));

  ebTimeDigiCollection_ = consumes<EcalTimeDigiCollection>(ps.getParameter<edm::InputTag>("EBTimeDigiCollection"));
  eeTimeDigiCollection_ = consumes<EcalTimeDigiCollection>(ps.getParameter<edm::InputTag>("EETimeDigiCollection"));

  EBDetailedTimeRecHitCollection_ = ps.getParameter<std::string>("EBDetailedTimeRecHitCollection");
  EEDetailedTimeRecHitCollection_ = ps.getParameter<std::string>("EEDetailedTimeRecHitCollection");

  correctForVertexZPosition_ = ps.getParameter<bool>("correctForVertexZPosition");
  useMCTruthVertex_ = ps.getParameter<bool>("useMCTruthVertex");
  if (correctForVertexZPosition_) {
    if (not useMCTruthVertex_) {
      recoVertex_ = consumes<reco::VertexCollection>(ps.getParameter<edm::InputTag>("recoVertex"));
    } else {
      simVertex_ = consumes<edm::SimVertexContainer>(ps.getParameter<edm::InputTag>("simVertex"));
    }
  }

  ebTimeLayer_ = ps.getParameter<int>("EBTimeLayer");
  eeTimeLayer_ = ps.getParameter<int>("EETimeLayer");

  produces<EBRecHitCollection>(EBDetailedTimeRecHitCollection_);
  produces<EERecHitCollection>(EEDetailedTimeRecHitCollection_);
}

EcalDetailedTimeRecHitProducer::~EcalDetailedTimeRecHitProducer() {}

void EcalDetailedTimeRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
  using namespace edm;
  using namespace reco;

  edm::ESHandle<CaloGeometry> hGeometry;
  es.get<CaloGeometryRecord>().get(hGeometry);

  m_geometry = hGeometry.product();

  Handle<EBRecHitCollection> pEBRecHits;
  Handle<EERecHitCollection> pEERecHits;

  const EBRecHitCollection* EBRecHits = nullptr;
  const EERecHitCollection* EERecHits = nullptr;

  evt.getByToken(EBRecHitCollection_, pEBRecHits);
  if (pEBRecHits.isValid()) {
    EBRecHits = pEBRecHits.product();  // get a ptr to the product
#ifdef DEBUG
    LogDebug("EcalRecHitDebug") << "total # EB rechits to be re-calibrated: " << EBRecHits->size();
#endif
  }

  evt.getByToken(EERecHitCollection_, pEERecHits);
  if (pEERecHits.isValid()) {
    EERecHits = pEERecHits.product();  // get a ptr to the product
#ifdef DEBUG
    LogDebug("EcalRecHitDebug") << "total # EE uncalibrated rechits to be re-calibrated: " << EERecHits->size();
#endif
  }

  Handle<EcalTimeDigiCollection> pEBTimeDigis;
  Handle<EcalTimeDigiCollection> pEETimeDigis;

  const EcalTimeDigiCollection* ebTimeDigis = nullptr;
  const EcalTimeDigiCollection* eeTimeDigis = nullptr;

  evt.getByToken(ebTimeDigiCollection_, pEBTimeDigis);
  //evt.getByToken( digiProducer_, pEBTimeDigis);
  if (pEBTimeDigis.isValid()) {
    ebTimeDigis = pEBTimeDigis.product();  // get a ptr to the produc
    edm::LogInfo("EcalDetailedTimeRecHitInfo") << "total # ebTimeDigis: " << ebTimeDigis->size();
  }

  evt.getByToken(eeTimeDigiCollection_, pEETimeDigis);
  //evt.getByToken( digiProducer_, pEETimeDigis);
  if (pEETimeDigis.isValid()) {
    eeTimeDigis = pEETimeDigis.product();  // get a ptr to the product
    edm::LogInfo("EcalDetailedTimeRecHitInfo") << "total # eeTimeDigis: " << eeTimeDigis->size();
  }
  // collection of rechits to put in the event
  std::unique_ptr<EBRecHitCollection> EBDetailedTimeRecHits(new EBRecHitCollection);
  std::unique_ptr<EERecHitCollection> EEDetailedTimeRecHits(new EERecHitCollection);

  std::unique_ptr<GlobalPoint> vertex;

  if (correctForVertexZPosition_) {
    if (!useMCTruthVertex_) {
      //Get the first reco vertex
      // get primary vertices

      edm::Handle<VertexCollection> VertexHandle;
      evt.getByToken(recoVertex_, VertexHandle);

      if (VertexHandle.isValid()) {
        if (!(*VertexHandle).empty())  //at least 1 vertex
        {
          const reco::Vertex* myVertex = &(*VertexHandle)[0];
          vertex.reset(new GlobalPoint(myVertex->x(), myVertex->y(), myVertex->z()));
        }
      }

    } else {
      edm::Handle<SimVertexContainer> VertexHandle;
      evt.getByToken(simVertex_, VertexHandle);

      if (VertexHandle.isValid()) {
        if (!(*VertexHandle).empty())  //at least 1 vertex
        {
          assert((*VertexHandle)[0].vertexId() == 0);
          const SimVertex* myVertex = &(*VertexHandle)[0];
          vertex.reset(new GlobalPoint(myVertex->position().x(), myVertex->position().y(), myVertex->position().z()));
        }
      }
    }
  }

  if (EBRecHits && ebTimeDigis) {
    // loop over uncalibrated rechits to make calibrated ones
    for (EBRecHitCollection::const_iterator it = EBRecHits->begin(); it != EBRecHits->end(); ++it) {
      EcalRecHit aHit((*it));
      EcalTimeDigiCollection::const_iterator timeDigi = ebTimeDigis->find((*it).id());
      if (timeDigi != ebTimeDigis->end()) {
        if (timeDigi->sampleOfInterest() >= 0) {
          float myTime = (*timeDigi)[timeDigi->sampleOfInterest()];
          //Vertex corrected ToF
          if (vertex) {
            aHit.setTime(myTime + deltaTimeOfFlight(*vertex, (*it).id(), ebTimeLayer_));
          } else
            //Uncorrected ToF
            aHit.setTime(myTime);
        }
      }
      // leave standard time if no timeDigi is associated (e.g. noise recHits)
      EBDetailedTimeRecHits->push_back(aHit);
    }
  }

  if (EERecHits && eeTimeDigis) {
    // loop over uncalibrated rechits to make calibrated ones
    for (EERecHitCollection::const_iterator it = EERecHits->begin(); it != EERecHits->end(); ++it) {
      EcalRecHit aHit(*it);
      EcalTimeDigiCollection::const_iterator timeDigi = eeTimeDigis->find((*it).id());
      if (timeDigi != eeTimeDigis->end()) {
        if (timeDigi->sampleOfInterest() >= 0) {
          float myTime = (*timeDigi)[timeDigi->sampleOfInterest()];
          //Vertex corrected ToF
          if (vertex) {
            aHit.setTime(myTime + deltaTimeOfFlight(*vertex, (*it).id(), eeTimeLayer_));
          } else
            //Uncorrected ToF
            aHit.setTime(myTime);
        }
      }
      EEDetailedTimeRecHits->push_back(aHit);
    }
  }
  // put the collection of recunstructed hits in the event
  LogInfo("EcalDetailedTimeRecHitInfo") << "total # EB rechits: " << EBDetailedTimeRecHits->size();
  LogInfo("EcalDetailedTimeRecHitInfo") << "total # EE rechits: " << EEDetailedTimeRecHits->size();

  evt.put(std::move(EBDetailedTimeRecHits), EBDetailedTimeRecHitCollection_);
  evt.put(std::move(EEDetailedTimeRecHits), EEDetailedTimeRecHitCollection_);
}

double EcalDetailedTimeRecHitProducer::deltaTimeOfFlight(GlobalPoint& vertex, const DetId& detId, int layer) const {
  auto cellGeometry(m_geometry->getGeometry(detId));
  assert(nullptr != cellGeometry);
  GlobalPoint layerPos =
      cellGeometry->getPosition(double(layer) + 0.5);  //depth in mm in the middle of the layer position
  GlobalVector tofVector = layerPos - vertex;
  return (layerPos.mag() * cm - tofVector.mag() * cm) / (float)c_light;
}

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