Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-23 02:42:51 

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   int count_thisbx = 0;
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       ++count_thisbx;
0472       if (PhitRefs_bx0[id].empty()) {
0473         PhitRefs_bx0[id].emplace_back(charge, charge, tof);
0474       } else if (tof > std::get<2>(PhitRefs_bx0[id].back())) {
0475         PhitRefs_bx0[id].emplace_back(charge, charge + std::get<1>(PhitRefs_bx0[id].back()), tof);
0476       } else if (tof == std::get<2>(PhitRefs_bx0[id].back())) {
0477         std::get<0>(PhitRefs_bx0[id].back()) += charge;
0478         std::get<1>(PhitRefs_bx0[id].back()) += charge;
0479       } else {
0480         //find position to insert new entry preserving time sorting
0481         auto findPos = std::upper_bound(PhitRefs_bx0[id].begin(),
0482                                         PhitRefs_bx0[id].end(),
0483                                         hit_timeStamp(charge, 0.f, tof),
0484                                         [](const auto& i, const auto& j) { return std::get<2>(i) <= std::get<2>(j); });
0485 
0486         auto insertedPos = findPos;
0487 
0488         if (tof == std::get<2>(*(findPos - 1))) {
0489           std::get<0>(*(findPos - 1)) += charge;
0490           std::get<1>(*(findPos - 1)) += charge;
0491 
0492         } else {
0493           insertedPos = PhitRefs_bx0[id].insert(findPos,
0494                                                 (findPos == PhitRefs_bx0[id].begin())
0495                                                     ? hit_timeStamp(charge, charge, tof)
0496                                                     : hit_timeStamp(charge, charge + std::get<1>(*(findPos - 1)), tof));
0497         }
0498         //cumulate the charge of new entry for all elements that follow in the sorted list
0499         //and resize list accounting for cases when the inserted element itself crosses the threshold
0500 
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     bool orderChanged = false;
0616     if (itime == (int)thisBx_) {
0617       //if start empty => just add charge and time
0618       if (hitRefs_bx0[id].empty()) {
0619         hitRefs_bx0[id].emplace_back(charge, tof);
0620 
0621       } else if (tof <= hitRefs_bx0[id].back().second) {
0622         //find position to insert new entry preserving time sorting
0623         std::vector<std::pair<float, float>>::iterator findPos =
0624             std::upper_bound(hitRefs_bx0[id].begin(),
0625                              hitRefs_bx0[id].end(),
0626                              std::pair<float, float>(0.f, tof),
0627                              [](const auto& i, const auto& j) { return i.second <= j.second; });
0628 
0629         std::vector<std::pair<float, float>>::iterator insertedPos = findPos;
0630         if (findPos->second == tof) {
0631           //just merge timestamps with exact timing
0632           findPos->first += charge;
0633         } else {
0634           //insert new element cumulating the charge
0635           insertedPos = hitRefs_bx0[id].insert(findPos,
0636                                                (findPos == hitRefs_bx0[id].begin())
0637                                                    ? std::pair<float, float>(charge, tof)
0638                                                    : std::pair<float, float>((findPos - 1)->first + charge, tof));
0639         }
0640 
0641         //cumulate the charge of new entry for all elements that follow in the sorted list
0642         //and resize list accounting for cases when the inserted element itself crosses the threshold
0643         for (std::vector<std::pair<float, float>>::iterator step = insertedPos; step != hitRefs_bx0[id].end(); ++step) {
0644           if (step != insertedPos)
0645             step->first += charge;
0646           // resize the list stopping with the first timeStamp with cumulative charge above threshold
0647           if (step->first > tdcForToAOnset[waferThickness - 1] && step->second != hitRefs_bx0[id].back().second) {
0648             hitRefs_bx0[id].resize(std::upper_bound(hitRefs_bx0[id].begin(),
0649                                                     hitRefs_bx0[id].end(),
0650                                                     std::pair<float, float>(0.f, step->second),
0651                                                     [](const auto& i, const auto& j) { return i.second < j.second; }) -
0652                                    hitRefs_bx0[id].begin());
0653             for (auto stepEnd = step + 1; stepEnd != hitRefs_bx0[id].end(); ++stepEnd)
0654               stepEnd->first += charge;
0655             break;
0656           }
0657         }
0658 
0659         orderChanged = true;
0660       } else {
0661         //add new entry at the end of the list
0662         if (hitRefs_bx0[id].back().first <= tdcForToAOnset[waferThickness - 1]) {
0663           hitRefs_bx0[id].emplace_back(hitRefs_bx0[id].back().first + charge, tof);
0664         }
0665       }
0666     }
0667     float accChargeForToA = hitRefs_bx0[id].empty() ? 0.f : hitRefs_bx0[id].back().first;
0668     //now compute the firing ToA through the interpolation of the consecutive time-stamps at threshold
0669     if (weightToAbyEnergy)
0670       (simHitIt->second).hit_info[1][itime] += charge * tof;
0671     else if (accChargeForToA > tdcForToAOnset[waferThickness - 1] &&
0672              ((simHitIt->second).hit_info[1][itime] == 0 || orderChanged == true)) {
0673       float fireTDC = hitRefs_bx0[id].back().second;
0674       if (hitRefs_bx0[id].size() > 1) {
0675         float chargeBeforeThr = (hitRefs_bx0[id].end() - 2)->first;
0676         float tofchargeBeforeThr = (hitRefs_bx0[id].end() - 2)->second;
0677 
0678         float deltaQ = accChargeForToA - chargeBeforeThr;
0679         float deltaTOF = fireTDC - tofchargeBeforeThr;
0680         fireTDC = (tdcForToAOnset[waferThickness - 1] - chargeBeforeThr) * deltaTOF / deltaQ + tofchargeBeforeThr;
0681       }
0682       (simHitIt->second).hit_info[1][itime] = fireTDC;
0683     }
0684   }
0685   hitRefs.clear();
0686 }
0687 void HGCDigitizer::accumulate_forPreMix(const PHGCSimAccumulator& simAccumulator, const bool minbiasFlag) {
0688   //configuration to apply for the computation of time-of-flight
0689   auto weightToAbyEnergy = theDigitizer_->toaModeByEnergy();
0690   auto tdcForToAOnset = theDigitizer_->tdcForToAOnset();
0691 
0692   if (nullptr != gHGCal_) {
0693     loadSimHitAccumulator_forPreMix(*simHitAccumulator_,
0694                                     *pusimHitAccumulator_,
0695                                     gHGCal_,
0696                                     hitRefs_bx0,
0697                                     simAccumulator,
0698                                     premixStage1MinCharge_,
0699                                     premixStage1MaxCharge_,
0700                                     !weightToAbyEnergy,
0701                                     tdcForToAOnset,
0702                                     minbiasFlag,
0703                                     hitOrder_monitor,
0704                                     thisBx_);
0705   }
0706 }
0707 
0708 //
0709 void HGCDigitizer::resetSimHitDataAccumulator() {
0710   for (HGCSimHitDataAccumulator::iterator it = simHitAccumulator_->begin(); it != simHitAccumulator_->end(); it++) {
0711     it->second.hit_info[0].fill(0.);
0712     it->second.hit_info[1].fill(0.);
0713   }
0714 }
0715 
0716 uint32_t HGCDigitizer::getType() const {
0717   uint32_t idx = std::numeric_limits<unsigned>::max();
0718   switch (theDigitizer_->det()) {
0719     case DetId::HGCalEE:
0720       idx = 0;
0721       break;
0722     case DetId::HGCalHSi:
0723       idx = 1;
0724       break;
0725     case DetId::HGCalHSc:
0726       idx = 2;
0727       break;
0728     case DetId::Forward:
0729       idx = 3;
0730       break;
0731     default:
0732       break;
0733   }
0734   return idx;
0735 }
0736 
0737 void HGCDigitizer::checkPosition(const HGCalDigiCollection* digis) const {
0738   const double tol(0.5);
0739   if (nullptr != gHGCal_) {
0740     for (const auto& digi : *(digis)) {
0741       const DetId& id = digi.id();
0742       const GlobalPoint& global = gHGCal_->getPosition(id);
0743       double r = global.perp();
0744       double z = std::abs(global.z());
0745       std::pair<double, double> zrange = gHGCal_->topology().dddConstants().rangeZ(true);
0746       std::pair<double, double> rrange = gHGCal_->topology().dddConstants().rangeR(z, true);
0747       bool ok = ((r >= rrange.first) && (r <= rrange.second) && (z >= zrange.first) && (z <= zrange.second));
0748       std::string ck = (((r < rrange.first - tol) || (r > rrange.second + tol) || (z < zrange.first - tol) ||
0749                          (z > zrange.second + tol))
0750                             ? "***** ERROR *****"
0751                             : "");
0752       bool val = gHGCal_->topology().valid(id);
0753       if ((!ok) || (!val)) {
0754         if (id.det() == DetId::HGCalEE || id.det() == DetId::HGCalHSi) {
0755           edm::LogVerbatim("HGCDigitizer") << "Check " << HGCSiliconDetId(id) << " " << global << " R " << r << ":"
0756                                            << rrange.first << ":" << rrange.second << " Z " << z << ":" << zrange.first
0757                                            << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0758         } else if (id.det() == DetId::HGCalHSc) {
0759           edm::LogVerbatim("HGCDigitizer") << "Check " << HGCScintillatorDetId(id) << " " << global << " R " << r << ":"
0760                                            << rrange.first << ":" << rrange.second << " Z " << z << ":" << zrange.first
0761                                            << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0762         } else if ((id.det() == DetId::Forward) && (id.subdetId() == static_cast<int>(HFNose))) {
0763           edm::LogVerbatim("HGCDigitizer") << "Check " << HFNoseDetId(id) << " " << global << " R " << r << ":"
0764                                            << rrange.first << ":" << rrange.second << " Z " << z << ":" << zrange.first
0765                                            << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0766         } else {
0767           edm::LogVerbatim("HGCDigitizer")
0768               << "Check " << std::hex << id.rawId() << std::dec << " " << id.det() << ":" << id.subdetId() << " "
0769               << global << " R " << r << ":" << rrange.first << ":" << rrange.second << " Z " << z << ":"
0770               << zrange.first << ":" << zrange.second << " Flag " << ok << ":" << val << " " << ck;
0771         }
0772       }
0773     }
0774   }
0775 }
[ 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!