HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.cc

0001 #include "RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.h"
0002
0003 #include <memory>
0004
0005 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0006 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0007 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0008 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h"
0013
0014 HGCalRecHitWorkerSimple::HGCalRecHitWorkerSimple(const edm::ParameterSet& ps, edm::ConsumesCollector iC)
0015     : HGCalRecHitWorkerBaseClass(ps, iC),
0016       caloGeomToken_(iC.esConsumes<CaloGeometry, CaloGeometryRecord>()),
0017       ee_geometry_token_(iC.esConsumes(edm::ESInputTag("", "HGCalEESensitive"))),
0018       hef_geometry_token_(iC.esConsumes(edm::ESInputTag("", "HGCalHESiliconSensitive"))),
0019       hfnose_geometry_token_(iC.esConsumes(edm::ESInputTag("", "HGCalHFNoseSensitive"))) {
0020   rechitMaker_ = std::make_unique<HGCalRecHitSimpleAlgo>();
0021   tools_ = std::make_unique<hgcal::RecHitTools>();
0022   constexpr float keV2GeV = 1e-6;
0023
0024   // HGCee constants
0025   hgcEE_keV2DIGI_ = ps.getParameter<double>("HGCEE_keV2DIGI");
0026   hgcEE_fCPerMIP_ = ps.getParameter<std::vector<double> >("HGCEE_fCPerMIP");
0027   hgcEE_isSiFE_ = ps.getParameter<bool>("HGCEE_isSiFE");
0028   hgceeUncalib2GeV_ = keV2GeV / hgcEE_keV2DIGI_;
0029
0030   // HGChef constants
0031   hgcHEF_keV2DIGI_ = ps.getParameter<double>("HGCHEF_keV2DIGI");
0032   hgcHEF_fCPerMIP_ = ps.getParameter<std::vector<double> >("HGCHEF_fCPerMIP");
0033   hgcHEF_isSiFE_ = ps.getParameter<bool>("HGCHEF_isSiFE");
0034   hgchefUncalib2GeV_ = keV2GeV / hgcHEF_keV2DIGI_;
0035
0036   // HGCheb constants
0037   hgcHEB_keV2DIGI_ = ps.getParameter<double>("HGCHEB_keV2DIGI");
0038   hgcHEB_isSiFE_ = ps.getParameter<bool>("HGCHEB_isSiFE");
0039   hgchebUncalib2GeV_ = keV2GeV / hgcHEB_keV2DIGI_;
0040
0041   // HGChfnose constants
0042   hgcHFNose_keV2DIGI_ = ps.getParameter<double>("HGCHFNose_keV2DIGI");
0043   hgcHFNose_fCPerMIP_ = ps.getParameter<std::vector<double> >("HGCHFNose_fCPerMIP");
0044   hgcHFNose_isSiFE_ = ps.getParameter<bool>("HGCHFNose_isSiFE");
0045   hgchfnoseUncalib2GeV_ = keV2GeV / hgcHFNose_keV2DIGI_;
0046
0047   // layer weights (from Valeri/Arabella)
0048   const auto& dweights = ps.getParameter<std::vector<double> >("layerWeights");
0049   for (auto weight : dweights) {
0050     weights_.push_back(weight);
0051   }
0052   const auto& weightnose = ps.getParameter<std::vector<double> >("layerNoseWeights");
0053   for (auto const& weight : weightnose)
0054     weightsNose_.emplace_back(weight);
0055
0056   rechitMaker_->setLayerWeights(weights_);
0057   rechitMaker_->setNoseLayerWeights(weightsNose_);
0058
0059   // residual correction for cell thickness
0060   // first for silicon
0061   const auto& rcorr = ps.getParameter<std::vector<double> >("thicknessCorrection");
0062   rcorr_.clear();
0063   rcorr_.push_back(1.f);
0064   for (auto corr : rcorr) {
0065     rcorr_.push_back(1.0 / corr);
0066   }
0067   // here for scintillator
0068   rcorrscint_ = 1.0 / ps.getParameter<double>("sciThicknessCorrection");
0069
0070   //This is for the index position in CE_H silicon thickness cases
0071   deltasi_index_regemfac_ = ps.getParameter<int>("deltasi_index_regemfac");
0072   const auto& rcorrnose = ps.getParameter<std::vector<double> >("thicknessNoseCorrection");
0073   rcorrNose_.clear();
0074   rcorrNose_.push_back(1.f);
0075   for (auto corr : rcorrnose) {
0076     rcorrNose_.push_back(1.0 / corr);
0077   }
0078
0079   hgcEE_noise_fC_ = ps.getParameter<edm::ParameterSet>("HGCEE_noise_fC").getParameter<std::vector<double> >("values");
0080   hgcEE_cce_ = ps.getParameter<edm::ParameterSet>("HGCEE_cce").getParameter<std::vector<double> >("values");
0081   hgcHEF_noise_fC_ = ps.getParameter<edm::ParameterSet>("HGCHEF_noise_fC").getParameter<std::vector<double> >("values");
0082   hgcHEF_cce_ = ps.getParameter<edm::ParameterSet>("HGCHEF_cce").getParameter<std::vector<double> >("values");
0083   hgcHEB_noise_MIP_ = ps.getParameter<edm::ParameterSet>("HGCHEB_noise_MIP").getParameter<double>("noise_MIP");
0084   hgcHFNose_noise_fC_ =
0085       ps.getParameter<edm::ParameterSet>("HGCHFNose_noise_fC").getParameter<std::vector<double> >("values");
0086   hgcHFNose_cce_ = ps.getParameter<edm::ParameterSet>("HGCHFNose_cce").getParameter<std::vector<double> >("values");
0087
0088   // don't produce rechit if detid is a ghost one
0089   rangeMatch_ = ps.getParameter<uint32_t>("rangeMatch");
0090   rangeMask_ = ps.getParameter<uint32_t>("rangeMask");
0091
0092   // error for recHit time
0093   timeEstimatorSi_ = hgcalsimclustertime::ComputeClusterTime(ps.getParameter<double>("minValSiPar"),
0094                                                              ps.getParameter<double>("maxValSiPar"),
0095                                                              ps.getParameter<double>("constSiPar"),
0096                                                              ps.getParameter<double>("noiseSiPar"));
0097 }
0098
0099 void HGCalRecHitWorkerSimple::set(const edm::EventSetup& es) {
0100   const CaloGeometry& geom = es.getData(caloGeomToken_);
0101   tools_->setGeometry(geom);
0102   rechitMaker_->set(geom);
0103   if (hgcEE_isSiFE_) {
0104     const HGCalGeometry& hgceeGeoHandle = es.getData(ee_geometry_token_);
0105     ddds_[0] = &(hgceeGeoHandle.topology().dddConstants());
0106   } else {
0107     ddds_[0] = nullptr;
0108   }
0109   if (hgcHEF_isSiFE_) {
0110     edm::ESHandle<HGCalGeometry> hgchefGeoHandle = es.getHandle(hef_geometry_token_);
0111     ddds_[1] = &(hgchefGeoHandle->topology().dddConstants());
0112   } else {
0113     ddds_[1] = nullptr;
0114   }
0115   ddds_[2] = nullptr;
0116   if (hgcHFNose_isSiFE_) {
0117     edm::ESHandle<HGCalGeometry> hgchfnoseGeoHandle = es.getHandle(hfnose_geometry_token_);
0118     ddds_[3] = &(hgchfnoseGeoHandle->topology().dddConstants());
0119   } else {
0120     ddds_[3] = nullptr;
0121   }
0122 }
0123
0124 void HGCalRecHitWorkerSimple::run(const edm::Event& evt,
0125                                   const HGCUncalibratedRecHitCollection& uncalibRHColl,
0126                                   HGCRecHitCollection& result) {
0127   for (const auto& uncalibRH : uncalibRHColl) {
0128     DetId detid = uncalibRH.id();
0129     // don't produce rechit if detid is a ghost one
0130
0131     if (detid.det() == DetId::Forward || detid.det() == DetId::Hcal) {
0132       if ((detid & rangeMask_) == rangeMatch_)
0133         continue;
0134     }
0135
0136     int thickness = -1;
0137     float sigmaNoiseGeV = 0.f;
0138     unsigned int layer = tools_->getLayerWithOffset(detid);
0139     float cce_correction = 1.0f;
0140     int idtype(0);
0141
0142     switch (detid.det()) {
0143       case DetId::HGCalEE:
0144         idtype = hgcee;
0145         thickness = 1 + HGCSiliconDetId(detid).type();
0146         break;
0147       case DetId::HGCalHSi:
0148         idtype = hgcfh;
0149         thickness = 1 + HGCSiliconDetId(detid).type();
0150         break;
0151       case DetId::HGCalHSc:
0152         idtype = hgcbh;
0153         break;
0154       default:
0155         switch (detid.subdetId()) {
0156           case HGCEE:
0157             idtype = hgcee;
0158             thickness = ddds_[detid.subdetId() - 3]->waferTypeL(HGCalDetId(detid).wafer());
0159             break;
0160           case HGCHEF:
0161             idtype = hgcfh;
0162             thickness = ddds_[detid.subdetId() - 3]->waferTypeL(HGCalDetId(detid).wafer());
0163             break;
0164           case HcalEndcap:
0165             [[fallthrough]];
0166           case HGCHEB:
0167             idtype = hgcbh;
0168             break;
0169           case HFNose:
0170             idtype = hgchfnose;
0171             thickness = 1 + HFNoseDetId(detid).type();
0172             break;
0173           default:
0174             break;
0175         }
0176     }
0177
0178     switch (idtype) {
0179       case hgcee:
0180         rechitMaker_->setADCToGeVConstant(float(hgceeUncalib2GeV_));
0181         cce_correction = hgcEE_cce_[thickness - 1];
0182         sigmaNoiseGeV = 1e-3 * weights_[layer] * rcorr_[thickness] * hgcEE_noise_fC_[thickness - 1] /
0183                         hgcEE_fCPerMIP_[thickness - 1];
0184         break;
0185       case hgcfh:
0186         rechitMaker_->setADCToGeVConstant(float(hgchefUncalib2GeV_));
0187         cce_correction = hgcHEF_cce_[thickness - 1];
0188         sigmaNoiseGeV = 1e-3 * weights_[layer] * rcorr_[thickness + deltasi_index_regemfac_] *
0189                         hgcHEF_noise_fC_[thickness - 1] / hgcHEF_fCPerMIP_[thickness - 1];
0190         break;
0191       case hgcbh:
0192         rechitMaker_->setADCToGeVConstant(float(hgchebUncalib2GeV_));
0193         sigmaNoiseGeV = 1e-3 * hgcHEB_noise_MIP_ * weights_[layer] * rcorrscint_;
0194         break;
0195       case hgchfnose:
0196         rechitMaker_->setADCToGeVConstant(float(hgchfnoseUncalib2GeV_));
0197         cce_correction = hgcHFNose_cce_[thickness - 1];
0198         sigmaNoiseGeV = 1e-3 * weightsNose_[layer] * rcorrNose_[thickness] * hgcHFNose_noise_fC_[thickness - 1] /
0199                         hgcHFNose_fCPerMIP_[thickness - 1];
0200         break;
0201       default:
0202         throw cms::Exception("NonHGCRecHit") << "Rechit with detid = " << detid.rawId() << " is not HGC!";
0203     }
0204
0205     // make the rechit and put in the output collection
0206
0207     HGCRecHit myrechit(rechitMaker_->makeRecHit(uncalibRH, 0));
0208     double new_E = myrechit.energy();
0209     if (detid.det() == DetId::Forward && detid.subdetId() == ForwardSubdetector::HFNose) {
0210       new_E *= (thickness == -1 ? 1.0 : rcorrNose_[thickness]) / cce_correction;
0211     }  //regional factors for silicon in CE_H
0212     else if (idtype == hgcfh) {
0213       new_E *= rcorr_[thickness + deltasi_index_regemfac_] / cce_correction;
0214     }  //regional factors for scintillator and silicon in CE_E
0215     else {
0216       new_E *= (thickness == -1 ? rcorrscint_ : rcorr_[thickness]) / cce_correction;
0217     }
0218
0219     myrechit.setEnergy(new_E);
0220     float SoN = new_E / sigmaNoiseGeV;
0221     myrechit.setSignalOverSigmaNoise(SoN);
0222
0223     if (detid.det() == DetId::HGCalHSc || myrechit.time() == -99.) {
0224       myrechit.setTimeError(-1.);
0225     } else {
0226       float timeError = timeEstimatorSi_.getTimeError("recHit", SoN);
0227       myrechit.setTimeError(timeError);
0228     }
0229
0230     result.push_back(myrechit);
0231   }
0232 }
0233
0234 HGCalRecHitWorkerSimple::~HGCalRecHitWorkerSimple() {}
0235
0236 #include "FWCore/Framework/interface/MakerMacros.h"
0237 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalRecHitWorkerFactory.h"
0238 DEFINE_EDM_PLUGIN(HGCalRecHitWorkerFactory, HGCalRecHitWorkerSimple, "HGCalRecHitWorkerSimple");