FastTimingCommon/interface/MTDDigitizer.h

0001 #ifndef FastTimingSimProducers_FastTimingCommon_MTDDigitizer_h
0002 #define FastTimingSimProducers_FastTimingCommon_MTDDigitizer_h
0003
0004 #include "SimFastTiming/FastTimingCommon/interface/MTDDigitizerBase.h"
0005 #include "SimFastTiming/FastTimingCommon/interface/MTDDigitizerTraits.h"
0006
0007 #include "DataFormats/DetId/interface/DetId.h"
0008 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0009 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011
0012 #include "FWCore/Framework/interface/ConsumesCollector.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/ProducesCollector.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
0018 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
0019 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0020 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0021 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0022 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0023
0024 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0025
0026 #include "DataFormats/Math/interface/liblogintpack.h"
0027
0028 #include <vector>
0029 #include <unordered_map>
0030 #include <unordered_set>
0031 #include <memory>
0032 #include <tuple>
0033
0034 namespace mtd_digitizer {
0035
0036   namespace MTDHelpers {
0037     // index , det id, time
0038     typedef std::tuple<int, uint32_t, float> MTDCaloHitTuple_t;
0039
0040     inline bool orderByDetIdThenTime(const MTDCaloHitTuple_t& a, const MTDCaloHitTuple_t& b) {
0041       unsigned int detId_a(std::get<1>(a)), detId_b(std::get<1>(b));
0042
0043       if (detId_a < detId_b)
0044         return true;
0045       if (detId_a > detId_b)
0046         return false;
0047
0048       double time_a(std::get<2>(a)), time_b(std::get<2>(b));
0049       if (time_a < time_b)
0050         return true;
0051
0052       return false;
0053     }
0054   }  // namespace MTDHelpers
0055
0056   inline void saveSimHitAccumulator(PMTDSimAccumulator& simResult,
0057                                     const MTDSimHitDataAccumulator& simData,
0058                                     const float minCharge,
0059                                     const float maxCharge) {
0060     constexpr auto nEnergies = std::tuple_size<decltype(MTDCellInfo().hit_info)>::value;
0061     static_assert(nEnergies == PMTDSimAccumulator::Data::energyMask + 1,
0062                   "PMTDSimAccumulator bit pattern needs to be updated");
0063     static_assert(nSamples == PMTDSimAccumulator::Data::sampleMask,
0064                   "PMTDSimAccumulator bit pattern needs to be updated");
0065
0066     const float minPackChargeLog = minCharge > 0.f ? std::log(minCharge) : -2;
0067     const float maxPackChargeLog = std::log(maxCharge);
0068     constexpr uint16_t base = PMTDSimAccumulator::Data::dataMask;
0069
0070     simResult.reserve(simData.size());
0071     // mimicking the digitization
0072     for (const auto& elem : simData) {
0073       // store only non-zero
0074       for (size_t iEn = 0; iEn < nEnergies; ++iEn) {
0075         const auto& samples = elem.second.hit_info[iEn];
0076         for (size_t iSample = 0; iSample < nSamples; ++iSample) {
0077           if (samples[iSample] > minCharge) {
0078             unsigned short packed;
0079             if (iEn == 1 || iEn == 3) {
0080               // assuming linear range for tof of 0..25
0081               packed = samples[iSample] / PREMIX_MAX_TOF * base;
0082             } else {
0083               packed = logintpack::pack16log(samples[iSample], minPackChargeLog, maxPackChargeLog, base);
0084             }
0085             simResult.emplace_back(elem.first.detid_, elem.first.row_, elem.first.column_, iEn, iSample, packed);
0086           }
0087         }
0088       }
0089     }
0090   }
0091
0092   inline void loadSimHitAccumulator(MTDSimHitDataAccumulator& simData,
0093                                     const PMTDSimAccumulator& simAccumulator,
0094                                     const float minCharge,
0095                                     const float maxCharge) {
0096     const float minPackChargeLog = minCharge > 0.f ? std::log(minCharge) : -2;
0097     const float maxPackChargeLog = std::log(maxCharge);
0098     constexpr uint16_t base = PMTDSimAccumulator::Data::dataMask;
0099
0100     for (const auto& detIdIndexHitInfo : simAccumulator) {
0101       auto foo = simData.emplace(
0102           MTDCellId(detIdIndexHitInfo.detId(), detIdIndexHitInfo.row(), detIdIndexHitInfo.column()), MTDCellInfo());
0103       auto simIt = foo.first;
0104       auto& hit_info = simIt->second.hit_info;
0105
0106       size_t iEn = detIdIndexHitInfo.energyIndex();
0107       size_t iSample = detIdIndexHitInfo.sampleIndex();
0108
0109       if (iEn > PMTDSimAccumulator::Data::energyMask + 1 || iSample > PMTDSimAccumulator::Data::sampleMask)
0110         throw cms::Exception("MTDDigitixer::loadSimHitAccumulator")
0111             << "Index out of range: iEn = " << iEn << " iSample = " << iSample << std::endl;
0112
0113       float value;
0114       if (iEn == 1 || iEn == 3) {
0115         value = static_cast<float>(detIdIndexHitInfo.data()) / base * PREMIX_MAX_TOF;
0116       } else {
0117         value = logintpack::unpack16log(detIdIndexHitInfo.data(), minPackChargeLog, maxPackChargeLog, base);
0118       }
0119
0120       if (iEn == 0 || iEn == 2) {
0121         hit_info[iEn][iSample] += value;
0122       } else if (hit_info[iEn][iSample] == 0 || value < hit_info[iEn][iSample]) {
0123         // For iEn==1 the digitizers just set the TOF of the first SimHit
0124         hit_info[iEn][iSample] = value;
0125       }
0126     }
0127   }
0128
0129   template <class Traits>
0130   class MTDDigitizer : public MTDDigitizerBase {
0131   public:
0132     typedef typename Traits::DeviceSim DeviceSim;
0133     typedef typename Traits::ElectronicsSim ElectronicsSim;
0134     typedef typename Traits::DigiCollection DigiCollection;
0135
0136     MTDDigitizer(const edm::ParameterSet& config, edm::ProducesCollector producesCollector, edm::ConsumesCollector& iC)
0137         : MTDDigitizerBase(config, producesCollector, iC),
0138           geomToken_(iC.esConsumes()),
0139           geom_(nullptr),
0140           deviceSim_(config.getParameterSet("DeviceSimulation"), iC),
0141           electronicsSim_(config.getParameterSet("ElectronicsSimulation"), iC),
0142           maxSimHitsAccTime_(config.getParameter<uint32_t>("maxSimHitsAccTime")) {}
0143
0144     ~MTDDigitizer() override {}
0145
0146     /**
0147        @short handle SimHit accumulation
0148     */
0149     void accumulate(edm::Event const& e, edm::EventSetup const& c, CLHEP::HepRandomEngine* hre) override;
0150     void accumulate(PileUpEventPrincipal const& e, edm::EventSetup const& c, CLHEP::HepRandomEngine* hre) override;
0151     void accumulate(edm::Handle<edm::PSimHitContainer> const& hits,
0152                     int bxCrossing,
0153                     CLHEP::HepRandomEngine* hre) override;
0154     // for premixing
0155     void accumulate(const PMTDSimAccumulator& simAccumulator) override;
0156
0157     /**
0158        @short actions at the start/end of event
0159     */
0160     void initializeEvent(edm::Event const& e, edm::EventSetup const& c) override;
0161     void finalizeEvent(edm::Event& e, edm::EventSetup const& c, CLHEP::HepRandomEngine* hre) override;
0162
0163   private:
0164     void resetSimHitDataAccumulator() { MTDSimHitDataAccumulator().swap(simHitAccumulator_); }
0165
0166     const edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> geomToken_;
0167     const MTDGeometry* geom_;
0168
0169     // implementations
0170     DeviceSim deviceSim_;            // processes a given simhit into an entry in a MTDSimHitDataAccumulator
0171     ElectronicsSim electronicsSim_;  // processes a MTDSimHitDataAccumulator into a BTLDigiCollection/ETLDigiCollection
0172
0173     //handle sim hits
0174     const int maxSimHitsAccTime_;
0175     MTDSimHitDataAccumulator simHitAccumulator_;
0176   };
0177
0178   template <class Traits>
0179   void MTDDigitizer<Traits>::accumulate(edm::Event const& e, edm::EventSetup const& c, CLHEP::HepRandomEngine* hre) {
0180     edm::Handle<edm::PSimHitContainer> simHits;
0181     e.getByLabel(inputSimHits_, simHits);
0182     if (simHits.isValid())
0183       accumulate(simHits, 0, hre);
0184     else
0185       edm::LogWarning("MTDMix") << "MTDDigitizer:: Cannot find hits for " << inputSimHits_;
0186   }
0187
0188   template <class Traits>
0189   void MTDDigitizer<Traits>::accumulate(PileUpEventPrincipal const& e,
0190                                         edm::EventSetup const& c,
0191                                         CLHEP::HepRandomEngine* hre) {
0192     edm::Handle<edm::PSimHitContainer> simHits;
0193     e.getByLabel(inputSimHits_, simHits);
0194     if (simHits.isValid())
0195       accumulate(simHits, e.bunchCrossing(), hre);
0196     else
0197       edm::LogWarning("MTDMix") << "MTDDigitizer:: Cannot find hits for " << inputSimHits_;
0198   }
0199
0200   template <class Traits>
0201   void MTDDigitizer<Traits>::accumulate(edm::Handle<edm::PSimHitContainer> const& hits,
0202                                         int bxCrossing,
0203                                         CLHEP::HepRandomEngine* hre) {
0204     using namespace MTDHelpers;
0205
0206     //create list of tuples (pos in container, RECO DetId, time) to be sorted first
0207     int nchits = (int)hits->size();
0208     std::vector<MTDCaloHitTuple_t> hitRefs;
0209     hitRefs.reserve(nchits);
0210     for (int i = 0; i < nchits; ++i) {
0211       const auto& the_hit = hits->at(i);
0212
0213       DetId id = the_hit.detUnitId();
0214
0215       if (verbosity_ > 0) {
0216         edm::LogInfo("MTDDigitizer") << " i/p " << std::hex << the_hit.detUnitId() << std::dec << " o/p " << id.rawId()
0217                                      << std::endl;
0218       }
0219
0220       if (0 != id.rawId()) {
0221         hitRefs.emplace_back(i, id.rawId(), the_hit.tof());
0222       }
0223     }
0224     std::sort(hitRefs.begin(), hitRefs.end(), MTDHelpers::orderByDetIdThenTime);
0225
0226     deviceSim_.getHitsResponse(hitRefs, hits, &simHitAccumulator_, hre);
0227
0228     hitRefs.clear();
0229   }
0230
0231   template <class Traits>
0232   void MTDDigitizer<Traits>::accumulate(const PMTDSimAccumulator& simAccumulator) {
0233     loadSimHitAccumulator(simHitAccumulator_, simAccumulator, premixStage1MinCharge_, premixStage1MaxCharge_);
0234   }
0235
0236   template <class Traits>
0237   void MTDDigitizer<Traits>::initializeEvent(edm::Event const& e, edm::EventSetup const& c) {
0238     geom_ = &c.getData(geomToken_);
0239
0240     deviceSim_.getEvent(e);
0241     deviceSim_.getEventSetup(c);
0242     if (not premixStage1_) {
0243       electronicsSim_.getEvent(e);
0244       electronicsSim_.getEventSetup(c);
0245     }
0246   }
0247
0248   template <class Traits>
0249   void MTDDigitizer<Traits>::finalizeEvent(edm::Event& e, edm::EventSetup const& c, CLHEP::HepRandomEngine* hre) {
0250     if (premixStage1_) {
0251       auto simResult = std::make_unique<PMTDSimAccumulator>();
0252       saveSimHitAccumulator(*simResult, simHitAccumulator_, premixStage1MinCharge_, premixStage1MaxCharge_);
0253       e.put(std::move(simResult), digiCollection_);
0254     } else {
0255       auto digiCollection = std::make_unique<DigiCollection>();
0256       electronicsSim_.run(simHitAccumulator_, *digiCollection, hre);
0257       e.put(std::move(digiCollection), digiCollection_);
0258     }
0259
0260     //release memory for next event
0261     resetSimHitDataAccumulator();
0262   }
0263 }  // namespace mtd_digitizer
0264
0265 #endif