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 bool HGCalRecHitWorkerSimple::run(const edm::Event& evt,
0125                                   const HGCUncalibratedRecHit& uncalibRH,
0126                                   HGCRecHitCollection& result) {
0127   DetId detid = uncalibRH.id();
0128   // don't produce rechit if detid is a ghost one
0129
0130   if (detid.det() == DetId::Forward || detid.det() == DetId::Hcal) {
0131     if ((detid & rangeMask_) == rangeMatch_)
0132       return false;
0133   }
0134
0135   int thickness = -1;
0136   float sigmaNoiseGeV = 0.f;
0137   unsigned int layer = tools_->getLayerWithOffset(detid);
0138   float cce_correction = 1.0;
0139   int idtype(0);
0140
0141   switch (detid.det()) {
0142     case DetId::HGCalEE:
0143       idtype = hgcee;
0144       thickness = 1 + HGCSiliconDetId(detid).type();
0145       break;
0146     case DetId::HGCalHSi:
0147       idtype = hgcfh;
0148       thickness = 1 + HGCSiliconDetId(detid).type();
0149       break;
0150     case DetId::HGCalHSc:
0151       idtype = hgcbh;
0152       break;
0153     default:
0154       switch (detid.subdetId()) {
0155         case HGCEE:
0156           idtype = hgcee;
0157           thickness = ddds_[detid.subdetId() - 3]->waferTypeL(HGCalDetId(detid).wafer());
0158           break;
0159         case HGCHEF:
0160           idtype = hgcfh;
0161           thickness = ddds_[detid.subdetId() - 3]->waferTypeL(HGCalDetId(detid).wafer());
0162           break;
0163         case HcalEndcap:
0164           [[fallthrough]];
0165         case HGCHEB:
0166           idtype = hgcbh;
0167           break;
0168         case HFNose:
0169           idtype = hgchfnose;
0170           thickness = 1 + HFNoseDetId(detid).type();
0171           break;
0172         default:
0173           break;
0174       }
0175   }
0176
0177   switch (idtype) {
0178     case hgcee:
0179       rechitMaker_->setADCToGeVConstant(float(hgceeUncalib2GeV_));
0180       cce_correction = hgcEE_cce_[thickness - 1];
0181       sigmaNoiseGeV =
0182           1e-3 * weights_[layer] * rcorr_[thickness] * hgcEE_noise_fC_[thickness - 1] / hgcEE_fCPerMIP_[thickness - 1];
0183       break;
0184     case hgcfh:
0185       rechitMaker_->setADCToGeVConstant(float(hgchefUncalib2GeV_));
0186       cce_correction = hgcHEF_cce_[thickness - 1];
0187       sigmaNoiseGeV = 1e-3 * weights_[layer] * rcorr_[thickness + deltasi_index_regemfac_] *
0188                       hgcHEF_noise_fC_[thickness - 1] / hgcHEF_fCPerMIP_[thickness - 1];
0189       break;
0190     case hgcbh:
0191       rechitMaker_->setADCToGeVConstant(float(hgchebUncalib2GeV_));
0192       sigmaNoiseGeV = 1e-3 * hgcHEB_noise_MIP_ * weights_[layer] * rcorrscint_;
0193       break;
0194     case hgchfnose:
0195       rechitMaker_->setADCToGeVConstant(float(hgchfnoseUncalib2GeV_));
0196       cce_correction = hgcHFNose_cce_[thickness - 1];
0197       sigmaNoiseGeV = 1e-3 * weightsNose_[layer] * rcorrNose_[thickness] * hgcHFNose_noise_fC_[thickness - 1] /
0198                       hgcHFNose_fCPerMIP_[thickness - 1];
0199       break;
0200     default:
0201       throw cms::Exception("NonHGCRecHit") << "Rechit with detid = " << detid.rawId() << " is not HGC!";
0202   }
0203
0204   // make the rechit and put in the output collection
0205
0206   HGCRecHit myrechit(rechitMaker_->makeRecHit(uncalibRH, 0));
0207   double new_E = myrechit.energy();
0208   if (detid.det() == DetId::Forward && detid.subdetId() == ForwardSubdetector::HFNose) {
0209     new_E *= (thickness == -1 ? 1.0 : rcorrNose_[thickness]) / cce_correction;
0210   }  //regional factors for silicon in CE_H
0211   else if (idtype == hgcfh) {
0212     new_E *= rcorr_[thickness + deltasi_index_regemfac_] / cce_correction;
0213   }  //regional factors for scintillator and silicon in CE_E
0214   else {
0215     new_E *= (thickness == -1 ? rcorrscint_ : rcorr_[thickness]) / cce_correction;
0216   }
0217
0218   myrechit.setEnergy(new_E);
0219   float SoN = new_E / sigmaNoiseGeV;
0220   myrechit.setSignalOverSigmaNoise(SoN);
0221
0222   if (detid.det() == DetId::HGCalHSc || myrechit.time() < 0.) {
0223     myrechit.setTimeError(-1.);
0224   } else {
0225     float timeError = timeEstimatorSi_.getTimeError("recHit", SoN);
0226     myrechit.setTimeError(timeError);
0227   }
0228
0229   result.push_back(myrechit);
0230
0231   return true;
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");