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   assert(hgcEE_cce_.size() == hgcEE_noise_fC_.size());
0082   assert(hgcEE_cce_.size() == hgcEE_fCPerMIP_.size());
0083   hgcHEF_noise_fC_ = ps.getParameter<edm::ParameterSet>("HGCHEF_noise_fC").getParameter<std::vector<double> >("values");
0084   hgcHEF_cce_ = ps.getParameter<edm::ParameterSet>("HGCHEF_cce").getParameter<std::vector<double> >("values");
0085   assert(hgcHEF_cce_.size() == hgcHEF_noise_fC_.size());
0086   assert(hgcHEF_cce_.size() == hgcHEF_fCPerMIP_.size());
0087   hgcHEB_noise_MIP_ = ps.getParameter<edm::ParameterSet>("HGCHEB_noise_MIP").getParameter<double>("noise_MIP");
0088   hgcHFNose_noise_fC_ =
0089       ps.getParameter<edm::ParameterSet>("HGCHFNose_noise_fC").getParameter<std::vector<double> >("values");
0090   hgcHFNose_cce_ = ps.getParameter<edm::ParameterSet>("HGCHFNose_cce").getParameter<std::vector<double> >("values");
0091   assert(hgcHFNose_cce_.size() == hgcHFNose_noise_fC_.size());
0092   assert(hgcHFNose_cce_.size() == hgcHFNose_fCPerMIP_.size());
0093
0094   // don't produce rechit if detid is a ghost one
0095   rangeMatch_ = ps.getParameter<uint32_t>("rangeMatch");
0096   rangeMask_ = ps.getParameter<uint32_t>("rangeMask");
0097
0098   // error for recHit time
0099   timeEstimatorSi_ = hgcalsimclustertime::ComputeClusterTime(ps.getParameter<double>("minValSiPar"),
0100                                                              ps.getParameter<double>("maxValSiPar"),
0101                                                              ps.getParameter<double>("constSiPar"),
0102                                                              ps.getParameter<double>("noiseSiPar"));
0103 }
0104
0105 void HGCalRecHitWorkerSimple::set(const edm::EventSetup& es) {
0106   const CaloGeometry& geom = es.getData(caloGeomToken_);
0107   tools_->setGeometry(geom);
0108   rechitMaker_->set(geom);
0109   if (hgcEE_isSiFE_) {
0110     const HGCalGeometry& hgceeGeoHandle = es.getData(ee_geometry_token_);
0111     ddds_[0] = &(hgceeGeoHandle.topology().dddConstants());
0112   } else {
0113     ddds_[0] = nullptr;
0114   }
0115   if (hgcHEF_isSiFE_) {
0116     edm::ESHandle<HGCalGeometry> hgchefGeoHandle = es.getHandle(hef_geometry_token_);
0117     ddds_[1] = &(hgchefGeoHandle->topology().dddConstants());
0118   } else {
0119     ddds_[1] = nullptr;
0120   }
0121   ddds_[2] = nullptr;
0122   if (hgcHFNose_isSiFE_) {
0123     edm::ESHandle<HGCalGeometry> hgchfnoseGeoHandle = es.getHandle(hfnose_geometry_token_);
0124     ddds_[3] = &(hgchfnoseGeoHandle->topology().dddConstants());
0125   } else {
0126     ddds_[3] = nullptr;
0127   }
0128 }
0129
0130 void HGCalRecHitWorkerSimple::run(const edm::Event& evt,
0131                                   const HGCUncalibratedRecHitCollection& uncalibRHColl,
0132                                   HGCRecHitCollection& result) {
0133   for (const auto& uncalibRH : uncalibRHColl) {
0134     DetId detid = uncalibRH.id();
0135     // don't produce rechit if detid is a ghost one
0136
0137     if (detid.det() == DetId::Forward || detid.det() == DetId::Hcal) {
0138       if ((detid & rangeMask_) == rangeMatch_)
0139         continue;
0140     }
0141
0142     int thickness = -1;
0143     float sigmaNoiseGeV = 0.f;
0144     unsigned int layer = tools_->getLayerWithOffset(detid);
0145     float cce_correction = 1.0f;
0146     int idtype(0);
0147
0148     switch (detid.det()) {
0149       case DetId::HGCalEE:
0150         idtype = hgcee;
0151         thickness = 1 + HGCSiliconDetId(detid).type();
0152         break;
0153       case DetId::HGCalHSi:
0154         idtype = hgcfh;
0155         thickness = 1 + HGCSiliconDetId(detid).type();
0156         break;
0157       case DetId::HGCalHSc:
0158         idtype = hgcbh;
0159         break;
0160       default:
0161         switch (detid.subdetId()) {
0162           case HGCEE:
0163             idtype = hgcee;
0164             thickness = ddds_[detid.subdetId() - 3]->waferTypeL(HGCalDetId(detid).wafer());
0165             break;
0166           case HGCHEF:
0167             idtype = hgcfh;
0168             thickness = ddds_[detid.subdetId() - 3]->waferTypeL(HGCalDetId(detid).wafer());
0169             break;
0170           case HcalEndcap:
0171             [[fallthrough]];
0172           case HGCHEB:
0173             idtype = hgcbh;
0174             break;
0175           case HFNose:
0176             idtype = hgchfnose;
0177             thickness = 1 + HFNoseDetId(detid).type();
0178             break;
0179           default:
0180             break;
0181         }
0182     }
0183
0184     switch (idtype) {
0185       case hgcee:
0186         rechitMaker_->setADCToGeVConstant(float(hgceeUncalib2GeV_));
0187         assert(thickness - 1 < static_cast<int>(hgcEE_cce_.size()));
0188         assert(thickness < static_cast<int>(rcorr_.size()));
0189         assert(layer < weights_.size());
0190         cce_correction = hgcEE_cce_[thickness - 1];
0191         sigmaNoiseGeV = 1e-3 * weights_[layer] * rcorr_[thickness] * hgcEE_noise_fC_[thickness - 1] /
0192                         hgcEE_fCPerMIP_[thickness - 1];
0193         break;
0194       case hgcfh:
0195         rechitMaker_->setADCToGeVConstant(float(hgchefUncalib2GeV_));
0196         assert(thickness - 1 < static_cast<int>(hgcHEF_cce_.size()));
0197         assert(thickness + deltasi_index_regemfac_ < static_cast<int>(rcorr_.size()));
0198         assert(layer < weights_.size());
0199         cce_correction = hgcHEF_cce_[thickness - 1];
0200         sigmaNoiseGeV = 1e-3 * weights_[layer] * rcorr_[thickness + deltasi_index_regemfac_] *
0201                         hgcHEF_noise_fC_[thickness - 1] / hgcHEF_fCPerMIP_[thickness - 1];
0202         break;
0203       case hgcbh:
0204         assert(layer < weights_.size());
0205         rechitMaker_->setADCToGeVConstant(float(hgchebUncalib2GeV_));
0206         sigmaNoiseGeV = 1e-3 * hgcHEB_noise_MIP_ * weights_[layer] * rcorrscint_;
0207         break;
0208       case hgchfnose:
0209         rechitMaker_->setADCToGeVConstant(float(hgchfnoseUncalib2GeV_));
0210         assert(thickness - 1 < static_cast<int>(hgcHFNose_cce_.size()));
0211         assert(thickness < static_cast<int>(rcorrNose_.size()));
0212         assert(layer < weightsNose_.size());
0213         cce_correction = hgcHFNose_cce_[thickness - 1];
0214         sigmaNoiseGeV = 1e-3 * weightsNose_[layer] * rcorrNose_[thickness] * hgcHFNose_noise_fC_[thickness - 1] /
0215                         hgcHFNose_fCPerMIP_[thickness - 1];
0216         break;
0217       default:
0218         throw cms::Exception("NonHGCRecHit") << "Rechit with detid = " << detid.rawId() << " is not HGC!";
0219     }
0220
0221     // make the rechit and put in the output collection
0222
0223     HGCRecHit myrechit(rechitMaker_->makeRecHit(uncalibRH, 0));
0224     double new_E = myrechit.energy();
0225     if (detid.det() == DetId::Forward && detid.subdetId() == ForwardSubdetector::HFNose) {
0226       new_E *= (thickness == -1 ? 1.0 : rcorrNose_[thickness]) / cce_correction;
0227     }  //regional factors for silicon in CE_H
0228     else if (idtype == hgcfh) {
0229       new_E *= rcorr_[thickness + deltasi_index_regemfac_] / cce_correction;
0230     }  //regional factors for scintillator and silicon in CE_E
0231     else {
0232       new_E *= (thickness == -1 ? rcorrscint_ : rcorr_[thickness]) / cce_correction;
0233     }
0234
0235     myrechit.setEnergy(new_E);
0236     float SoN = new_E / sigmaNoiseGeV;
0237     myrechit.setSignalOverSigmaNoise(SoN);
0238
0239     if (detid.det() == DetId::HGCalHSc || myrechit.time() == -99.) {
0240       myrechit.setTimeError(-1.);
0241     } else {
0242       float timeError = timeEstimatorSi_.getTimeError("recHit", SoN);
0243       myrechit.setTimeError(timeError);
0244     }
0245
0246     result.push_back(myrechit);
0247   }
0248 }
0249
0250 HGCalRecHitWorkerSimple::~HGCalRecHitWorkerSimple() {}
0251
0252 #include "FWCore/Framework/interface/MakerMacros.h"
0253 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalRecHitWorkerFactory.h"
0254 DEFINE_EDM_PLUGIN(HGCalRecHitWorkerFactory, HGCalRecHitWorkerSimple, "HGCalRecHitWorkerSimple");