File indexing completed on 2022-05-31 02:17:46
0001 #include "SimCalorimetry/HGCalSimProducers/interface/HGCDigitizer.h"
0002 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0003 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0004 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0005 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0006 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0007 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0008 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0009 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0010 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0011 #include "SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerPluginFactory.h"
0012 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0017 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0018 #include "DataFormats/Math/interface/liblogintpack.h"
0019 #include <algorithm>
0020 #include <boost/foreach.hpp>
0021 #include "FWCore/Utilities/interface/transform.h"
0022
0023
0024 using namespace std;
0025 using namespace hgc_digi;
0026
0027 typedef std::unordered_map<uint32_t, std::vector<std::pair<float, float>>> IdHit_Map;
0028 typedef std::tuple<float, float, float> hit_timeStamp;
0029 typedef std::unordered_map<uint32_t, std::vector<hit_timeStamp>> hitRec_container;
0030
0031 namespace {
0032
0033 constexpr std::array<double, 4> occupancyGuesses = {{0.5, 0.2, 0.2, 0.8}};
0034
0035 float getPositionDistance(const HGCalGeometry* geom, const DetId& id) { return geom->getPosition(id).mag(); }
0036
0037 int getCellThickness(const HGCalGeometry* geom, const DetId& detid) {
0038 const auto& dddConst = geom->topology().dddConstants();
0039 return (1 + dddConst.waferType(detid));
0040 }
0041
0042 void getValidDetIds(const HGCalGeometry* geom, std::unordered_set<DetId>& valid) {
0043 const std::vector<DetId>& ids = geom->getValidDetIds();
0044 valid.reserve(ids.size());
0045 valid.insert(ids.begin(), ids.end());
0046 }
0047
0048 DetId simToReco(const HGCalGeometry* geom, unsigned simId) {
0049 DetId result(0);
0050 const auto& topo = geom->topology();
0051 const auto& dddConst = topo.dddConstants();
0052
0053 if (dddConst.waferHexagon8() || dddConst.tileTrapezoid()) {
0054 result = DetId(simId);
0055 } else {
0056 int subdet(DetId(simId).subdetId()), layer, cell, sec, subsec, zp;
0057 HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell);
0058
0059
0060 auto recoLayerCell = dddConst.simToReco(cell, layer, sec, topo.detectorType());
0061 cell = recoLayerCell.first;
0062 layer = recoLayerCell.second;
0063 if (layer < 0 || cell < 0) {
0064 return result;
0065 } else {
0066
0067 result = HGCalDetId((ForwardSubdetector)subdet, zp, layer, subsec, sec, cell);
0068 }
0069 }
0070 return result;
0071 }
0072
0073 void saveSimHitAccumulator_forPreMix(PHGCSimAccumulator& simResult,
0074 const hgc::HGCPUSimHitDataAccumulator& simData,
0075 const std::unordered_set<DetId>& validIds,
0076 const float minCharge,
0077 const float maxCharge) {
0078 constexpr auto nEnergies = std::tuple_size<decltype(hgc_digi::HGCCellHitInfo().PUhit_info)>::value;
0079 static_assert(nEnergies <= PHGCSimAccumulator::SimHitCollection::energyMask + 1,
0080 "PHGCSimAccumulator bit pattern needs to updated");
0081 static_assert(hgc_digi::nSamples <= PHGCSimAccumulator::SimHitCollection::sampleMask + 1,
0082 "PHGCSimAccumulator bit pattern needs to updated");
0083 const float minPackChargeLog = minCharge > 0.f ? std::log(minCharge) : -2;
0084 const float maxPackChargeLog = std::log(maxCharge);
0085 constexpr uint16_t base = 1 << PHGCSimAccumulator::SimHitCollection::sampleOffset;
0086
0087 simResult.reserve(simData.size());
0088
0089 for (const auto& id : validIds) {
0090 auto found = simData.find(id);
0091 if (found == simData.end())
0092 continue;
0093
0094 const hgc_digi::PUSimHitData& accCharge_across_bxs = found->second.PUhit_info[0];
0095 const hgc_digi::PUSimHitData& timing_across_bxs = found->second.PUhit_info[1];
0096 for (size_t iSample = 0; iSample < hgc_digi::nSamples; ++iSample) {
0097 const std::vector<float>& accCharge_inthis_bx = accCharge_across_bxs[iSample];
0098 const std::vector<float>& timing_inthis_bx = timing_across_bxs[iSample];
0099 std::vector<unsigned short> vc, vt;
0100 size_t nhits = accCharge_inthis_bx.size();
0101
0102 for (size_t ihit = 0; ihit < nhits; ++ihit) {
0103 if (accCharge_inthis_bx[ihit] > minCharge) {
0104 unsigned short c =
0105 logintpack::pack16log(accCharge_inthis_bx[ihit], minPackChargeLog, maxPackChargeLog, base);
0106 unsigned short t = logintpack::pack16log(timing_inthis_bx[ihit], minPackChargeLog, maxPackChargeLog, base);
0107 vc.emplace_back(c);
0108 vt.emplace_back(t);
0109 }
0110 }
0111 simResult.emplace_back(id.rawId(), iSample, vc, vt);
0112 }
0113 }
0114 simResult.shrink_to_fit();
0115 }
0116
0117 void loadSimHitAccumulator_forPreMix(hgc::HGCSimHitDataAccumulator& simData,
0118 hgc::HGCPUSimHitDataAccumulator& PUSimData,
0119 const HGCalGeometry* geom,
0120 IdHit_Map& hitRefs_bx0,
0121 const PHGCSimAccumulator& simAccumulator,
0122 const float minCharge,
0123 const float maxCharge,
0124 bool setIfZero,
0125 const std::array<float, 3>& tdcForToAOnset,
0126 const bool minbiasFlag,
0127 std::unordered_map<uint32_t, bool>& hitOrder_monitor,
0128 const unsigned int thisBx) {
0129 const float minPackChargeLog = minCharge > 0.f ? std::log(minCharge) : -2;
0130 const float maxPackChargeLog = std::log(maxCharge);
0131 constexpr uint16_t base = 1 << PHGCSimAccumulator::SimHitCollection::sampleOffset;
0132 for (const auto& detIdIndexHitInfo : simAccumulator) {
0133 unsigned int detId = detIdIndexHitInfo.detId();
0134
0135 auto simIt = simData.emplace(detId, HGCCellInfo()).first;
0136 size_t nhits = detIdIndexHitInfo.nhits();
0137
0138 hitOrder_monitor[detId] = false;
0139 if (nhits > 0) {
0140 unsigned short iSample = detIdIndexHitInfo.sampleIndex();
0141
0142 const auto& unsigned_charge_array = detIdIndexHitInfo.chargeArray();
0143 const auto& unsigned_time_array = detIdIndexHitInfo.timeArray();
0144
0145 float p_charge, p_time;
0146 unsigned short unsigned_charge, unsigned_time;
0147
0148 for (size_t ihit = 0; ihit < nhits; ++ihit) {
0149 unsigned_charge = (unsigned_charge_array[ihit] & PHGCSimAccumulator::SimHitCollection::dataMask);
0150 unsigned_time = (unsigned_time_array[ihit] & PHGCSimAccumulator::SimHitCollection::dataMask);
0151 p_time = logintpack::unpack16log(unsigned_time, minPackChargeLog, maxPackChargeLog, base);
0152 p_charge = logintpack::unpack16log(unsigned_charge, minPackChargeLog, maxPackChargeLog, base);
0153
0154 (simIt->second).hit_info[0][iSample] += p_charge;
0155 if (iSample == (unsigned short)thisBx) {
0156 if (hitRefs_bx0[detId].empty()) {
0157 hitRefs_bx0[detId].emplace_back(p_charge, p_time);
0158 } else {
0159 if (p_time < hitRefs_bx0[detId].back().second) {
0160 auto findPos = std::upper_bound(hitRefs_bx0[detId].begin(),
0161 hitRefs_bx0[detId].end(),
0162 std::make_pair(0.f, p_time),
0163 [](const auto& i, const auto& j) { return i.second < j.second; });
0164
0165 auto insertedPos = findPos;
0166 if (findPos == hitRefs_bx0[detId].begin()) {
0167 insertedPos = hitRefs_bx0[detId].emplace(insertedPos, p_charge, p_time);
0168 } else {
0169 auto prevPos = findPos - 1;
0170 if (prevPos->second == p_time) {
0171 prevPos->first = prevPos->first + p_charge;
0172 insertedPos = prevPos;
0173 } else if (prevPos->second < p_time) {
0174 insertedPos = hitRefs_bx0[detId].emplace(findPos, (prevPos)->first + p_charge, p_time);
0175 }
0176 }
0177
0178 for (auto step = insertedPos; step != hitRefs_bx0[detId].end(); ++step) {
0179 if (step != insertedPos)
0180 step->first += p_charge;
0181 }
0182
0183 hitOrder_monitor[detId] = true;
0184
0185 } else if (p_time == hitRefs_bx0[detId].back().second) {
0186 hitRefs_bx0[detId].back().first += p_charge;
0187 } else if (p_time > hitRefs_bx0[detId].back().second) {
0188 hitRefs_bx0[detId].emplace_back(hitRefs_bx0[detId].back().first + p_charge, p_time);
0189 }
0190 }
0191 }
0192 }
0193 }
0194 }
0195
0196 if (minbiasFlag) {
0197 for (const auto& hitmapIterator : hitRefs_bx0) {
0198 unsigned int detectorId = hitmapIterator.first;
0199 auto simIt = simData.emplace(detectorId, HGCCellInfo()).first;
0200 const bool orderChanged = hitOrder_monitor[detectorId];
0201 int waferThickness = getCellThickness(geom, detectorId);
0202 float cell_threshold = tdcForToAOnset[waferThickness - 1];
0203 const auto& hitRec = hitmapIterator.second;
0204 float accChargeForToA(0.f), fireTDC(0.f);
0205 const auto aboveThrPos = std::upper_bound(
0206 hitRec.begin(), hitRec.end(), std::make_pair(cell_threshold, 0.f), [](const auto& i, const auto& j) {
0207 return i.first < j.first;
0208 });
0209
0210 if (aboveThrPos == hitRec.end()) {
0211 accChargeForToA = hitRec.back().first;
0212 fireTDC = 0.f;
0213 } else if (hitRec.end() - aboveThrPos > 0 || orderChanged) {
0214 accChargeForToA = aboveThrPos->first;
0215 fireTDC = aboveThrPos->second;
0216 if (aboveThrPos - hitRec.begin() >= 1) {
0217 const auto& belowThrPos = aboveThrPos - 1;
0218 float chargeBeforeThr = belowThrPos->first;
0219 float timeBeforeThr = belowThrPos->second;
0220 float deltaQ = accChargeForToA - chargeBeforeThr;
0221 float deltaTOF = fireTDC - timeBeforeThr;
0222 fireTDC = (cell_threshold - chargeBeforeThr) * deltaTOF / deltaQ + timeBeforeThr;
0223 }
0224 }
0225 (simIt->second).hit_info[1][9] = fireTDC;
0226 }
0227 }
0228 }
0229 }
0230
0231 HGCDigitizer::HGCDigitizer(const edm::ParameterSet& ps, edm::ConsumesCollector& iC)
0232 : simHitAccumulator_(new HGCSimHitDataAccumulator()),
0233 pusimHitAccumulator_(new HGCPUSimHitDataAccumulator()),
0234 digiCollection_(ps.getParameter<std::string>("digiCollection")),
0235 digitizationType_(ps.getParameter<uint32_t>("digitizationType")),
0236 premixStage1_(ps.getParameter<bool>("premixStage1")),
0237 premixStage1MinCharge_(ps.getParameter<double>("premixStage1MinCharge")),
0238 premixStage1MaxCharge_(ps.getParameter<double>("premixStage1MaxCharge")),
0239 maxSimHitsAccTime_(ps.getParameter<uint32_t>("maxSimHitsAccTime")),
0240 bxTime_(ps.getParameter<double>("bxTime")),
0241 hitsProducer_(ps.getParameter<std::string>("hitsProducer")),
0242 hitCollection_(ps.getParameter<std::string>("hitCollection")),
0243 hitToken_(iC.consumes<std::vector<PCaloHit>>(edm::InputTag(hitsProducer_, hitCollection_))),
0244 geomToken_(iC.esConsumes()),
0245 verbosity_(ps.getUntrackedParameter<uint32_t>("verbosity", 0)),
0246 refSpeed_(0.1 * CLHEP::c_light),
0247 tofDelay_(ps.getParameter<double>("tofDelay")),
0248 averageOccupancies_(occupancyGuesses),
0249 nEvents_(1) {
0250
0251
0252 const auto& myCfg_ = ps.getParameter<edm::ParameterSet>("digiCfg");
0253
0254 if (myCfg_.existsAs<edm::ParameterSet>("chargeCollectionEfficiencies")) {
0255 cce_.clear();
0256 const auto& temp = myCfg_.getParameter<edm::ParameterSet>("chargeCollectionEfficiencies")
0257 .getParameter<std::vector<double>>("values");
0258 for (double cce : temp) {
0259 cce_.emplace_back(cce);
0260 }
0261 } else {
0262 std::vector<float>().swap(cce_);
0263 }
0264
0265 auto const& pluginName = ps.getParameter<std::string>("digitizer");
0266 theDigitizer_ = HGCDigitizerPluginFactory::get()->create(pluginName, ps);
0267 }
0268
0269
0270 void HGCDigitizer::initializeEvent(edm::Event const& e, edm::EventSetup const& es) {
0271 if (geomWatcher_.check(es)) {
0272 std::unordered_set<DetId>().swap(validIds_);
0273
0274
0275 CaloGeometry const& geom = es.getData(geomToken_);
0276
0277 gHGCal_ =
0278 dynamic_cast<const HGCalGeometry*>(geom.getSubdetectorGeometry(theDigitizer_->det(), theDigitizer_->subdet()));
0279
0280 int nadded(0);
0281
0282 if (nullptr != gHGCal_) {
0283 getValidDetIds(gHGCal_, validIds_);
0284 } else {
0285 throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0286 }
0287
0288 if (verbosity_ > 0)
0289 edm::LogInfo("HGCDigitizer") << "Added " << nadded << ":" << validIds_.size() << " detIds without "
0290 << hitCollection_ << " in first event processed" << std::endl;
0291 }
0292
0293
0294 unsigned idx = getType();
0295 simHitAccumulator_->reserve(averageOccupancies_[idx] * validIds_.size());
0296 pusimHitAccumulator_->reserve(averageOccupancies_[idx] * validIds_.size());
0297 }
0298
0299
0300 void HGCDigitizer::finalizeEvent(edm::Event& e, edm::EventSetup const& es, CLHEP::HepRandomEngine* hre) {
0301 hitRefs_bx0.clear();
0302 PhitRefs_bx0.clear();
0303 hitOrder_monitor.clear();
0304
0305 const CaloSubdetectorGeometry* theGeom = static_cast<const CaloSubdetectorGeometry*>(gHGCal_);
0306
0307 ++nEvents_;
0308
0309 unsigned idx = getType();
0310
0311 if (validIds_.size() * averageOccupancies_[idx] > simHitAccumulator_->size()) {
0312 simHitAccumulator_->reserve(simHitAccumulator_->size());
0313 pusimHitAccumulator_->reserve(simHitAccumulator_->size());
0314 }
0315
0316 const double thisOcc = simHitAccumulator_->size() / ((double)validIds_.size());
0317 averageOccupancies_[idx] = (averageOccupancies_[idx] * (nEvents_ - 1) + thisOcc) / nEvents_;
0318
0319 if (premixStage1_) {
0320 auto simRecord = std::make_unique<PHGCSimAccumulator>();
0321
0322 if (!pusimHitAccumulator_->empty()) {
0323 saveSimHitAccumulator_forPreMix(
0324 *simRecord, *pusimHitAccumulator_, validIds_, premixStage1MinCharge_, premixStage1MaxCharge_);
0325 }
0326
0327 e.put(std::move(simRecord), digiCollection());
0328
0329 } else {
0330 auto digiResult = std::make_unique<HGCalDigiCollection>();
0331 theDigitizer_->run(digiResult, *simHitAccumulator_, theGeom, validIds_, digitizationType_, hre);
0332 edm::LogVerbatim("HGCDigitizer") << "HGCDigitizer:: finalize event - produced " << digiResult->size()
0333 << " hits in det/subdet " << theDigitizer_->det() << "/"
0334 << theDigitizer_->subdet();
0335 #ifdef EDM_ML_DEBUG
0336 checkPosition(&(*digiResult));
0337 #endif
0338 e.put(std::move(digiResult), digiCollection());
0339 }
0340
0341 hgc::HGCSimHitDataAccumulator().swap(*simHitAccumulator_);
0342 hgc::HGCPUSimHitDataAccumulator().swap(*pusimHitAccumulator_);
0343 }
0344 void HGCDigitizer::accumulate_forPreMix(edm::Event const& e,
0345 edm::EventSetup const& eventSetup,
0346 CLHEP::HepRandomEngine* hre) {
0347
0348
0349 const edm::Handle<edm::PCaloHitContainer>& hits = e.getHandle(hitToken_);
0350 if (!hits.isValid()) {
0351 edm::LogError("HGCDigitizer") << " @ accumulate_minbias : can't find " << hitCollection_ << " collection of "
0352 << hitsProducer_;
0353 return;
0354 }
0355
0356
0357 if (nullptr != gHGCal_) {
0358 accumulate_forPreMix(hits, 0, gHGCal_, hre);
0359 } else {
0360 throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0361 }
0362 }
0363
0364
0365 void HGCDigitizer::accumulate(edm::Event const& e, edm::EventSetup const& eventSetup, CLHEP::HepRandomEngine* hre) {
0366
0367 const edm::Handle<edm::PCaloHitContainer>& hits = e.getHandle(hitToken_);
0368 if (!hits.isValid()) {
0369 edm::LogError("HGCDigitizer") << " @ accumulate : can't find " << hitCollection_ << " collection of "
0370 << hitsProducer_;
0371 return;
0372 }
0373
0374
0375 if (nullptr != gHGCal_) {
0376 accumulate(hits, 0, gHGCal_, hre);
0377 } else {
0378 throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0379 }
0380 }
0381
0382
0383 void HGCDigitizer::accumulate_forPreMix(PileUpEventPrincipal const& e,
0384 edm::EventSetup const& eventSetup,
0385 CLHEP::HepRandomEngine* hre) {
0386 edm::InputTag hitTag(hitsProducer_, hitCollection_);
0387 edm::Handle<edm::PCaloHitContainer> hits;
0388 e.getByLabel(hitTag, hits);
0389
0390 if (!hits.isValid()) {
0391 edm::LogError("HGCDigitizer") << " @ accumulate : can't find " << hitCollection_ << " collection of "
0392 << hitsProducer_;
0393 return;
0394 }
0395
0396 if (nullptr != gHGCal_) {
0397 accumulate_forPreMix(hits, e.bunchCrossing(), gHGCal_, hre);
0398 } else {
0399 throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0400 }
0401 }
0402
0403
0404 void HGCDigitizer::accumulate(PileUpEventPrincipal const& e,
0405 edm::EventSetup const& eventSetup,
0406 CLHEP::HepRandomEngine* hre) {
0407
0408 edm::InputTag hitTag(hitsProducer_, hitCollection_);
0409 edm::Handle<edm::PCaloHitContainer> hits;
0410 e.getByLabel(hitTag, hits);
0411
0412 if (!hits.isValid()) {
0413 edm::LogError("HGCDigitizer") << " @ accumulate : can't find " << hitCollection_ << " collection of "
0414 << hitsProducer_;
0415 return;
0416 }
0417
0418
0419 if (nullptr != gHGCal_) {
0420 accumulate(hits, e.bunchCrossing(), gHGCal_, hre);
0421 } else {
0422 throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0423 }
0424 }
0425
0426
0427 void HGCDigitizer::accumulate_forPreMix(edm::Handle<edm::PCaloHitContainer> const& hits,
0428 int bxCrossing,
0429 const HGCalGeometry* geom,
0430 CLHEP::HepRandomEngine* hre) {
0431 if (nullptr == geom)
0432 return;
0433
0434 auto keV2fC = theDigitizer_->keV2fC();
0435 auto tdcForToAOnset = theDigitizer_->tdcForToAOnset();
0436
0437 int nchits = (int)hits->size();
0438 int count_thisbx = 0;
0439 std::vector<HGCCaloHitTuple_t> hitRefs;
0440 hitRefs.reserve(nchits);
0441 for (int i = 0; i < nchits; ++i) {
0442 const auto& the_hit = hits->at(i);
0443 DetId id = simToReco(geom, the_hit.id());
0444
0445 if (id.rawId() != 0) {
0446 hitRefs.emplace_back(i, id.rawId(), (float)the_hit.time());
0447 }
0448 }
0449 std::sort(hitRefs.begin(), hitRefs.end(), this->orderByDetIdThenTime);
0450
0451 nchits = hitRefs.size();
0452 for (int i = 0; i < nchits; ++i) {
0453 const int hitidx = std::get<0>(hitRefs[i]);
0454 const uint32_t id = std::get<1>(hitRefs[i]);
0455 if (!validIds_.count(id))
0456 continue;
0457
0458 if (id == 0)
0459 continue;
0460
0461 const float toa = std::get<2>(hitRefs[i]);
0462 const PCaloHit& hit = hits->at(hitidx);
0463 const float charge = hit.energy() * 1e6 * keV2fC;
0464
0465 const float dist2center(getPositionDistance(geom, id));
0466 const float tof = toa - dist2center / refSpeed_ + tofDelay_;
0467 const int itime = std::floor(tof / bxTime_) + 9;
0468
0469 if (itime < 0 || itime > (int)maxBx_)
0470 continue;
0471
0472 if (itime >= (int)(maxBx_ + 1))
0473 continue;
0474
0475 int waferThickness = getCellThickness(geom, id);
0476 if (itime == (int)thisBx_) {
0477 ++count_thisbx;
0478 if (PhitRefs_bx0[id].empty()) {
0479 PhitRefs_bx0[id].emplace_back(charge, charge, tof);
0480 } else if (tof > std::get<2>(PhitRefs_bx0[id].back())) {
0481 PhitRefs_bx0[id].emplace_back(charge, charge + std::get<1>(PhitRefs_bx0[id].back()), tof);
0482 } else if (tof == std::get<2>(PhitRefs_bx0[id].back())) {
0483 std::get<0>(PhitRefs_bx0[id].back()) += charge;
0484 std::get<1>(PhitRefs_bx0[id].back()) += charge;
0485 } else {
0486
0487 auto findPos = std::upper_bound(PhitRefs_bx0[id].begin(),
0488 PhitRefs_bx0[id].end(),
0489 hit_timeStamp(charge, 0.f, tof),
0490 [](const auto& i, const auto& j) { return std::get<2>(i) <= std::get<2>(j); });
0491
0492 auto insertedPos = findPos;
0493
0494 if (tof == std::get<2>(*(findPos - 1))) {
0495 std::get<0>(*(findPos - 1)) += charge;
0496 std::get<1>(*(findPos - 1)) += charge;
0497
0498 } else {
0499 insertedPos = PhitRefs_bx0[id].insert(findPos,
0500 (findPos == PhitRefs_bx0[id].begin())
0501 ? hit_timeStamp(charge, charge, tof)
0502 : hit_timeStamp(charge, charge + std::get<1>(*(findPos - 1)), tof));
0503 }
0504
0505
0506
0507 for (auto step = insertedPos; step != PhitRefs_bx0[id].end(); ++step) {
0508 if (step != insertedPos)
0509 std::get<1>(*(step)) += charge;
0510
0511
0512 if (std::get<1>(*step) > tdcForToAOnset[waferThickness - 1] &&
0513 std::get<2>(*step) != std::get<2>(PhitRefs_bx0[id].back())) {
0514 PhitRefs_bx0[id].resize(
0515 std::upper_bound(PhitRefs_bx0[id].begin(),
0516 PhitRefs_bx0[id].end(),
0517 hit_timeStamp(charge, 0.f, std::get<2>(*step)),
0518 [](const auto& i, const auto& j) { return std::get<2>(i) < std::get<2>(j); }) -
0519 PhitRefs_bx0[id].begin());
0520 for (auto stepEnd = step + 1; stepEnd != PhitRefs_bx0[id].end(); ++stepEnd)
0521 std::get<1>(*stepEnd) += charge;
0522 break;
0523 }
0524 }
0525 }
0526 }
0527 }
0528
0529 for (const auto& hitCollection : PhitRefs_bx0) {
0530 const uint32_t detectorId = hitCollection.first;
0531 auto simHitIt = pusimHitAccumulator_->emplace(detectorId, HGCCellHitInfo()).first;
0532
0533 for (const auto& hit_timestamp : PhitRefs_bx0[detectorId]) {
0534 (simHitIt->second).PUhit_info[1][thisBx_].push_back(std::get<2>(hit_timestamp));
0535 (simHitIt->second).PUhit_info[0][thisBx_].push_back(std::get<0>(hit_timestamp));
0536 }
0537 }
0538
0539 if (nchits == 0) {
0540 HGCPUSimHitDataAccumulator::iterator simHitIt = pusimHitAccumulator_->emplace(0, HGCCellHitInfo()).first;
0541 (simHitIt->second).PUhit_info[1][9].push_back(0.0);
0542 (simHitIt->second).PUhit_info[0][9].push_back(0.0);
0543 }
0544 hitRefs.clear();
0545 PhitRefs_bx0.clear();
0546 }
0547
0548
0549 void HGCDigitizer::accumulate(edm::Handle<edm::PCaloHitContainer> const& hits,
0550 int bxCrossing,
0551 const HGCalGeometry* geom,
0552 CLHEP::HepRandomEngine* hre) {
0553 if (nullptr == geom)
0554 return;
0555
0556
0557 auto weightToAbyEnergy = theDigitizer_->toaModeByEnergy();
0558 auto tdcForToAOnset = theDigitizer_->tdcForToAOnset();
0559 auto keV2fC = theDigitizer_->keV2fC();
0560
0561
0562 int nchits = (int)hits->size();
0563
0564 std::vector<HGCCaloHitTuple_t> hitRefs;
0565 hitRefs.reserve(nchits);
0566 for (int i = 0; i < nchits; ++i) {
0567 const auto& the_hit = hits->at(i);
0568 DetId id = simToReco(geom, the_hit.id());
0569
0570 if (verbosity_ > 0) {
0571 edm::LogVerbatim("HGCDigitizer") << "HGCDigitizer::i/p " << std::hex << the_hit.id() << " o/p " << id.rawId()
0572 << std::dec;
0573 }
0574
0575 if (0 != id.rawId()) {
0576 hitRefs.emplace_back(i, id.rawId(), (float)the_hit.time());
0577 }
0578 }
0579
0580 std::sort(hitRefs.begin(), hitRefs.end(), this->orderByDetIdThenTime);
0581
0582 nchits = hitRefs.size();
0583 for (int i = 0; i < nchits; ++i) {
0584 const int hitidx = std::get<0>(hitRefs[i]);
0585 const uint32_t id = std::get<1>(hitRefs[i]);
0586
0587
0588
0589 if (!validIds_.count(id))
0590 continue;
0591 HGCSimHitDataAccumulator::iterator simHitIt = simHitAccumulator_->emplace(id, HGCCellInfo()).first;
0592
0593 if (id == 0)
0594 continue;
0595
0596 const float toa = std::get<2>(hitRefs[i]);
0597 const PCaloHit& hit = hits->at(hitidx);
0598 const float charge = hit.energy() * 1e6 * keV2fC;
0599
0600
0601 const float dist2center(getPositionDistance(geom, id));
0602
0603
0604
0605 const float tof = toa - dist2center / refSpeed_ + tofDelay_;
0606 const int itime = std::floor(tof / bxTime_) + 9;
0607
0608
0609
0610
0611
0612 if (itime < 0 || itime > (int)maxBx_)
0613 continue;
0614
0615
0616 if (itime >= (int)simHitIt->second.hit_info[0].size())
0617 continue;
0618
0619 (simHitIt->second).hit_info[0][itime] += charge;
0620
0621
0622
0623 int waferThickness = getCellThickness(geom, id);
0624 bool orderChanged = false;
0625 if (itime == (int)thisBx_) {
0626
0627 if (hitRefs_bx0[id].empty()) {
0628 hitRefs_bx0[id].emplace_back(charge, tof);
0629
0630 } else if (tof <= hitRefs_bx0[id].back().second) {
0631
0632 std::vector<std::pair<float, float>>::iterator findPos =
0633 std::upper_bound(hitRefs_bx0[id].begin(),
0634 hitRefs_bx0[id].end(),
0635 std::pair<float, float>(0.f, tof),
0636 [](const auto& i, const auto& j) { return i.second <= j.second; });
0637
0638 std::vector<std::pair<float, float>>::iterator insertedPos = findPos;
0639 if (findPos->second == tof) {
0640
0641 findPos->first += charge;
0642 } else {
0643
0644 insertedPos = hitRefs_bx0[id].insert(findPos,
0645 (findPos == hitRefs_bx0[id].begin())
0646 ? std::pair<float, float>(charge, tof)
0647 : std::pair<float, float>((findPos - 1)->first + charge, tof));
0648 }
0649
0650
0651
0652 for (std::vector<std::pair<float, float>>::iterator step = insertedPos; step != hitRefs_bx0[id].end(); ++step) {
0653 if (step != insertedPos)
0654 step->first += charge;
0655
0656 if (step->first > tdcForToAOnset[waferThickness - 1] && step->second != hitRefs_bx0[id].back().second) {
0657 hitRefs_bx0[id].resize(std::upper_bound(hitRefs_bx0[id].begin(),
0658 hitRefs_bx0[id].end(),
0659 std::pair<float, float>(0.f, step->second),
0660 [](const auto& i, const auto& j) { return i.second < j.second; }) -
0661 hitRefs_bx0[id].begin());
0662 for (auto stepEnd = step + 1; stepEnd != hitRefs_bx0[id].end(); ++stepEnd)
0663 stepEnd->first += charge;
0664 break;
0665 }
0666 }
0667
0668 orderChanged = true;
0669 } else {
0670
0671 if (hitRefs_bx0[id].back().first <= tdcForToAOnset[waferThickness - 1]) {
0672 hitRefs_bx0[id].emplace_back(hitRefs_bx0[id].back().first + charge, tof);
0673 }
0674 }
0675 }
0676 float accChargeForToA = hitRefs_bx0[id].empty() ? 0.f : hitRefs_bx0[id].back().first;
0677
0678 if (weightToAbyEnergy)
0679 (simHitIt->second).hit_info[1][itime] += charge * tof;
0680 else if (accChargeForToA > tdcForToAOnset[waferThickness - 1] &&
0681 ((simHitIt->second).hit_info[1][itime] == 0 || orderChanged == true)) {
0682 float fireTDC = hitRefs_bx0[id].back().second;
0683 if (hitRefs_bx0[id].size() > 1) {
0684 float chargeBeforeThr = (hitRefs_bx0[id].end() - 2)->first;
0685 float tofchargeBeforeThr = (hitRefs_bx0[id].end() - 2)->second;
0686
0687 float deltaQ = accChargeForToA - chargeBeforeThr;
0688 float deltaTOF = fireTDC - tofchargeBeforeThr;
0689 fireTDC = (tdcForToAOnset[waferThickness - 1] - chargeBeforeThr) * deltaTOF / deltaQ + tofchargeBeforeThr;
0690 }
0691 (simHitIt->second).hit_info[1][itime] = fireTDC;
0692 }
0693 }
0694 hitRefs.clear();
0695 }
0696 void HGCDigitizer::accumulate_forPreMix(const PHGCSimAccumulator& simAccumulator, const bool minbiasFlag) {
0697
0698 auto weightToAbyEnergy = theDigitizer_->toaModeByEnergy();
0699 auto tdcForToAOnset = theDigitizer_->tdcForToAOnset();
0700
0701 if (nullptr != gHGCal_) {
0702 loadSimHitAccumulator_forPreMix(*simHitAccumulator_,
0703 *pusimHitAccumulator_,
0704 gHGCal_,
0705 hitRefs_bx0,
0706 simAccumulator,
0707 premixStage1MinCharge_,
0708 premixStage1MaxCharge_,
0709 !weightToAbyEnergy,
0710 tdcForToAOnset,
0711 minbiasFlag,
0712 hitOrder_monitor,
0713 thisBx_);
0714 }
0715 }
0716
0717
0718 void HGCDigitizer::resetSimHitDataAccumulator() {
0719 for (HGCSimHitDataAccumulator::iterator it = simHitAccumulator_->begin(); it != simHitAccumulator_->end(); it++) {
0720 it->second.hit_info[0].fill(0.);
0721 it->second.hit_info[1].fill(0.);
0722 }
0723 }
0724
0725 uint32_t HGCDigitizer::getType() const {
0726 uint32_t idx = std::numeric_limits<unsigned>::max();
0727 switch (theDigitizer_->det()) {
0728 case DetId::HGCalEE:
0729 idx = 0;
0730 break;
0731 case DetId::HGCalHSi:
0732 idx = 1;
0733 break;
0734 case DetId::HGCalHSc:
0735 idx = 2;
0736 break;
0737 case DetId::Forward:
0738 idx = 3;
0739 break;
0740 default:
0741 break;
0742 }
0743 return idx;
0744 }
0745
0746 void HGCDigitizer::checkPosition(const HGCalDigiCollection* digis) const {
0747 const double tol(0.5);
0748 if (nullptr != gHGCal_) {
0749 for (const auto& digi : *(digis)) {
0750 const DetId& id = digi.id();
0751 const GlobalPoint& global = gHGCal_->getPosition(id);
0752 double r = global.perp();
0753 double z = std::abs(global.z());
0754 std::pair<double, double> zrange = gHGCal_->topology().dddConstants().rangeZ(true);
0755 std::pair<double, double> rrange = gHGCal_->topology().dddConstants().rangeR(z, true);
0756 bool ok = ((r >= rrange.first) && (r <= rrange.second) && (z >= zrange.first) && (z <= zrange.second));
0757 std::string ck = (((r < rrange.first - tol) || (r > rrange.second + tol) || (z < zrange.first - tol) ||
0758 (z > zrange.second + tol))
0759 ? "***** ERROR *****"
0760 : "");
0761 bool val = gHGCal_->topology().valid(id);
0762 if ((!ok) || (!val)) {
0763 if (id.det() == DetId::HGCalEE || id.det() == DetId::HGCalHSi) {
0764 edm::LogVerbatim("HGCDigitizer") << "Check " << HGCSiliconDetId(id) << " " << global << " R " << r << ":"
0765 << rrange.first << ":" << rrange.second << " Z " << z << ":" << zrange.first
0766 << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0767 } else if (id.det() == DetId::HGCalHSc) {
0768 edm::LogVerbatim("HGCDigitizer") << "Check " << HGCScintillatorDetId(id) << " " << global << " R " << r << ":"
0769 << rrange.first << ":" << rrange.second << " Z " << z << ":" << zrange.first
0770 << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0771 } else if ((id.det() == DetId::Forward) && (id.subdetId() == static_cast<int>(HFNose))) {
0772 edm::LogVerbatim("HGCDigitizer") << "Check " << HFNoseDetId(id) << " " << global << " R " << r << ":"
0773 << rrange.first << ":" << rrange.second << " Z " << z << ":" << zrange.first
0774 << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0775 } else {
0776 edm::LogVerbatim("HGCDigitizer")
0777 << "Check " << std::hex << id.rawId() << std::dec << " " << id.det() << ":" << id.subdetId() << " "
0778 << global << " R " << r << ":" << rrange.first << ":" << rrange.second << " Z " << z << ":"
0779 << zrange.first << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0780 }
0781 }
0782 }
0783 }
0784 }