HcalRecProducers/plugins/ZdcHitReconstructor_Run3.cc

0001 #include "ZdcHitReconstructor_Run3.h"
0002 #include "DataFormats/Common/interface/EDCollection.h"
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0007 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
0008 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0009 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0010 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0011 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0012 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0013 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0014
0015 #include <iostream>
0016
0017 #include <Eigen/Dense>
0018
0019 #include <vector>
0020 namespace zdchelper {
0021   void setZDCSaturation(ZDCRecHit rh, QIE10DataFrame& digi, int maxValue) {
0022     for (auto it = digi.begin(); it != digi.end(); it++) {
0023       QIE10DataFrame::Sample sample = QIE10DataFrame::Sample(*it);
0024       if (sample.adc() >= maxValue) {
0025         rh.setFlagField(1, HcalCaloFlagLabels::ADCSaturationBit);
0026         break;
0027       }
0028     }
0029   }
0030
0031 }  // namespace zdchelper
0032
0033 /*  Zdc Hit reconstructor allows for CaloRecHits with status words */
0034 ZdcHitReconstructor_Run3::ZdcHitReconstructor_Run3(edm::ParameterSet const& conf)
0035
0036     : reco_(conf.getParameter<int>("recoMethod")),
0037       saturationFlagSetter_(nullptr),
0038       det_(DetId::Hcal),
0039       correctionMethodEM_(conf.getParameter<int>("correctionMethodEM")),
0040       correctionMethodHAD_(conf.getParameter<int>("correctionMethodHAD")),
0041       correctionMethodRPD_(conf.getParameter<int>("correctionMethodRPD")),
0042       ootpuRatioEM_(conf.getParameter<double>("ootpuRatioEM")),
0043       ootpuRatioHAD_(conf.getParameter<double>("ootpuRatioHAD")),
0044       ootpuRatioRPD_(conf.getParameter<double>("ootpuRatioRPD")),
0045       ootpuFracEM_(conf.getParameter<double>("ootpuFracEM")),
0046       ootpuFracHAD_(conf.getParameter<double>("ootpuFracHAD")),
0047       ootpuFracRPD_(conf.getParameter<double>("ootpuFracRPD")),
0048       chargeRatiosEM_(conf.getParameter<std::vector<double>>("chargeRatiosEM")),
0049       chargeRatiosHAD_(conf.getParameter<std::vector<double>>("chargeRatiosHAD")),
0050       chargeRatiosRPD_(conf.getParameter<std::vector<double>>("chargeRatiosRPD")),
0051       bxTs_(conf.getParameter<std::vector<unsigned int>>("bxTs")),
0052       nTs_(conf.getParameter<int>("nTs")),
0053       forceSOI_(conf.getParameter<bool>("forceSOI")),
0054       signalSOI_(conf.getParameter<std::vector<unsigned int>>("signalSOI")),
0055       noiseSOI_(conf.getParameter<std::vector<unsigned int>>("noiseSOI")),
0056       setSaturationFlags_(conf.getParameter<bool>("setSaturationFlags")),
0057       dropZSmarkedPassed_(conf.getParameter<bool>("dropZSmarkedPassed")),
0058       skipRPD_(conf.getParameter<bool>("skipRPD")) {
0059   tok_input_QIE10 = consumes<QIE10DigiCollection>(conf.getParameter<edm::InputTag>("digiLabelQIE10ZDC"));
0060
0061   std::string subd = conf.getParameter<std::string>("Subdetector");
0062
0063   if (setSaturationFlags_) {
0064     const edm::ParameterSet& pssat = conf.getParameter<edm::ParameterSet>("saturationParameters");
0065     maxADCvalue_ = pssat.getParameter<int>("maxADCvalue");
0066   }
0067   if (!strcasecmp(subd.c_str(), "ZDC")) {
0068     det_ = DetId::Calo;
0069     subdet_ = HcalZDCDetId::SubdetectorId;
0070     produces<ZDCRecHitCollection>();
0071   } else if (!strcasecmp(subd.c_str(), "CALIB")) {
0072     subdet_ = HcalOther;
0073     subdetOther_ = HcalCalibration;
0074     produces<HcalCalibRecHitCollection>();
0075   } else {
0076     std::cout << "ZdcHitReconstructor_Run3 is not associated with a specific subdetector!" << std::endl;
0077   }
0078   reco_.initCorrectionMethod(correctionMethodEM_, 1);
0079   reco_.initCorrectionMethod(correctionMethodHAD_, 2);
0080   reco_.initCorrectionMethod(correctionMethodRPD_, 4);
0081   reco_.initTemplateFit(bxTs_, chargeRatiosEM_, nTs_, 1);
0082   reco_.initTemplateFit(bxTs_, chargeRatiosHAD_, nTs_, 2);
0083   reco_.initTemplateFit(bxTs_, chargeRatiosRPD_, nTs_, 4);
0084   reco_.initRatioSubtraction(ootpuRatioEM_, ootpuFracEM_, 1);
0085   reco_.initRatioSubtraction(ootpuRatioHAD_, ootpuFracHAD_, 2);
0086   reco_.initRatioSubtraction(ootpuRatioRPD_, ootpuFracRPD_, 4);
0087   // ES tokens
0088   htopoToken_ = esConsumes<HcalTopology, HcalRecNumberingRecord, edm::Transition::BeginRun>();
0089   paramsToken_ = esConsumes<HcalLongRecoParams, HcalLongRecoParamsRcd, edm::Transition::BeginRun>();
0090   conditionsToken_ = esConsumes<HcalDbService, HcalDbRecord>();
0091   qualToken_ = esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"));
0092   sevToken_ = esConsumes<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd>();
0093 }
0094
0095 ZdcHitReconstructor_Run3::~ZdcHitReconstructor_Run3() { delete saturationFlagSetter_; }
0096
0097 void ZdcHitReconstructor_Run3::beginRun(edm::Run const& r, edm::EventSetup const& es) {
0098   const HcalTopology& htopo = es.getData(htopoToken_);
0099   const HcalLongRecoParams& p = es.getData(paramsToken_);
0100   longRecoParams_ = std::make_unique<HcalLongRecoParams>(p);
0101   longRecoParams_->setTopo(&htopo);
0102 }
0103
0104 void ZdcHitReconstructor_Run3::endRun(edm::Run const& r, edm::EventSetup const& es) {}
0105
0106 void ZdcHitReconstructor_Run3::produce(edm::Event& e, const edm::EventSetup& eventSetup) {
0107   // get conditions
0108   const HcalDbService* conditions = &eventSetup.getData(conditionsToken_);
0109   const HcalChannelQuality* myqual = &eventSetup.getData(qualToken_);
0110   const HcalSeverityLevelComputer* mySeverity = &eventSetup.getData(sevToken_);
0111
0112   // define vectors to pass noiseTS and signalTS
0113   std::vector<unsigned int> mySignalTS;
0114   std::vector<unsigned int> myNoiseTS;
0115
0116   if (det_ == DetId::Calo && subdet_ == HcalZDCDetId::SubdetectorId) {
0117     edm::Handle<QIE10DigiCollection> digi;
0118     e.getByToken(tok_input_QIE10, digi);
0119
0120     // create empty output
0121     auto rec = std::make_unique<ZDCRecHitCollection>();
0122     rec->reserve(digi->size());
0123
0124     // testing QEI10 conditions
0125     for (auto it = digi->begin(); it != digi->end(); it++) {
0126       QIE10DataFrame QIE10_i = static_cast<QIE10DataFrame>(*it);
0127       HcalZDCDetId cell = QIE10_i.id();
0128       bool isRPD = cell.section() == 4;
0129       if (isRPD && skipRPD_)
0130         continue;
0131       if (cell.section() == 1 && cell.channel() > 5)
0132         continue;  // ignore extra EM channels
0133
0134       DetId detcell = (DetId)cell;
0135
0136       // check on cells to be ignored and dropped: (rof,20.Feb.09)
0137       const HcalChannelStatus* mydigistatus = myqual->getValues(detcell.rawId());
0138       if (mySeverity->dropChannel(mydigistatus->getValue()))
0139         continue;
0140       if (dropZSmarkedPassed_)
0141         if (QIE10_i.zsMarkAndPass())
0142           continue;
0143
0144       const HcalCalibrations& calibrations = conditions->getHcalCalibrations(cell);
0145       const HcalQIECoder* channelCoder = conditions->getHcalCoder(cell);
0146       const HcalQIEShape* shape = conditions->getHcalShape(channelCoder);
0147       HcalCoderDb coder(*channelCoder, *shape);
0148
0149       // pass the effective pedestals to rec hit since both ped value and width used in subtraction of pedestals
0150       const HcalPedestal* effPeds = conditions->getEffectivePedestal(cell);
0151
0152       if (forceSOI_)
0153         rec->push_back(reco_.reconstruct(QIE10_i, noiseSOI_, signalSOI_, coder, calibrations, *effPeds));
0154
0155       else {
0156         const HcalLongRecoParam* myParams = longRecoParams_->getValues(detcell);
0157         mySignalTS.clear();
0158         myNoiseTS.clear();
0159         mySignalTS = myParams->signalTS();
0160         myNoiseTS = myParams->noiseTS();
0161
0162         rec->push_back(reco_.reconstruct(QIE10_i, myNoiseTS, mySignalTS, coder, calibrations, *effPeds));
0163       }
0164       // saturationFlagSetter_ doesn't work with QIE10
0165       // created new function zdchelper::setZDCSaturation to work with QIE10
0166       (rec->back()).setFlags(0);
0167       if (setSaturationFlags_)
0168         zdchelper::setZDCSaturation(rec->back(), QIE10_i, maxADCvalue_);
0169     }
0170     // return result
0171     e.put(std::move(rec));
0172   }  // else if (det_==DetId::Calo...)
0173
0174 }  // void HcalHitReconstructor::produce(...)
0175
0176 void ZdcHitReconstructor_Run3::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0177   // zdcreco
0178   edm::ParameterSetDescription desc;
0179   desc.add<edm::InputTag>("digiLabelQIE10ZDC", edm::InputTag("hcalDigis", "ZDC"));
0180   desc.add<std::string>("Subdetector", "ZDC");
0181   desc.add<bool>("dropZSmarkedPassed", true);
0182   desc.add<bool>("skipRPD", true);
0183   desc.add<int>("recoMethod", 1);
0184   desc.add<int>("correctionMethodEM", 1);
0185   desc.add<int>("correctionMethodHAD", 1);
0186   desc.add<int>("correctionMethodRPD", 0);
0187   desc.add<double>("ootpuRatioEM", 3.0);
0188   desc.add<double>("ootpuRatioHAD", 3.0);
0189   desc.add<double>("ootpuRatioRPD", -1.0);
0190   desc.add<double>("ootpuFracEM", 1.0);
0191   desc.add<double>("ootpuFracHAD", 1.0);
0192   desc.add<double>("ootpuFracRPD", 0.0);
0193   desc.add<std::vector<double>>("chargeRatiosEM",
0194                                 {
0195                                     1.0,
0196                                     0.23157,
0197                                     0.10477,
0198                                     0.06312,
0199                                 });
0200   desc.add<std::vector<double>>("chargeRatiosHAD",
0201                                 {
0202                                     1.0,
0203                                     0.23157,
0204                                     0.10477,
0205                                     0.06312,
0206                                 });
0207   desc.add<std::vector<double>>("chargeRatiosRPD",
0208                                 {
0209                                     1.0,
0210                                     0.23157,
0211                                     0.10477,
0212                                     0.06312,
0213                                 });
0214   desc.add<std::vector<unsigned int>>("bxTs",
0215                                       {
0216                                           0,
0217                                           2,
0218                                           4,
0219                                       });
0220   desc.add<int>("nTs", 6);
0221   desc.add<bool>("forceSOI", false);
0222   desc.add<std::vector<unsigned int>>("signalSOI",
0223                                       {
0224                                           2,
0225                                       });
0226   desc.add<std::vector<unsigned int>>("noiseSOI",
0227                                       {
0228                                           1,
0229                                       });
0230   desc.add<bool>("setSaturationFlags", false);
0231   {
0232     edm::ParameterSetDescription psd0;
0233     psd0.add<int>("maxADCvalue", 255);
0234     desc.add<edm::ParameterSetDescription>("saturationParameters", psd0);
0235   }
0236   descriptions.add("zdcrecoRun3", desc);
0237   // or use the following to generate the label from the module's C++ type
0238   //descriptions.addWithDefaultLabel(desc);
0239 }
0240
0241 //define this as a plug-in
0242 DEFINE_FWK_MODULE(ZdcHitReconstructor_Run3);