FastTimingCommon/src/ETLDeviceSim.cc

0001 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0002 #include "SimFastTiming/FastTimingCommon/interface/ETLDeviceSim.h"
0003 #include "DataFormats/Math/interface/GeantUnits.h"
0004
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "DataFormats/ForwardDetId/interface/MTDDetId.h"
0007 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009
0010 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0011 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0012 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0013 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0014
0015 ETLDeviceSim::ETLDeviceSim(const edm::ParameterSet& pset, edm::ConsumesCollector iC)
0016     : geomToken_(iC.esConsumes()),
0017       geom_(nullptr),
0018       MIPPerMeV_(1.0 / pset.getParameter<double>("meVPerMIP")),
0019       integratedLum_(pset.getParameter<double>("IntegratedLuminosity")),
0020       fluence_(pset.getParameter<std::string>("FluenceVsRadius")),
0021       lgadGain_(pset.getParameter<std::string>("LGADGainVsFluence")),
0022       lgadGainDegradation_(pset.getParameter<std::string>("LGADGainDegradation")),
0023       applyDegradation_(pset.getParameter<bool>("applyDegradation")),
0024       bxTime_(pset.getParameter<double>("bxTime")),
0025       tofDelay_(pset.getParameter<double>("tofDelay")),
0026       MPVMuon_(pset.getParameter<std::string>("MPVMuon")),
0027       MPVPion_(pset.getParameter<std::string>("MPVPion")),
0028       MPVKaon_(pset.getParameter<std::string>("MPVKaon")),
0029       MPVElectron_(pset.getParameter<std::string>("MPVElectron")),
0030       MPVProton_(pset.getParameter<std::string>("MPVProton")) {}
0031
0032 void ETLDeviceSim::getEventSetup(const edm::EventSetup& evs) { geom_ = &evs.getData(geomToken_); }
0033
0034 void ETLDeviceSim::getHitsResponse(const std::vector<std::tuple<int, uint32_t, float> >& hitRefs,
0035                                    const edm::Handle<edm::PSimHitContainer>& hits,
0036                                    mtd_digitizer::MTDSimHitDataAccumulator* simHitAccumulator,
0037                                    CLHEP::HepRandomEngine* hre) {
0038   using namespace geant_units::operators;
0039
0040   std::vector<double> emptyV;
0041   std::vector<double> radius(1);
0042   std::vector<double> fluence(1);
0043   std::vector<double> gain(1);
0044   std::vector<double> param(2);
0045   std::vector<double> momentum(1);
0046
0047   //loop over sorted hits
0048   const int nchits = hitRefs.size();
0049   LogTrace("ETLDeviceSim") << "Processing " << nchits << " SIM hits";
0050
0051   for (int i = 0; i < nchits; ++i) {
0052     const int hitidx = std::get<0>(hitRefs[i]);
0053     const uint32_t id = std::get<1>(hitRefs[i]);
0054     const MTDDetId detId(id);
0055
0056     // Safety check (this should never happen, it should be an exception
0057     if (detId.det() != DetId::Forward || detId.mtdSubDetector() != 2) {
0058       throw cms::Exception("ETLDeviceSim")
0059           << "got a DetId that was not ETL: Det = " << detId.det() << "  subDet = " << detId.mtdSubDetector();
0060     }
0061
0062     if (id == 0)
0063       continue;  // to be ignored at RECO level
0064
0065     ETLDetId etlid(detId);
0066     DetId geoId = ETLDetId(etlid);
0067     const MTDGeomDet* thedet = geom_->idToDet(geoId);
0068     if (thedet == nullptr) {
0069       throw cms::Exception("ETLDeviceSim") << "GeographicalID: " << std::hex << geoId.rawId() << " (" << detId.rawId()
0070                                            << ") is invalid!" << std::dec << std::endl;
0071     }
0072     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0073     const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0074
0075     const float toa = std::get<2>(hitRefs[i]) + tofDelay_;
0076     const PSimHit& hit = hits->at(hitidx);
0077     float charge = convertGeVToMeV(hit.energyLoss()) * MIPPerMeV_;
0078
0079     momentum[0] = hit.pabs();
0080
0081     // particle type
0082     int particleType = abs(hit.particleType());
0083     float MPV_ = 0;
0084     if (particleType == 11) {
0085       MPV_ = MPVElectron_.evaluate(momentum, emptyV);
0086     } else if (particleType == 13) {
0087       MPV_ = MPVMuon_.evaluate(momentum, emptyV);
0088     } else if (particleType == 211) {
0089       MPV_ = MPVPion_.evaluate(momentum, emptyV);
0090     } else if (particleType == 321) {
0091       MPV_ = MPVKaon_.evaluate(momentum, emptyV);
0092     } else {
0093       MPV_ = MPVProton_.evaluate(momentum, emptyV);
0094     }
0095     float MPV_charge = convertGeVToMeV(MPV_) * MIPPerMeV_;
0096
0097     // calculate the simhit row and column
0098     const auto& position = hit.localPosition();
0099
0100     // ETL is already in module-local coordinates so just scale to cm from mm
0101     Local3DPoint simscaled(convertMmToCm(position.x()), convertMmToCm(position.y()), convertMmToCm(position.z()));
0102     const auto& global_point = thedet->toGlobal(simscaled);
0103
0104     radius[0] = global_point.perp();
0105     fluence[0] = integratedLum_ * fluence_.evaluate(radius, emptyV);
0106     gain[0] = lgadGain_.evaluate(fluence, emptyV);
0107
0108     //The following lines check whether the pixel point is actually out of the active area.
0109     if (topo.isInPixel(simscaled)) {
0110       charge *= gain[0];
0111       MPV_charge *= gain[0];
0112     } else {
0113       if (applyDegradation_) {
0114         double dGapCenter = TMath::Max(TMath::Abs(simscaled.x()), TMath::Abs(simscaled.y()));
0115         param[0] = gain[0];
0116         param[1] = dGapCenter;
0117         gain[0] = lgadGainDegradation_.evaluate(param, emptyV);
0118         charge *= gain[0];
0119         MPV_charge *= gain[0];
0120       }
0121     }
0122
0123     const auto& thepixel = topo.pixelIndex(simscaled);
0124     const uint8_t row = static_cast<uint8_t>(thepixel.first);
0125     const uint8_t col = static_cast<uint8_t>(thepixel.second);
0126     LogTrace("ETLDeviceSim") << "Processing hit in pixel # " << hitidx << " DetId " << etlid.rawId() << " row/col "
0127                              << (uint32_t)row << " " << (uint32_t)col << " inPixel " << topo.isInPixel(simscaled)
0128                              << " tof " << toa << " ene " << hit.energyLoss() << " MIP " << MIPPerMeV_ << " gain "
0129                              << gain[0] << " charge " << charge << " MPV " << MPV_charge;
0130
0131     auto simHitIt =
0132         simHitAccumulator->emplace(mtd_digitizer::MTDCellId(id, row, col), mtd_digitizer::MTDCellInfo()).first;
0133
0134     // Accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples)
0135     const int itime = std::floor(toa / bxTime_) + 9;
0136     if (itime < 0 || itime > 14)
0137       continue;
0138
0139     // Check if time index is ok and store energy
0140     if (itime >= (int)simHitIt->second.hit_info[0].size())
0141       continue;
0142     (simHitIt->second).hit_info[0][itime] += charge;
0143     // Store the time of the first SimHit in the right DataFrame bucket
0144     const float tof = toa - (itime - 9) * bxTime_;
0145
0146     if ((simHitIt->second).hit_info[1][itime] == 0. || tof < (simHitIt->second).hit_info[1][itime]) {
0147       (simHitIt->second).hit_info[1][itime] = tof;
0148     }
0149     (simHitIt->second).hit_info[2][itime] += MPV_charge;
0150   }
0151 }