PFClusterProducer/interface/PFHFRecHitCreator.h

0001 #ifndef RecoParticleFlow_PFClusterProducer_PFHFRecHitCreator_h
0002 #define RecoParticleFlow_PFClusterProducer_PFHFRecHitCreator_h
0003
0004 #include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitCreatorBase.h"
0005
0006 #include "DataFormats/HcalRecHit/interface/HORecHit.h"
0007 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
0008 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0009 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0010 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0011
0012 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0013 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0014 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0015 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0016 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0017 #include "Geometry/CaloGeometry/interface/IdealZPrism.h"
0018
0019 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0020
0021 class PFHFRecHitCreator final : public PFRecHitCreatorBase {
0022 public:
0023   PFHFRecHitCreator(const edm::ParameterSet& iConfig, edm::ConsumesCollector& cc)
0024       : PFRecHitCreatorBase(iConfig, cc),
0025         recHitToken_(cc.consumes<edm::SortedCollection<HFRecHit> >(iConfig.getParameter<edm::InputTag>("src"))),
0026         EM_Depth_(iConfig.getParameter<double>("EMDepthCorrection")),
0027         HAD_Depth_(iConfig.getParameter<double>("HADDepthCorrection")),
0028         shortFibre_Cut(iConfig.getParameter<double>("ShortFibre_Cut")),
0029         longFibre_Fraction(iConfig.getParameter<double>("LongFibre_Fraction")),
0030         longFibre_Cut(iConfig.getParameter<double>("LongFibre_Cut")),
0031         shortFibre_Fraction(iConfig.getParameter<double>("ShortFibre_Fraction")),
0032         thresh_HF_(iConfig.getParameter<double>("thresh_HF")),
0033         HFCalib_(iConfig.getParameter<double>("HFCalib29")),
0034         geomToken_(cc.esConsumes()) {}
0035
0036   void importRecHits(std::unique_ptr<reco::PFRecHitCollection>& out,
0037                      std::unique_ptr<reco::PFRecHitCollection>& cleaned,
0038                      const edm::Event& iEvent,
0039                      const edm::EventSetup& iSetup) override {
0040     reco::PFRecHitCollection tmpOut;
0041
0042     beginEvent(iEvent, iSetup);
0043
0044     edm::Handle<edm::SortedCollection<HFRecHit> > recHitHandle;
0045
0046     edm::ESHandle<CaloGeometry> geoHandle = iSetup.getHandle(geomToken_);
0047
0048     // get the ecal geometry
0049     const CaloSubdetectorGeometry* hcalGeo = geoHandle->getSubdetectorGeometry(DetId::Hcal, HcalForward);
0050
0051     iEvent.getByToken(recHitToken_, recHitHandle);
0052     for (const auto& erh : *recHitHandle) {
0053       const HcalDetId& detid = (HcalDetId)erh.detid();
0054       auto depth = detid.depth();
0055
0056       // ATTN: skip dual anode in HF for now (should be fixed in upstream changes)
0057       if (depth > 2)
0058         continue;
0059
0060       auto energy = erh.energy();
0061       auto time = erh.time();
0062
0063       std::shared_ptr<const CaloCellGeometry> thisCell = hcalGeo->getGeometry(detid);
0064       auto zp = dynamic_cast<IdealZPrism const*>(thisCell.get());
0065       assert(zp);
0066       thisCell = zp->forPF();
0067
0068       // find rechit geometry
0069       if (!thisCell) {
0070         edm::LogError("PFHFRecHitCreator")
0071             << "warning detid " << detid.rawId() << " not found in geometry" << std::endl;
0072         continue;
0073       }
0074
0075       PFLayer::Layer layer = depth == 1 ? PFLayer::HF_EM : PFLayer::HF_HAD;
0076
0077       reco::PFRecHit rh(thisCell, detid.rawId(), layer, energy);
0078       rh.setTime(time);
0079       rh.setDepth(depth);
0080
0081       bool rcleaned = false;
0082       bool keep = true;
0083
0084       //Apply Q tests
0085       for (const auto& qtest : qualityTests_) {
0086         if (!qtest->test(rh, erh, rcleaned)) {
0087           keep = false;
0088         }
0089       }
0090
0091       if (keep) {
0092         tmpOut.push_back(std::move(rh));
0093       } else if (rcleaned)
0094         cleaned->push_back(std::move(rh));
0095     }
0096     //Sort by DetID the collection
0097     DetIDSorter sorter;
0098     if (!tmpOut.empty())
0099       std::sort(tmpOut.begin(), tmpOut.end(), sorter);
0100
0101     /////////////////////HF DUAL READOUT/////////////////////////
0102
0103     for (auto& hit : tmpOut) {
0104       reco::PFRecHit newHit = hit;
0105       const HcalDetId& detid = (HcalDetId)hit.detId();
0106       if (detid.depth() == 1) {
0107         double lONG = hit.energy();
0108         //find the short hit
0109         HcalDetId shortID(HcalForward, detid.ieta(), detid.iphi(), 2);
0110         auto found_hit =
0111             std::lower_bound(tmpOut.begin(), tmpOut.end(), shortID, [](const reco::PFRecHit& a, HcalDetId b) {
0112               return a.detId() < b.rawId();
0113             });
0114         if (found_hit != tmpOut.end() && found_hit->detId() == shortID.rawId()) {
0115           double sHORT = found_hit->energy();
0116           //Ask for fraction
0117           double energy = lONG - sHORT;
0118
0119           if (abs(detid.ieta()) <= 32)
0120             energy *= HFCalib_;
0121           newHit.setEnergy(energy);
0122           if (!(lONG > longFibre_Cut && (sHORT / lONG < shortFibre_Fraction)))
0123             if (energy > thresh_HF_)
0124               out->push_back(newHit);
0125         } else {
0126           //make only long hit
0127           double energy = lONG;
0128           if (abs(detid.ieta()) <= 32)
0129             energy *= HFCalib_;
0130           newHit.setEnergy(energy);
0131
0132           if (energy > thresh_HF_)
0133             out->push_back(newHit);
0134         }
0135
0136       } else {
0137         double sHORT = hit.energy();
0138         HcalDetId longID(HcalForward, detid.ieta(), detid.iphi(), 1);
0139         auto found_hit =
0140             std::lower_bound(tmpOut.begin(), tmpOut.end(), longID, [](const reco::PFRecHit& a, HcalDetId b) {
0141               return a.detId() < b.rawId();
0142             });
0143         double energy = sHORT;
0144         if (found_hit != tmpOut.end() && found_hit->detId() == longID.rawId()) {
0145           double lONG = found_hit->energy();
0146           //Ask for fraction
0147
0148           //If in this case lONG-sHORT<0 add the energy to the sHORT
0149           if ((lONG - sHORT) < thresh_HF_)
0150             energy += lONG;
0151           else
0152             energy += sHORT;
0153
0154           if (abs(detid.ieta()) <= 32)
0155             energy *= HFCalib_;
0156
0157           newHit.setEnergy(energy);
0158           if (!(sHORT > shortFibre_Cut && (lONG / sHORT < longFibre_Fraction)))
0159             if (energy > thresh_HF_)
0160               out->push_back(newHit);
0161
0162         } else {
0163           //only short hit!
0164           if (abs(detid.ieta()) <= 32)
0165             energy *= HFCalib_;
0166           newHit.setEnergy(energy);
0167           if (energy > thresh_HF_)
0168             out->push_back(newHit);
0169         }
0170       }
0171     }
0172   }
0173
0174 protected:
0175   edm::EDGetTokenT<edm::SortedCollection<HFRecHit> > recHitToken_;
0176   double EM_Depth_;
0177   double HAD_Depth_;
0178   // Don't allow large energy in short fibres if there is no energy in long fibres
0179   double shortFibre_Cut;
0180   double longFibre_Fraction;
0181
0182   // Don't allow large energy in long fibres if there is no energy in short fibres
0183   double longFibre_Cut;
0184   double shortFibre_Fraction;
0185   double thresh_HF_;
0186   double HFCalib_;
0187
0188   class DetIDSorter {
0189   public:
0190     DetIDSorter() = default;
0191     ~DetIDSorter() = default;
0192
0193     bool operator()(const reco::PFRecHit& a, const reco::PFRecHit& b) {
0194       if (DetId(a.detId()).det() == DetId::Hcal || DetId(b.detId()).det() == DetId::Hcal)
0195         return (HcalDetId)(a.detId()) < (HcalDetId)(b.detId());
0196       else
0197         return a.detId() < b.detId();
0198     }
0199   };
0200
0201 private:
0202   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0203 };
0204 #endif