Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-11 03:00:22 

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/AbstractServices/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         assert(waferThickness - 1 < static_cast<int>(tdcForToAOnset.size()));
0200         float cell_threshold = tdcForToAOnset[waferThickness - 1];
0201         const auto& hitRec = hitmapIterator.second;
0202         float fireTDC(0.f);
0203         const auto aboveThrPos = std::upper_bound(
0204             hitRec.begin(), hitRec.end(), std::make_pair(cell_threshold, 0.f), [](const auto& i, const auto& j) {
0205               return i.first < j.first;
0206             });
0207 
0208         if (aboveThrPos != hitRec.end()) {
0209           if (hitRec.end() - aboveThrPos > 0 || orderChanged) {
0210             fireTDC = aboveThrPos->second;
0211             if (aboveThrPos - hitRec.begin() >= 1) {
0212               float accChargeForToA = aboveThrPos->first;
0213               const auto& belowThrPos = aboveThrPos - 1;
0214               float chargeBeforeThr = belowThrPos->first;
0215               float timeBeforeThr = belowThrPos->second;
0216               float deltaQ = accChargeForToA - chargeBeforeThr;
0217               float deltaTOF = fireTDC - timeBeforeThr;
0218               fireTDC = (cell_threshold - chargeBeforeThr) * deltaTOF / deltaQ + timeBeforeThr;
0219             }
0220           }
0221         }
0222         (simIt->second).hit_info[1][9] = fireTDC;
0223       }
0224     }
0225   }
0226 }  //namespace
0227 
0228 HGCDigitizer::HGCDigitizer(const edm::ParameterSet& ps, edm::ConsumesCollector& iC)
0229     : simHitAccumulator_(new HGCSimHitDataAccumulator()),
0230       pusimHitAccumulator_(new HGCPUSimHitDataAccumulator()),
0231       digiCollection_(ps.getParameter<std::string>("digiCollection")),
0232       digitizationType_(ps.getParameter<uint32_t>("digitizationType")),
0233       premixStage1_(ps.getParameter<bool>("premixStage1")),
0234       premixStage1MinCharge_(ps.getParameter<double>("premixStage1MinCharge")),
0235       premixStage1MaxCharge_(ps.getParameter<double>("premixStage1MaxCharge")),
0236       maxSimHitsAccTime_(ps.getParameter<uint32_t>("maxSimHitsAccTime")),
0237       bxTime_(ps.getParameter<double>("bxTime")),
0238       hitsProducer_(ps.getParameter<std::string>("hitsProducer")),
0239       hitCollection_(ps.getParameter<std::string>("hitCollection")),
0240       hitToken_(iC.consumes<std::vector<PCaloHit>>(edm::InputTag(hitsProducer_, hitCollection_))),
0241       geomToken_(iC.esConsumes()),
0242       verbosity_(ps.getUntrackedParameter<uint32_t>("verbosity", 0)),
0243       tofDelay_(ps.getParameter<double>("tofDelay")),
0244       averageOccupancies_(occupancyGuesses),
0245       nEvents_(1) {
0246   //configure from cfg
0247 
0248   const auto& myCfg_ = ps.getParameter<edm::ParameterSet>("digiCfg");
0249 
0250   if (myCfg_.existsAs<edm::ParameterSet>("chargeCollectionEfficiencies")) {
0251     cce_.clear();
0252     const auto& temp = myCfg_.getParameter<edm::ParameterSet>("chargeCollectionEfficiencies")
0253                            .getParameter<std::vector<double>>("values");
0254     for (double cce : temp) {
0255       cce_.emplace_back(cce);
0256     }
0257   } else {
0258     std::vector<float>().swap(cce_);
0259   }
0260 
0261   auto const& pluginName = ps.getParameter<std::string>("digitizer");
0262   theDigitizer_ = HGCDigitizerPluginFactory::get()->create(pluginName, ps);
0263 }
0264 
0265 //
0266 void HGCDigitizer::initializeEvent(edm::Event const& e, edm::EventSetup const& es) {
0267   if (geomWatcher_.check(es)) {
0268     std::unordered_set<DetId>().swap(validIds_);
0269 
0270     //get geometry
0271     CaloGeometry const& geom = es.getData(geomToken_);
0272 
0273     gHGCal_ =
0274         dynamic_cast<const HGCalGeometry*>(geom.getSubdetectorGeometry(theDigitizer_->det(), theDigitizer_->subdet()));
0275 
0276     int nadded(0);
0277     //valid ID lists
0278     if (nullptr != gHGCal_) {
0279       getValidDetIds(gHGCal_, validIds_);
0280     } else {
0281       throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0282     }
0283 
0284     if (verbosity_ > 0)
0285       edm::LogInfo("HGCDigitizer") << "Added " << nadded << ":" << validIds_.size() << " detIds without "
0286                                    << hitCollection_ << " in first event processed" << std::endl;
0287   }
0288 
0289   // reserve memory for a full detector
0290   unsigned idx = getType();
0291   simHitAccumulator_->reserve(averageOccupancies_[idx] * validIds_.size());
0292   pusimHitAccumulator_->reserve(averageOccupancies_[idx] * validIds_.size());
0293 }
0294 
0295 //
0296 void HGCDigitizer::finalizeEvent(edm::Event& e, edm::EventSetup const& es, CLHEP::HepRandomEngine* hre) {
0297   hitRefs_bx0.clear();
0298   PhitRefs_bx0.clear();
0299   hitOrder_monitor.clear();
0300 
0301   const CaloSubdetectorGeometry* theGeom = static_cast<const CaloSubdetectorGeometry*>(gHGCal_);
0302 
0303   ++nEvents_;
0304 
0305   unsigned idx = getType();
0306   // release memory for unfilled parts of hash table
0307   if (validIds_.size() * averageOccupancies_[idx] > simHitAccumulator_->size()) {
0308     simHitAccumulator_->reserve(simHitAccumulator_->size());
0309     pusimHitAccumulator_->reserve(simHitAccumulator_->size());
0310   }
0311   //update occupancy guess
0312   const double thisOcc = simHitAccumulator_->size() / ((double)validIds_.size());
0313   averageOccupancies_[idx] = (averageOccupancies_[idx] * (nEvents_ - 1) + thisOcc) / nEvents_;
0314 
0315   if (premixStage1_) {
0316     auto simRecord = std::make_unique<PHGCSimAccumulator>();
0317 
0318     if (!pusimHitAccumulator_->empty()) {
0319       saveSimHitAccumulator_forPreMix(
0320           *simRecord, *pusimHitAccumulator_, validIds_, premixStage1MinCharge_, premixStage1MaxCharge_);
0321     }
0322 
0323     e.put(std::move(simRecord), digiCollection());
0324 
0325   } else {
0326     auto digiResult = std::make_unique<HGCalDigiCollection>();
0327     theDigitizer_->run(digiResult, *simHitAccumulator_, theGeom, validIds_, digitizationType_, hre);
0328     edm::LogVerbatim("HGCDigitizer") << "HGCDigitizer:: finalize event - produced " << digiResult->size()
0329                                      << " hits in det/subdet " << theDigitizer_->det() << "/"
0330                                      << theDigitizer_->subdet();
0331 #ifdef EDM_ML_DEBUG
0332     checkPosition(&(*digiResult));
0333 #endif
0334     e.put(std::move(digiResult), digiCollection());
0335   }
0336 
0337   hgc::HGCSimHitDataAccumulator().swap(*simHitAccumulator_);
0338   hgc::HGCPUSimHitDataAccumulator().swap(*pusimHitAccumulator_);
0339 }
0340 void HGCDigitizer::accumulate_forPreMix(edm::Event const& e,
0341                                         edm::EventSetup const& eventSetup,
0342                                         CLHEP::HepRandomEngine* hre) {
0343   //get inputs
0344 
0345   const edm::Handle<edm::PCaloHitContainer>& hits = e.getHandle(hitToken_);
0346   if (!hits.isValid()) {
0347     edm::LogError("HGCDigitizer") << " @ accumulate_minbias : can't find " << hitCollection_ << " collection of "
0348                                   << hitsProducer_;
0349     return;
0350   }
0351 
0352   //accumulate in-time the main event
0353   if (nullptr != gHGCal_) {
0354     accumulate_forPreMix(hits, 0, gHGCal_, hre);
0355   } else {
0356     throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0357   }
0358 }
0359 
0360 //
0361 void HGCDigitizer::accumulate(edm::Event const& e, edm::EventSetup const& eventSetup, CLHEP::HepRandomEngine* hre) {
0362   //get inputs
0363   const edm::Handle<edm::PCaloHitContainer>& hits = e.getHandle(hitToken_);
0364   if (!hits.isValid()) {
0365     edm::LogError("HGCDigitizer") << " @ accumulate : can't find " << hitCollection_ << " collection of "
0366                                   << hitsProducer_;
0367     return;
0368   }
0369 
0370   //accumulate in-time the main event
0371   if (nullptr != gHGCal_) {
0372     accumulate(hits, 0, gHGCal_, hre);
0373   } else {
0374     throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0375   }
0376 }
0377 
0378 //
0379 void HGCDigitizer::accumulate_forPreMix(PileUpEventPrincipal const& e,
0380                                         edm::EventSetup const& eventSetup,
0381                                         CLHEP::HepRandomEngine* hre) {
0382   edm::InputTag hitTag(hitsProducer_, hitCollection_);
0383   edm::Handle<edm::PCaloHitContainer> hits;
0384   e.getByLabel(hitTag, hits);
0385 
0386   if (!hits.isValid()) {
0387     edm::LogError("HGCDigitizer") << " @ accumulate : can't find " << hitCollection_ << " collection of "
0388                                   << hitsProducer_;
0389     return;
0390   }
0391 
0392   if (nullptr != gHGCal_) {
0393     accumulate_forPreMix(hits, e.bunchCrossing(), gHGCal_, hre);
0394   } else {
0395     throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0396   }
0397 }
0398 
0399 //
0400 void HGCDigitizer::accumulate(PileUpEventPrincipal const& e,
0401                               edm::EventSetup const& eventSetup,
0402                               CLHEP::HepRandomEngine* hre) {
0403   //get inputs
0404   edm::InputTag hitTag(hitsProducer_, hitCollection_);
0405   edm::Handle<edm::PCaloHitContainer> hits;
0406   e.getByLabel(hitTag, hits);
0407 
0408   if (!hits.isValid()) {
0409     edm::LogError("HGCDigitizer") << " @ accumulate : can't find " << hitCollection_ << " collection of "
0410                                   << hitsProducer_;
0411     return;
0412   }
0413 
0414   //accumulate for the simulated bunch crossing
0415   if (nullptr != gHGCal_) {
0416     accumulate(hits, e.bunchCrossing(), gHGCal_, hre);
0417   } else {
0418     throw cms::Exception("BadConfiguration") << "HGCDigitizer is not producing EE, FH, or BH digis!";
0419   }
0420 }
0421 
0422 //
0423 void HGCDigitizer::accumulate_forPreMix(edm::Handle<edm::PCaloHitContainer> const& hits,
0424                                         int bxCrossing,
0425                                         const HGCalGeometry* geom,
0426                                         CLHEP::HepRandomEngine* hre) {
0427   if (nullptr == geom)
0428     return;
0429 
0430   auto keV2fC = theDigitizer_->keV2fC();
0431   auto tdcForToAOnset = theDigitizer_->tdcForToAOnset();
0432 
0433   int nchits = (int)hits->size();
0434   std::vector<HGCCaloHitTuple_t> hitRefs;
0435   hitRefs.reserve(nchits);
0436   for (int i = 0; i < nchits; ++i) {
0437     const auto& the_hit = hits->at(i);
0438     DetId id = simToReco(geom, the_hit.id());
0439     // to be written the verbosity block
0440     if (id.rawId() != 0) {
0441       hitRefs.emplace_back(i, id.rawId(), (float)the_hit.time());
0442     }
0443   }
0444   std::sort(hitRefs.begin(), hitRefs.end(), this->orderByDetIdThenTime);
0445 
0446   nchits = hitRefs.size();
0447   for (int i = 0; i < nchits; ++i) {
0448     const int hitidx = std::get<0>(hitRefs[i]);
0449     const uint32_t id = std::get<1>(hitRefs[i]);
0450     if (!validIds_.count(id))
0451       continue;
0452 
0453     if (id == 0)
0454       continue;
0455 
0456     const float toa = std::get<2>(hitRefs[i]);
0457     const PCaloHit& hit = hits->at(hitidx);
0458     const float charge = hit.energy() * 1e6 * keV2fC;  // * getCCE(geom, id, cce_);
0459 
0460     const float tof = toa + tofDelay_;
0461     const int itime = std::floor(tof / bxTime_) + 9;
0462 
0463     if (itime < 0 || itime > (int)maxBx_)
0464       continue;
0465 
0466     if (itime >= (int)(maxBx_ + 1))
0467       continue;
0468 
0469     int waferThickness = getCellThickness(geom, id);
0470     if (itime == (int)thisBx_) {
0471       if (PhitRefs_bx0[id].empty()) {
0472         PhitRefs_bx0[id].emplace_back(charge, charge, tof);
0473       } else if (tof > std::get<2>(PhitRefs_bx0[id].back())) {
0474         PhitRefs_bx0[id].emplace_back(charge, charge + std::get<1>(PhitRefs_bx0[id].back()), tof);
0475       } else if (tof == std::get<2>(PhitRefs_bx0[id].back())) {
0476         std::get<0>(PhitRefs_bx0[id].back()) += charge;
0477         std::get<1>(PhitRefs_bx0[id].back()) += charge;
0478       } else {
0479         //find position to insert new entry preserving time sorting
0480         auto findPos = std::upper_bound(PhitRefs_bx0[id].begin(),
0481                                         PhitRefs_bx0[id].end(),
0482                                         hit_timeStamp(charge, 0.f, tof),
0483                                         [](const auto& i, const auto& j) { return std::get<2>(i) <= std::get<2>(j); });
0484 
0485         auto insertedPos = findPos;
0486 
0487         if (tof == std::get<2>(*(findPos - 1))) {
0488           std::get<0>(*(findPos - 1)) += charge;
0489           std::get<1>(*(findPos - 1)) += charge;
0490 
0491         } else {
0492           insertedPos = PhitRefs_bx0[id].insert(findPos,
0493                                                 (findPos == PhitRefs_bx0[id].begin())
0494                                                     ? hit_timeStamp(charge, charge, tof)
0495                                                     : hit_timeStamp(charge, charge + std::get<1>(*(findPos - 1)), tof));
0496         }
0497         //cumulate the charge of new entry for all elements that follow in the sorted list
0498         //and resize list accounting for cases when the inserted element itself crosses the threshold
0499 
0500         assert(waferThickness - 1 < static_cast<int>(tdcForToAOnset.size()));
0501         for (auto step = insertedPos; step != PhitRefs_bx0[id].end(); ++step) {
0502           if (step != insertedPos)
0503             std::get<1>(*(step)) += charge;
0504 
0505           // resize the list stopping with the first timeStamp with cumulative charge above threshold
0506           if (std::get<1>(*step) > tdcForToAOnset[waferThickness - 1] &&
0507               std::get<2>(*step) != std::get<2>(PhitRefs_bx0[id].back())) {
0508             PhitRefs_bx0[id].resize(
0509                 std::upper_bound(PhitRefs_bx0[id].begin(),
0510                                  PhitRefs_bx0[id].end(),
0511                                  hit_timeStamp(charge, 0.f, std::get<2>(*step)),
0512                                  [](const auto& i, const auto& j) { return std::get<2>(i) < std::get<2>(j); }) -
0513                 PhitRefs_bx0[id].begin());
0514             for (auto stepEnd = step + 1; stepEnd != PhitRefs_bx0[id].end(); ++stepEnd)
0515               std::get<1>(*stepEnd) += charge;
0516             break;
0517           }
0518         }
0519       }
0520     }
0521   }
0522 
0523   for (const auto& hitCollection : PhitRefs_bx0) {
0524     const uint32_t detectorId = hitCollection.first;
0525     auto simHitIt = pusimHitAccumulator_->emplace(detectorId, HGCCellHitInfo()).first;
0526 
0527     for (const auto& hit_timestamp : PhitRefs_bx0[detectorId]) {
0528       (simHitIt->second).PUhit_info[1][thisBx_].push_back(std::get<2>(hit_timestamp));
0529       (simHitIt->second).PUhit_info[0][thisBx_].push_back(std::get<0>(hit_timestamp));
0530     }
0531   }
0532 
0533   if (nchits == 0) {
0534     HGCPUSimHitDataAccumulator::iterator simHitIt = pusimHitAccumulator_->emplace(0, HGCCellHitInfo()).first;
0535     (simHitIt->second).PUhit_info[1][9].push_back(0.0);
0536     (simHitIt->second).PUhit_info[0][9].push_back(0.0);
0537   }
0538   hitRefs.clear();
0539   PhitRefs_bx0.clear();
0540 }
0541 
0542 //
0543 void HGCDigitizer::accumulate(edm::Handle<edm::PCaloHitContainer> const& hits,
0544                               int bxCrossing,
0545                               const HGCalGeometry* geom,
0546                               CLHEP::HepRandomEngine* hre) {
0547   if (nullptr == geom)
0548     return;
0549 
0550   //configuration to apply for the computation of time-of-flight
0551   auto weightToAbyEnergy = theDigitizer_->toaModeByEnergy();
0552   auto tdcForToAOnset = theDigitizer_->tdcForToAOnset();
0553   auto keV2fC = theDigitizer_->keV2fC();
0554 
0555   //create list of tuples (pos in container, RECO DetId, time) to be sorted first
0556   int nchits = (int)hits->size();
0557 
0558   std::vector<HGCCaloHitTuple_t> hitRefs;
0559   hitRefs.reserve(nchits);
0560   for (int i = 0; i < nchits; ++i) {
0561     const auto& the_hit = hits->at(i);
0562     DetId id = simToReco(geom, the_hit.id());
0563 
0564     if (verbosity_ > 0) {
0565       edm::LogVerbatim("HGCDigitizer") << "HGCDigitizer::i/p " << std::hex << the_hit.id() << " o/p " << id.rawId()
0566                                        << std::dec;
0567     }
0568 
0569     if (0 != id.rawId()) {
0570       hitRefs.emplace_back(i, id.rawId(), (float)the_hit.time());
0571     }
0572   }
0573 
0574   std::sort(hitRefs.begin(), hitRefs.end(), this->orderByDetIdThenTime);
0575   //loop over sorted hits
0576   nchits = hitRefs.size();
0577   for (int i = 0; i < nchits; ++i) {
0578     const int hitidx = std::get<0>(hitRefs[i]);
0579     const uint32_t id = std::get<1>(hitRefs[i]);
0580 
0581     //get the data for this cell, if not available then we skip it
0582 
0583     if (!validIds_.count(id))
0584       continue;
0585     HGCSimHitDataAccumulator::iterator simHitIt = simHitAccumulator_->emplace(id, HGCCellInfo()).first;
0586 
0587     if (id == 0)
0588       continue;  // to be ignored at RECO level
0589 
0590     const float toa = std::get<2>(hitRefs[i]);
0591     const PCaloHit& hit = hits->at(hitidx);
0592     const float charge = hit.energy() * 1e6 * keV2fC;
0593 
0594     //hit time: [time()]=ns + delay
0595     //accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples)
0596     const float tof = toa + tofDelay_;
0597     const int itime = std::floor(tof / bxTime_) + 9;
0598 
0599     //no need to add bx crossing - tof comes already corrected from the mixing module
0600     //itime += bxCrossing;
0601     //itime += 9;
0602 
0603     if (itime < 0 || itime > (int)maxBx_)
0604       continue;
0605 
0606     //check if time index is ok and store energy
0607     if (itime >= (int)simHitIt->second.hit_info[0].size())
0608       continue;
0609 
0610     (simHitIt->second).hit_info[0][itime] += charge;
0611 
0612     //for time-of-arrival: save the time-sorted list of timestamps with cumulative charge just above threshold
0613     //working version with pileup only for in-time hits
0614     int waferThickness = getCellThickness(geom, id);
0615     assert(waferThickness - 1 < static_cast<int>(tdcForToAOnset.size()));
0616     bool orderChanged = false;
0617     if (itime == (int)thisBx_) {
0618       //if start empty => just add charge and time
0619       if (hitRefs_bx0[id].empty()) {
0620         hitRefs_bx0[id].emplace_back(charge, tof);
0621 
0622       } else if (tof <= hitRefs_bx0[id].back().second) {
0623         //find position to insert new entry preserving time sorting
0624         std::vector<std::pair<float, float>>::iterator findPos =
0625             std::upper_bound(hitRefs_bx0[id].begin(),
0626                              hitRefs_bx0[id].end(),
0627                              std::pair<float, float>(0.f, tof),
0628                              [](const auto& i, const auto& j) { return i.second <= j.second; });
0629 
0630         std::vector<std::pair<float, float>>::iterator insertedPos = findPos;
0631         if (findPos->second == tof) {
0632           //just merge timestamps with exact timing
0633           findPos->first += charge;
0634         } else {
0635           //insert new element cumulating the charge
0636           insertedPos = hitRefs_bx0[id].insert(findPos,
0637                                                (findPos == hitRefs_bx0[id].begin())
0638                                                    ? std::pair<float, float>(charge, tof)
0639                                                    : std::pair<float, float>((findPos - 1)->first + charge, tof));
0640         }
0641 
0642         //cumulate the charge of new entry for all elements that follow in the sorted list
0643         //and resize list accounting for cases when the inserted element itself crosses the threshold
0644         for (std::vector<std::pair<float, float>>::iterator step = insertedPos; step != hitRefs_bx0[id].end(); ++step) {
0645           if (step != insertedPos)
0646             step->first += charge;
0647           // resize the list stopping with the first timeStamp with cumulative charge above threshold
0648           if (step->first > tdcForToAOnset[waferThickness - 1] && step->second != hitRefs_bx0[id].back().second) {
0649             hitRefs_bx0[id].resize(std::upper_bound(hitRefs_bx0[id].begin(),
0650                                                     hitRefs_bx0[id].end(),
0651                                                     std::pair<float, float>(0.f, step->second),
0652                                                     [](const auto& i, const auto& j) { return i.second < j.second; }) -
0653                                    hitRefs_bx0[id].begin());
0654             for (auto stepEnd = step + 1; stepEnd != hitRefs_bx0[id].end(); ++stepEnd)
0655               stepEnd->first += charge;
0656             break;
0657           }
0658         }
0659 
0660         orderChanged = true;
0661       } else {
0662         //add new entry at the end of the list
0663         if (hitRefs_bx0[id].back().first <= tdcForToAOnset[waferThickness - 1]) {
0664           hitRefs_bx0[id].emplace_back(hitRefs_bx0[id].back().first + charge, tof);
0665         }
0666       }
0667     }
0668     float accChargeForToA = hitRefs_bx0[id].empty() ? 0.f : hitRefs_bx0[id].back().first;
0669     //now compute the firing ToA through the interpolation of the consecutive time-stamps at threshold
0670     if (weightToAbyEnergy)
0671       (simHitIt->second).hit_info[1][itime] += charge * tof;
0672     else if (accChargeForToA > tdcForToAOnset[waferThickness - 1] &&
0673              ((simHitIt->second).hit_info[1][itime] == 0 || orderChanged == true)) {
0674       float fireTDC = hitRefs_bx0[id].back().second;
0675       if (hitRefs_bx0[id].size() > 1) {
0676         float chargeBeforeThr = (hitRefs_bx0[id].end() - 2)->first;
0677         float tofchargeBeforeThr = (hitRefs_bx0[id].end() - 2)->second;
0678 
0679         float deltaQ = accChargeForToA - chargeBeforeThr;
0680         float deltaTOF = fireTDC - tofchargeBeforeThr;
0681         fireTDC = (tdcForToAOnset[waferThickness - 1] - chargeBeforeThr) * deltaTOF / deltaQ + tofchargeBeforeThr;
0682       }
0683       (simHitIt->second).hit_info[1][itime] = fireTDC;
0684     }
0685   }
0686   hitRefs.clear();
0687 }
0688 void HGCDigitizer::accumulate_forPreMix(const PHGCSimAccumulator& simAccumulator, const bool minbiasFlag) {
0689   //configuration to apply for the computation of time-of-flight
0690   auto weightToAbyEnergy = theDigitizer_->toaModeByEnergy();
0691   auto tdcForToAOnset = theDigitizer_->tdcForToAOnset();
0692 
0693   if (nullptr != gHGCal_) {
0694     loadSimHitAccumulator_forPreMix(*simHitAccumulator_,
0695                                     *pusimHitAccumulator_,
0696                                     gHGCal_,
0697                                     hitRefs_bx0,
0698                                     simAccumulator,
0699                                     premixStage1MinCharge_,
0700                                     premixStage1MaxCharge_,
0701                                     !weightToAbyEnergy,
0702                                     tdcForToAOnset,
0703                                     minbiasFlag,
0704                                     hitOrder_monitor,
0705                                     thisBx_);
0706   }
0707 }
0708 
0709 //
0710 void HGCDigitizer::resetSimHitDataAccumulator() {
0711   for (HGCSimHitDataAccumulator::iterator it = simHitAccumulator_->begin(); it != simHitAccumulator_->end(); it++) {
0712     it->second.hit_info[0].fill(0.);
0713     it->second.hit_info[1].fill(0.);
0714   }
0715 }
0716 
0717 uint32_t HGCDigitizer::getType() const {
0718   uint32_t idx = std::numeric_limits<unsigned>::max();
0719   switch (theDigitizer_->det()) {
0720     case DetId::HGCalEE:
0721       idx = 0;
0722       break;
0723     case DetId::HGCalHSi:
0724       idx = 1;
0725       break;
0726     case DetId::HGCalHSc:
0727       idx = 2;
0728       break;
0729     case DetId::Forward:
0730       idx = 3;
0731       break;
0732     default:
0733       break;
0734   }
0735   return idx;
0736 }
0737 
0738 void HGCDigitizer::checkPosition(const HGCalDigiCollection* digis) const {
0739   const double tol(0.5);
0740   if (nullptr != gHGCal_) {
0741     for (const auto& digi : *(digis)) {
0742       const DetId& id = digi.id();
0743       const GlobalPoint& global = gHGCal_->getPosition(id);
0744       double r = global.perp();
0745       double z = std::abs(global.z());
0746       std::pair<double, double> zrange = gHGCal_->topology().dddConstants().rangeZ(true);
0747       std::pair<double, double> rrange = gHGCal_->topology().dddConstants().rangeR(z, true);
0748       bool ok = ((r >= rrange.first) && (r <= rrange.second) && (z >= zrange.first) && (z <= zrange.second));
0749       std::string ck = (((r < rrange.first - tol) || (r > rrange.second + tol) || (z < zrange.first - tol) ||
0750                          (z > zrange.second + tol))
0751                             ? "***** ERROR *****"
0752                             : "");
0753       bool val = gHGCal_->topology().valid(id);
0754       if ((!ok) || (!val)) {
0755         if (id.det() == DetId::HGCalEE || id.det() == DetId::HGCalHSi) {
0756           edm::LogVerbatim("HGCDigitizer") << "Check " << HGCSiliconDetId(id) << " " << global << " R " << r << ":"
0757                                            << rrange.first << ":" << rrange.second << " Z " << z << ":" << zrange.first
0758                                            << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0759         } else if (id.det() == DetId::HGCalHSc) {
0760           edm::LogVerbatim("HGCDigitizer") << "Check " << HGCScintillatorDetId(id) << " " << global << " R " << r << ":"
0761                                            << rrange.first << ":" << rrange.second << " Z " << z << ":" << zrange.first
0762                                            << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0763         } else if ((id.det() == DetId::Forward) && (id.subdetId() == static_cast<int>(HFNose))) {
0764           edm::LogVerbatim("HGCDigitizer") << "Check " << HFNoseDetId(id) << " " << global << " R " << r << ":"
0765                                            << rrange.first << ":" << rrange.second << " Z " << z << ":" << zrange.first
0766                                            << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0767         } else {
0768           edm::LogVerbatim("HGCDigitizer")
0769               << "Check " << std::hex << id.rawId() << std::dec << " " << id.det() << ":" << id.subdetId() << " "
0770               << global << " R " << r << ":" << rrange.first << ":" << rrange.second << " Z " << z << ":"
0771               << zrange.first << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0772         }
0773       }
0774     }
0775   }
0776 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!