PFRecHitCaloNavigatorWithTime.h

Go to the documentation of this file.

#ifndef RecoParticleFlow_PFClusterProducer_PFRecHitCaloNavigatorWithTime_h
#define RecoParticleFlow_PFClusterProducer_PFRecHitCaloNavigatorWithTime_h

#include <memory>

#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h"

#include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalDetId/interface/ESDetId.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"

#include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
#include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
#include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"

#include "Geometry/CaloTopology/interface/CaloTowerTopology.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"

#include "RecoParticleFlow/PFClusterProducer/interface/CaloRecHitResolutionProvider.h"

template <typename D, typename T, bool ownsTopo = true>
class PFRecHitCaloNavigatorWithTime : public PFRecHitNavigatorBase {
public:
  PFRecHitCaloNavigatorWithTime(const edm::ParameterSet& iConfig, edm::ConsumesCollector& cc) {
    sigmaCut2_ = pow(iConfig.getParameter<double>("sigmaCut"), 2);
    const edm::ParameterSet& timeResConf = iConfig.getParameterSet("timeResolutionCalc");
    _timeResolutionCalc = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
  }

  ~PFRecHitCaloNavigatorWithTime() override {
    if (!ownsTopo) {
      topology_.release();
    }
  }

  void associateNeighbours(reco::PFRecHit& hit,
                           std::unique_ptr<reco::PFRecHitCollection>& hits,
                           edm::RefProd<reco::PFRecHitCollection>& refProd) override {
    DetId detid(hit.detId());

    CaloNavigator<D> navigator(detid, topology_.get());

    DetId N(0);
    DetId E(0);
    DetId S(0);
    DetId W(0);
    DetId NW(0);
    DetId NE(0);
    DetId SW(0);
    DetId SE(0);

    N = navigator.north();
    associateNeighbour(N, hit, hits, refProd, 0, 1);

    if (N != DetId(0)) {
      NE = navigator.east();
    } else {
      navigator.home();
      E = navigator.east();
      NE = navigator.north();
    }
    associateNeighbour(NE, hit, hits, refProd, 1, 1);
    navigator.home();

    S = navigator.south();
    associateNeighbour(S, hit, hits, refProd, 0, -1);

    if (S != DetId(0)) {
      SW = navigator.west();
    } else {
      navigator.home();
      W = navigator.west();
      SW = navigator.south();
    }
    associateNeighbour(SW, hit, hits, refProd, -1, -1);
    navigator.home();

    E = navigator.east();
    associateNeighbour(E, hit, hits, refProd, 1, 0);

    if (E != DetId(0)) {
      SE = navigator.south();
    } else {
      navigator.home();
      S = navigator.south();
      SE = navigator.east();
    }
    associateNeighbour(SE, hit, hits, refProd, 1, -1);
    navigator.home();

    W = navigator.west();
    associateNeighbour(W, hit, hits, refProd, -1, 0);

    if (W != DetId(0)) {
      NW = navigator.north();
    } else {
      navigator.home();
      N = navigator.north();
      NW = navigator.west();
    }
    associateNeighbour(NW, hit, hits, refProd, -1, 1);
  }

protected:
  double sigmaCut2_;
  std::unique_ptr<const T> topology_;
  std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalc;

  void associateNeighbour(const DetId& id,
                          reco::PFRecHit& hit,
                          std::unique_ptr<reco::PFRecHitCollection>& hits,
                          edm::RefProd<reco::PFRecHitCollection>& refProd,
                          short eta,
                          short phi) {
    double sigma2 = 10000.0;

    auto found_hit = std::lower_bound(hits->begin(), hits->end(), id, [](const reco::PFRecHit& a, DetId b) {
      if (DetId(a.detId()).det() == DetId::Hcal || b.det() == DetId::Hcal)
        return (HcalDetId)(a.detId()) < (HcalDetId)(b);

      else
        return a.detId() < b.rawId();
    });

    if (found_hit != hits->end() &&
        ((found_hit->detId() == id.rawId()) ||
         ((id.det() == DetId::Hcal) && ((HcalDetId)(found_hit->detId()) == (HcalDetId)(id))))) {
      sigma2 = _timeResolutionCalc->timeResolution2(hit.energy()) +
               _timeResolutionCalc->timeResolution2(found_hit->energy());
      const double deltaTime = hit.time() - found_hit->time();
      if (deltaTime * deltaTime < sigmaCut2_ * sigma2) {
        hit.addNeighbour(eta, phi, 0, found_hit - hits->begin());
      }
    }
  }
};

#endif