Back to home page

Project CMSSW displayed by LXR

 
 

    


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 //#define EDM_ML_DEBUG
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       //sec is wafer and subsec is celltyp
0059       //skip this hit if after ganging it is not valid
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         //assign the RECO DetId
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 }  //namespace
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),  //[CLHEP::c_light]=mm/ns convert to cm/ns
0247       tofDelay_(ps.getParameter<double>("tofDelay")),
0248       averageOccupancies_(occupancyGuesses),
0249       nEvents_(1) {
0250   //configure from cfg
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     //get geometry
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     //valid ID lists
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   // reserve memory for a full detector
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   // release memory for unfilled parts of hash table
0311   if (validIds_.size() * averageOccupancies_[idx] > simHitAccumulator_->size()) {
0312     simHitAccumulator_->reserve(simHitAccumulator_->size());
0313     pusimHitAccumulator_->reserve(simHitAccumulator_->size());
0314   }
0315   //update occupancy guess
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   //get inputs
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   //accumulate in-time the main event
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   //get inputs
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   //accumulate in-time the main event
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   //get inputs
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   //accumulate for the simulated bunch crossing
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     // to be written the verbosity block
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;  // * getCCE(geom, id, cce_);
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         //find position to insert new entry preserving time sorting
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         //cumulate the charge of new entry for all elements that follow in the sorted list
0505         //and resize list accounting for cases when the inserted element itself crosses the threshold
0506 
0507         for (auto step = insertedPos; step != PhitRefs_bx0[id].end(); ++step) {
0508           if (step != insertedPos)
0509             std::get<1>(*(step)) += charge;
0510 
0511           // resize the list stopping with the first timeStamp with cumulative charge above threshold
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   //configuration to apply for the computation of time-of-flight
0557   auto weightToAbyEnergy = theDigitizer_->toaModeByEnergy();
0558   auto tdcForToAOnset = theDigitizer_->tdcForToAOnset();
0559   auto keV2fC = theDigitizer_->keV2fC();
0560 
0561   //create list of tuples (pos in container, RECO DetId, time) to be sorted first
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   //loop over sorted hits
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     //get the data for this cell, if not available then we skip it
0588 
0589     if (!validIds_.count(id))
0590       continue;
0591     HGCSimHitDataAccumulator::iterator simHitIt = simHitAccumulator_->emplace(id, HGCCellInfo()).first;
0592 
0593     if (id == 0)
0594       continue;  // to be ignored at RECO level
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     //distance to the center of the detector
0601     const float dist2center(getPositionDistance(geom, id));
0602 
0603     //hit time: [time()]=ns  [centerDist]=cm [refSpeed_]=cm/ns + delay by 1ns
0604     //accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples)
0605     const float tof = toa - dist2center / refSpeed_ + tofDelay_;
0606     const int itime = std::floor(tof / bxTime_) + 9;
0607 
0608     //no need to add bx crossing - tof comes already corrected from the mixing module
0609     //itime += bxCrossing;
0610     //itime += 9;
0611 
0612     if (itime < 0 || itime > (int)maxBx_)
0613       continue;
0614 
0615     //check if time index is ok and store energy
0616     if (itime >= (int)simHitIt->second.hit_info[0].size())
0617       continue;
0618 
0619     (simHitIt->second).hit_info[0][itime] += charge;
0620 
0621     //for time-of-arrival: save the time-sorted list of timestamps with cumulative charge just above threshold
0622     //working version with pileup only for in-time hits
0623     int waferThickness = getCellThickness(geom, id);
0624     bool orderChanged = false;
0625     if (itime == (int)thisBx_) {
0626       //if start empty => just add charge and time
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         //find position to insert new entry preserving time sorting
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           //just merge timestamps with exact timing
0641           findPos->first += charge;
0642         } else {
0643           //insert new element cumulating the charge
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         //cumulate the charge of new entry for all elements that follow in the sorted list
0651         //and resize list accounting for cases when the inserted element itself crosses the threshold
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           // resize the list stopping with the first timeStamp with cumulative charge above threshold
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         //add new entry at the end of the list
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     //now compute the firing ToA through the interpolation of the consecutive time-stamps at threshold
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   //configuration to apply for the computation of time-of-flight
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 }