FastTimingCommon/src/BTLDeviceSim.cc

0001 #include "SimFastTiming/FastTimingCommon/interface/BTLDeviceSim.h"
0002 #include "DataFormats/Math/interface/GeantUnits.h"
0003 #include "DataFormats/DetId/interface/DetId.h"
0004 #include "DataFormats/ForwardDetId/interface/MTDDetId.h"
0005 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0006 #include "FWCore/Framework/interface/ConsumesCollector.h"
0007
0008 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0009 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0010 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0011 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0012
0013 #include "CLHEP/Random/RandGaussQ.h"
0014
0015 BTLDeviceSim::BTLDeviceSim(const edm::ParameterSet& pset, edm::ConsumesCollector iC)
0016     : geomToken_(iC.esConsumes()),
0017       topoToken_(iC.esConsumes()),
0018       geom_(nullptr),
0019       topo_(nullptr),
0020       bxTime_(pset.getParameter<double>("bxTime")),
0021       LightYield_(pset.getParameter<double>("LightYield")),
0022       LightCollEff_(pset.getParameter<double>("LightCollectionEff")),
0023       LightCollSlope_(pset.getParameter<double>("LightCollectionSlope")),
0024       PDE_(pset.getParameter<double>("PhotonDetectionEff")),
0025       LCEpositionSlope_(pset.getParameter<double>("LCEpositionSlope")) {}
0026
0027 void BTLDeviceSim::getEventSetup(const edm::EventSetup& evs) {
0028   geom_ = &evs.getData(geomToken_);
0029   topo_ = &evs.getData(topoToken_);
0030 }
0031
0032 void BTLDeviceSim::getHitsResponse(const std::vector<std::tuple<int, uint32_t, float> >& hitRefs,
0033                                    const edm::Handle<edm::PSimHitContainer>& hits,
0034                                    mtd_digitizer::MTDSimHitDataAccumulator* simHitAccumulator,
0035                                    CLHEP::HepRandomEngine* hre) {
0036   using namespace geant_units::operators;
0037
0038   //loop over sorted simHits
0039   for (auto const& hitRef : hitRefs) {
0040     const int hitidx = std::get<0>(hitRef);
0041     const uint32_t id = std::get<1>(hitRef);
0042     const MTDDetId detId(id);
0043     const PSimHit& hit = hits->at(hitidx);
0044
0045     // --- Safety check on the detector ID
0046     if (detId.det() != DetId::Forward || detId.mtdSubDetector() != 1)
0047       continue;
0048
0049     if (id == 0)
0050       continue;  // to be ignored at RECO level
0051
0052     BTLDetId btlid(detId);
0053     DetId geoId = btlid.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo_->getMTDTopologyMode()));
0054     const MTDGeomDet* thedet = geom_->idToDet(geoId);
0055
0056     if (thedet == nullptr) {
0057       throw cms::Exception("BTLDeviceSim") << "GeographicalID: " << std::hex << geoId.rawId() << " (" << detId.rawId()
0058                                            << ") is invalid!" << std::dec << std::endl;
0059     }
0060     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0061     const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0062     // calculate the simhit row and column
0063     const auto& position = hit.localPosition();
0064     Local3DPoint simscaled(convertMmToCm(position.x()), convertMmToCm(position.y()), convertMmToCm(position.z()));
0065     // translate from crystal-local coordinates to module-local coordinates to get the row and column
0066     simscaled = topo.pixelToModuleLocalPoint(simscaled, btlid.row(topo.nrows()), btlid.column(topo.nrows()));
0067     const auto& thepixel = topo.pixel(simscaled);
0068     uint8_t row(thepixel.first), col(thepixel.second);
0069
0070     if (btlid.row(topo.nrows()) != row || btlid.column(topo.nrows()) != col) {
0071       edm::LogWarning("BTLDeviceSim") << "BTLDetId (row,column): (" << btlid.row(topo.nrows()) << ','
0072                                       << btlid.column(topo.nrows()) << ") is not equal to "
0073                                       << "topology (row,column): (" << uint32_t(row) << ',' << uint32_t(col)
0074                                       << "), overriding to detid";
0075       row = btlid.row(topo.nrows());
0076       col = btlid.column(topo.nrows());
0077     }
0078
0079     // --- Store the detector element ID as a key of the MTDSimHitDataAccumulator map
0080     auto simHitIt =
0081         simHitAccumulator->emplace(mtd_digitizer::MTDCellId(id, row, col), mtd_digitizer::MTDCellInfo()).first;
0082
0083     // --- Get the simHit energy and convert it from MeV to photo-electrons
0084     float Npe = convertGeVToMeV(hit.energyLoss()) * LightYield_ * LightCollEff_ * PDE_;
0085
0086     // --- Calculate the light propagation time to the crystal bases (labeled L and R)
0087     double distR = 0.5 * topo.pitch().first - convertMmToCm(hit.localPosition().x());
0088     double distL = 0.5 * topo.pitch().first + convertMmToCm(hit.localPosition().x());
0089
0090     double tR = std::get<2>(hitRef) + LightCollSlope_ * distR;
0091     double tL = std::get<2>(hitRef) + LightCollSlope_ * distL;
0092
0093     // --- Accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples)
0094     const int iBXR = std::floor(tR / bxTime_) + mtd_digitizer::kInTimeBX;
0095     const int iBXL = std::floor(tL / bxTime_) + mtd_digitizer::kInTimeBX;
0096
0097     // --- Right side
0098     if (iBXR > 0 && iBXR < mtd_digitizer::kNumberOfBX) {
0099       // Accumulate the energy of simHits in the same crystal
0100       (simHitIt->second).hit_info[0][iBXR] += Npe * (1. + LCEpositionSlope_ * convertMmToCm(hit.localPosition().x()));
0101
0102       // Store the time of the first SimHit in the i-th BX
0103       if ((simHitIt->second).hit_info[1][iBXR] == 0 || tR < (simHitIt->second).hit_info[1][iBXR])
0104         (simHitIt->second).hit_info[1][iBXR] = tR - (iBXR - mtd_digitizer::kInTimeBX) * bxTime_;
0105     }
0106
0107     // --- Left side
0108     if (iBXL > 0 && iBXL < mtd_digitizer::kNumberOfBX) {
0109       // Accumulate the energy of simHits in the same crystal
0110       (simHitIt->second).hit_info[2][iBXL] += Npe * (1. - LCEpositionSlope_ * convertMmToCm(hit.localPosition().x()));
0111
0112       // Store the time of the first SimHit in the i-th BX
0113       if ((simHitIt->second).hit_info[3][iBXL] == 0 || tL < (simHitIt->second).hit_info[3][iBXL])
0114         (simHitIt->second).hit_info[3][iBXL] = tL - (iBXL - mtd_digitizer::kInTimeBX) * bxTime_;
0115     }
0116
0117   }  // hitRef loop
0118 }