|
| 1 | +#include <vector> |
| 2 | +#include <map> |
| 3 | +#include <sstream> |
| 4 | +#include <string> |
| 5 | +#include <cmath> |
| 6 | +#include <utility> |
| 7 | +#include <algorithm> |
| 8 | +#include <fstream> |
| 9 | +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" |
| 10 | +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" |
| 11 | +#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h" |
| 12 | +#include "Geometry/Records/interface/CaloGeometryRecord.h" |
| 13 | +#include "Geometry/CaloTopology/interface/HcalTopology.h" |
| 14 | +#include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h" |
| 15 | +#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h" |
| 16 | +#include "DataFormats/HcalRecHit/interface/HcalRecHitDefs.h" |
| 17 | +#include "EventFilter/HcalRawToDigi/interface/HcalPacker.h" |
| 18 | +#include "DataFormats/FEDRawData/interface/FEDRawDataCollection.h" |
| 19 | +#include "DataFormats/HcalDetId/interface/HcalGenericDetId.h" |
| 20 | +#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h" |
| 21 | +#include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h" |
| 22 | + |
| 23 | +#include "CalibFormats/HcalObjects/interface/HcalDbService.h" |
| 24 | +#include "CalibFormats/HcalObjects/interface/HcalDbRecord.h" |
| 25 | +#include "CalibFormats/HcalObjects/interface/HcalCoderDb.h" |
| 26 | + |
| 27 | + |
| 28 | +#include "CalibFormats/CaloObjects/interface/CaloSamples.h" |
| 29 | +#include "CalibCalorimetry/HcalAlgos/interface/HcalSiPMnonlinearity.h" |
| 30 | +#include "FWCore/Framework/interface/ConsumesCollector.h" |
| 31 | + |
| 32 | +#include "SonicCMS/Core/interface/SonicEDProducer.h" |
| 33 | +#include "SonicCMS/TensorRT/interface/TRTClient.h" |
| 34 | +#include "FWCore/Framework/interface/Event.h" |
| 35 | +#include "FWCore/ParameterSet/interface/ParameterSet.h" |
| 36 | +#include "FWCore/Framework/interface/EventSetup.h" |
| 37 | +#include "FWCore/Framework/interface/ESHandle.h" |
| 38 | +#include "FWCore/Framework/interface/Frameworkfwd.h" |
| 39 | +#include "FWCore/Framework/interface/MakerMacros.h" |
| 40 | +#include "FWCore/MessageLogger/interface/MessageLogger.h" |
| 41 | +#include "FWCore/Framework/interface/stream/EDProducer.h" |
| 42 | +namespace { |
| 43 | + template<class DFrame> |
| 44 | + class RawChargeFromSample |
| 45 | + { |
| 46 | + public: |
| 47 | + inline RawChargeFromSample(const int sipmQTSShift, |
| 48 | + const int sipmQNTStoSum, |
| 49 | + const HcalDbService& cond, |
| 50 | + const HcalDetId id, |
| 51 | + const CaloSamples& cs, |
| 52 | + const int soi, |
| 53 | + const DFrame& frame, |
| 54 | + const int maxTS) {} |
| 55 | + |
| 56 | + inline double getRawCharge(const double decodedCharge, |
| 57 | + const double pedestal) const |
| 58 | + {return decodedCharge;} |
| 59 | + }; |
| 60 | + template<> |
| 61 | + class RawChargeFromSample<QIE11DataFrame> |
| 62 | + { |
| 63 | + public: |
| 64 | + inline RawChargeFromSample(const int sipmQTSShift, |
| 65 | + const int sipmQNTStoSum, |
| 66 | + const HcalDbService& cond, |
| 67 | + const HcalDetId id, |
| 68 | + const CaloSamples& cs, |
| 69 | + const int soi, |
| 70 | + const QIE11DataFrame& frame, |
| 71 | + const int maxTS) |
| 72 | + : siPMParameter_(*cond.getHcalSiPMParameter(id)), |
| 73 | + fcByPE_(siPMParameter_.getFCByPE()), |
| 74 | + corr_(cond.getHcalSiPMCharacteristics()->getNonLinearities(siPMParameter_.getType())) |
| 75 | + { |
| 76 | + if (fcByPE_ <= 0.0) |
| 77 | + throw cms::Exception("HBHEPhase1BadDB") |
| 78 | + << "Invalid fC/PE conversion factor for SiPM " << id |
| 79 | + << std::endl; |
| 80 | + |
| 81 | + const HcalCalibrations& calib = cond.getHcalCalibrations(id); |
| 82 | + const int firstTS = std::max(soi + sipmQTSShift, 0); |
| 83 | + const int lastTS = std::min(firstTS + sipmQNTStoSum, maxTS); |
| 84 | + double sipmQ = 0.0; |
| 85 | + |
| 86 | + for (int ts = firstTS; ts < lastTS; ++ts) |
| 87 | + { |
| 88 | + const double pedestal = calib.pedestal(frame[ts].capid()); |
| 89 | + sipmQ += (cs[ts] - pedestal); |
| 90 | + } |
| 91 | + |
| 92 | + const double effectivePixelsFired = sipmQ/fcByPE_; |
| 93 | + factor_ = corr_.getRecoCorrectionFactor(effectivePixelsFired); |
| 94 | + } |
| 95 | + |
| 96 | + inline double getRawCharge(const double decodedCharge, |
| 97 | + const double pedestal) const |
| 98 | + { |
| 99 | + return (decodedCharge - pedestal)*factor_ + pedestal; |
| 100 | + |
| 101 | + // Old version of TS-by-TS corrections looked as follows: |
| 102 | + // const double sipmQ = decodedCharge - pedestal; |
| 103 | + // const double nPixelsFired = sipmQ/fcByPE_; |
| 104 | + // return sipmQ*corr_.getRecoCorrectionFactor(nPixelsFired) + pedestal; |
| 105 | + } |
| 106 | + |
| 107 | + private: |
| 108 | + const HcalSiPMParameter& siPMParameter_; |
| 109 | + double fcByPE_; |
| 110 | + HcalSiPMnonlinearity corr_; |
| 111 | + double factor_; |
| 112 | + }; |
| 113 | + |
| 114 | +} |
| 115 | + |
| 116 | +template <typename Client> |
| 117 | +class HcalPhase1Reconstructor_FACILE : public SonicEDProducer<Client> |
| 118 | +{ |
| 119 | + public: |
| 120 | + //needed because base class has dependent scope |
| 121 | + using typename SonicEDProducer<Client>::Input; |
| 122 | + using typename SonicEDProducer<Client>::Output; |
| 123 | + explicit HcalPhase1Reconstructor_FACILE(edm::ParameterSet const& cfg) : |
| 124 | + SonicEDProducer<Client>(cfg), |
| 125 | + sipmQTSShift_(cfg.getParameter<unsigned>("sipmQTSShift")), |
| 126 | + sipmQNTStoSum_(cfg.getParameter<unsigned>("sipmQNTStoSum")), |
| 127 | + topN_(cfg.getParameter<unsigned>("topN")), |
| 128 | + fDigiName(cfg.getParameter<edm::InputTag>("digiLabelQIE11")), |
| 129 | + fRHName(cfg.getParameter<edm::InputTag>("edmRecHitName")), |
| 130 | + fChanInfoName(cfg.getParameter<edm::InputTag>("edmChanInfoName")), |
| 131 | + fTokRH(this->template consumes<edm::SortedCollection<HBHERecHit,edm::StrictWeakOrdering<HBHERecHit>> >(fRHName)), |
| 132 | + fTokChanInfo(this->template consumes<edm::SortedCollection<HBHEChannelInfo,edm::StrictWeakOrdering<HBHEChannelInfo>> >(fChanInfoName)), |
| 133 | + fTokDigis(this->template consumes<QIE11DigiCollection>(fDigiName)) |
| 134 | + { |
| 135 | + this->template produces<HBHERecHitCollection>(); |
| 136 | + this->setDebugName("HcalPhase1Reconstructor_FACILE"); |
| 137 | + } |
| 138 | + void acquire(edm::Event const& iEvent, edm::EventSetup const& iSetup, Input& iInput) override { |
| 139 | + |
| 140 | + auto ninput = client_.ninput(); |
| 141 | + auto batchSize = client_.batchSize(); |
| 142 | + iInput = Input(ninput*batchSize, 0.f); |
| 143 | + |
| 144 | + edm::Handle<QIE11DigiCollection> digis; |
| 145 | + iEvent.getByToken(fTokDigis, digis); |
| 146 | + |
| 147 | + edm::ESHandle<HcalDbService> conditions; |
| 148 | + iSetup.get<HcalDbRecord>().get(conditions); |
| 149 | + |
| 150 | + tmp->clear(); |
| 151 | + |
| 152 | + processData<QIE11DataFrame>(*digis, *conditions, iInput, ninput); |
| 153 | + } |
| 154 | + |
| 155 | + template<class DFrame, class Collection> |
| 156 | + void processData(const Collection& coll, |
| 157 | + const HcalDbService& cond, |
| 158 | + Input& iInput, |
| 159 | + auto ninput) |
| 160 | + { |
| 161 | + const bool skipDroppedChannels = false; |
| 162 | + |
| 163 | + unsigned int ib = 0; |
| 164 | + for (typename Collection::const_iterator it = coll.begin(); it != coll.end(); it++){ |
| 165 | + unsigned int ib = std::distance(coll.begin(),it); |
| 166 | + |
| 167 | + const DFrame& frame(*it); |
| 168 | + const HcalDetId cell(frame.id()); |
| 169 | + |
| 170 | + const HcalSubdetector subdet = cell.subdet(); |
| 171 | + if (!(subdet == HcalSubdetector::HcalBarrel || |
| 172 | + subdet == HcalSubdetector::HcalEndcap || |
| 173 | + subdet == HcalSubdetector::HcalOuter)) |
| 174 | + continue; |
| 175 | + |
| 176 | + const HcalCalibrations& calib = cond.getHcalCalibrations(cell); |
| 177 | + const HcalQIECoder* channelCoder = cond.getHcalCoder(cell); |
| 178 | + const HcalQIEShape* shape = cond.getHcalShape(channelCoder); |
| 179 | + const HcalCoderDb coder(*channelCoder, *shape); |
| 180 | + |
| 181 | + CaloSamples cs; |
| 182 | + coder.adc2fC(frame, cs); |
| 183 | + |
| 184 | + const int nRead = cs.size(); |
| 185 | + const int maxTS = std::min(nRead, static_cast<int>(HBHEChannelInfo::MAXSAMPLES)); |
| 186 | + |
| 187 | + const int soi = 3; |
| 188 | + const int nCycles = 8; |
| 189 | + const RawChargeFromSample<DFrame> rcfs(sipmQTSShift_, sipmQNTStoSum_, |
| 190 | + cond, cell, cs, soi, frame, maxTS); |
| 191 | + |
| 192 | + iInput[ib*ninput + 0] = (float)cell.iphi(); |
| 193 | + for (int inputTS = 0; inputTS < nCycles; ++inputTS){ |
| 194 | + auto s(frame[inputTS]); |
| 195 | + const uint8_t adc = s.adc(); |
| 196 | + const int capid = s.capid(); |
| 197 | + const double gain = calib.respcorrgain(capid); |
| 198 | + iInput[ib*ninput + 1] = (float)gain; |
| 199 | + const double rawCharge = rcfs.getRawCharge(cs[inputTS], calib.pedestal(capid)); |
| 200 | + iInput[ib*ninput+inputTS+2] = ((float)rawCharge); |
| 201 | + } |
| 202 | + for(unsigned int d = 1; d < 8; d++){ |
| 203 | + if(cell.depth() == (float)d) { iInput[ib*ninput + d + 9] = 1.; } |
| 204 | + else { iInput[ib*ninput + d + 9] = 0.; } |
| 205 | + } |
| 206 | + for(unsigned int d = 0; d < 30; d++){ |
| 207 | + if(std::abs(cell.ieta()) == (float)d) { iInput[ib*ninput + d + 17] = 1.; } |
| 208 | + else { iInput[ib*ninput + d + 17] = 0.; } |
| 209 | + } |
| 210 | + ib++; |
| 211 | + HBHERecHit rh = HBHERecHit(cell, 0.f,0.f,0.f); |
| 212 | + tmp->push_back(rh); |
| 213 | + } |
| 214 | + } |
| 215 | + void produce(edm::Event& iEvent, edm::EventSetup const& iSetup, Output const& iOutput) override { |
| 216 | + std::unique_ptr<HBHERecHitCollection> out; |
| 217 | + out = std::make_unique<HBHERecHitCollection>(); |
| 218 | + |
| 219 | + unsigned int ib = 0; |
| 220 | + for(HBHERecHitCollection::const_iterator it = tmp->begin(); it != tmp->end(); it++){ |
| 221 | + |
| 222 | + float rh_e = iOutput[ib]; |
| 223 | + if (rh_e < 0.01) rh_e = 0.01; |
| 224 | + else if (rh_e > 1000.) rh_e = 1000.; |
| 225 | + HBHERecHit rhout = HBHERecHit(it->id(),rh_e,0.f,0.f); |
| 226 | + out->push_back(rhout); ib++; |
| 227 | + } |
| 228 | + iEvent.put(std::move(out)); |
| 229 | + } |
| 230 | + ~HcalPhase1Reconstructor_FACILE() override {} |
| 231 | + |
| 232 | + private: |
| 233 | + int sipmQTSShift_; |
| 234 | + int sipmQNTStoSum_; |
| 235 | + unsigned topN_; |
| 236 | + edm::InputTag fDigiName; |
| 237 | + edm::InputTag fRHName; |
| 238 | + edm::InputTag fChanInfoName; |
| 239 | + edm::EDGetTokenT<edm::SortedCollection<HBHERecHit,edm::StrictWeakOrdering<HBHERecHit>>> fTokRH; |
| 240 | + edm::EDGetTokenT<edm::SortedCollection<HBHEChannelInfo,edm::StrictWeakOrdering<HBHEChannelInfo>>> fTokChanInfo; |
| 241 | + edm::EDGetTokenT<QIE11DigiCollection> fTokDigis; |
| 242 | + |
| 243 | + std::vector<HBHERecHit> tmprh; |
| 244 | + std::vector<HBHERecHit> *tmp = &tmprh; |
| 245 | + |
| 246 | + float depth, ieta, iphi; |
| 247 | + |
| 248 | + using SonicEDProducer<Client>::client_; |
| 249 | +}; |
| 250 | + |
| 251 | +typedef HcalPhase1Reconstructor_FACILE<TRTClientSync> HcalPhase1Reconstructor_FACILESync; |
| 252 | +typedef HcalPhase1Reconstructor_FACILE<TRTClientAsync> HcalPhase1Reconstructor_FACILEAsync; |
| 253 | +typedef HcalPhase1Reconstructor_FACILE<TRTClientPseudoAsync> HcalPhase1Reconstructor_FACILEPseudoAsync; |
| 254 | + |
| 255 | +DEFINE_FWK_MODULE(HcalPhase1Reconstructor_FACILESync); |
| 256 | +DEFINE_FWK_MODULE(HcalPhase1Reconstructor_FACILEAsync); |
| 257 | +DEFINE_FWK_MODULE(HcalPhase1Reconstructor_FACILEPseudoAsync); |
0 commit comments