HGCalSimProducers/plugins/HGCDigitizer.cc

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