File indexing completed on 2024-04-06 12:30:45
0001 #include "SimMuon/GEMDigitizer/interface/GEMSignalModel.h"
0002 #include "Geometry/GEMGeometry/interface/GEMEtaPartitionSpecs.h"
0003 #include "Geometry/CommonTopologies/interface/GEMStripTopology.h"
0004 #include "Geometry/GEMGeometry/interface/GEMGeometry.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 #include "CLHEP/Random/RandFlat.h"
0008 #include "CLHEP/Random/RandPoissonQ.h"
0009 #include "CLHEP/Random/RandGaussQ.h"
0010 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0011 #include "DataFormats/Math/interface/GeantUnits.h"
0012 #include <cmath>
0013 #include <utility>
0014 #include <map>
0015
0016 GEMSignalModel::GEMSignalModel(const edm::ParameterSet& config)
0017 : GEMDigiModel(config),
0018 averageEfficiency_(config.getParameter<double>("averageEfficiency")),
0019 averageShapingTime_(config.getParameter<double>("averageShapingTime")),
0020 timeResolution_(config.getParameter<double>("timeResolution")),
0021 timeJitter_(config.getParameter<double>("timeJitter")),
0022 signalPropagationSpeed_(config.getParameter<double>("signalPropagationSpeed")),
0023 resolutionX_(config.getParameter<double>("resolutionX")),
0024 cspeed(geant_units::operators::convertMmToCm(CLHEP::c_light)),
0025
0026 energyMinCut(28.1e-09) {}
0027
0028 GEMSignalModel::~GEMSignalModel() {}
0029
0030 void GEMSignalModel::simulate(const GEMEtaPartition* roll,
0031 const edm::PSimHitContainer& simHits,
0032 CLHEP::HepRandomEngine* engine,
0033 Strips& strips_,
0034 DetectorHitMap& detectorHitMap_) {
0035 const GEMStripTopology* top(dynamic_cast<const GEMStripTopology*>(&(roll->topology())));
0036 for (const auto& hit : simHits) {
0037 if (hit.energyLoss() < energyMinCut)
0038 continue;
0039 const int bx(getSimHitBx(&hit, engine));
0040 const std::vector<std::pair<int, int> >& cluster(simulateClustering(top, &hit, bx, engine));
0041 for (const auto& digi : cluster) {
0042 detectorHitMap_.emplace(digi, &hit);
0043 strips_.emplace(digi);
0044 }
0045 }
0046 }
0047
0048 int GEMSignalModel::getSimHitBx(const PSimHit* simhit, CLHEP::HepRandomEngine* engine) {
0049 int bx = -999;
0050 const LocalPoint simHitPos(simhit->localPosition());
0051 const float tof(simhit->timeOfFlight());
0052
0053 float randomJitterTime = CLHEP::RandGaussQ::shoot(engine, 0., timeJitter_);
0054 const GEMDetId id(simhit->detUnitId());
0055 const GEMEtaPartition* roll(geometry_->etaPartition(id));
0056 if (!roll) {
0057 throw cms::Exception("Geometry") << "GEMSignalModel::getSimHitBx() - GEM simhit id does not match any GEM roll id: "
0058 << id << "\n";
0059 return 999;
0060 }
0061 if (roll->id().region() == 0) {
0062 throw cms::Exception("Geometry")
0063 << "GEMSignalModel::getSimHitBx() - this GEM id is from barrel, which cannot happen: " << roll->id() << "\n";
0064 }
0065 const int nstrips = roll->nstrips();
0066 float middleStrip = nstrips / 2.;
0067 const LocalPoint& middleOfRoll = roll->centreOfStrip(middleStrip);
0068 const GlobalPoint& globMiddleRol = roll->toGlobal(middleOfRoll);
0069 double muRadius = sqrt(globMiddleRol.x() * globMiddleRol.x() + globMiddleRol.y() * globMiddleRol.y() +
0070 globMiddleRol.z() * globMiddleRol.z());
0071 double timeCalibrationOffset_ = muRadius / cspeed;
0072
0073 const GEMStripTopology* top(dynamic_cast<const GEMStripTopology*>(&(roll->topology())));
0074 const float halfStripLength(0.5 * top->stripLength());
0075 const float distanceFromEdge(halfStripLength - simHitPos.y());
0076
0077
0078 double signalPropagationSpeedTrue = signalPropagationSpeed_ * cspeed;
0079
0080
0081 const float averagePropagationTime(distanceFromEdge / signalPropagationSpeedTrue);
0082
0083 float randomResolutionTime = CLHEP::RandGaussQ::shoot(engine, 0., timeResolution_);
0084 const float simhitTime(tof + averageShapingTime_ + randomResolutionTime + averagePropagationTime + randomJitterTime);
0085 float referenceTime = 0.;
0086 referenceTime = timeCalibrationOffset_ + halfStripLength / signalPropagationSpeedTrue + averageShapingTime_;
0087 const float timeDifference(simhitTime - referenceTime);
0088
0089 bx = static_cast<int>(std::round((timeDifference) / 25.));
0090
0091
0092 LogDebug("GEMDigiProducer") << "checktime "
0093 << "bx = " << bx << "\tdeltaT = " << timeDifference << "\tsimT = " << simhitTime
0094 << "\trefT = " << referenceTime << "\ttof = " << tof
0095 << "\tavePropT = " << averagePropagationTime
0096 << "\taveRefPropT = " << halfStripLength / signalPropagationSpeedTrue << "\n";
0097 return bx;
0098 }
0099
0100 std::vector<std::pair<int, int> > GEMSignalModel::simulateClustering(const GEMStripTopology* top,
0101 const PSimHit* simHit,
0102 const int bx,
0103 CLHEP::HepRandomEngine* engine) {
0104 const LocalPoint& hit_entry(simHit->entryPoint());
0105 const LocalPoint& hit_exit(simHit->exitPoint());
0106
0107 LocalPoint start_point, end_point;
0108 if (hit_entry.x() < hit_exit.x()) {
0109 start_point = hit_entry;
0110 end_point = hit_exit;
0111 } else {
0112 start_point = hit_exit;
0113 end_point = hit_entry;
0114 }
0115
0116
0117 float smeared_start_x = start_point.x() - std::abs(CLHEP::RandGaussQ::shoot(engine, 0, resolutionX_));
0118 float smeared_end_x = end_point.x() + std::abs(CLHEP::RandGaussQ::shoot(engine, 0, resolutionX_));
0119
0120 LocalPoint smeared_start_point(smeared_start_x, start_point.y(), start_point.z());
0121 LocalPoint smeared_end_point(smeared_end_x, end_point.y(), end_point.z());
0122
0123 int cluster_start = top->channel(smeared_start_point);
0124 int cluster_end = top->channel(smeared_end_point);
0125
0126 std::vector<std::pair<int, int> > cluster;
0127 for (int strip = cluster_start; strip <= cluster_end; strip++) {
0128 cluster.emplace_back(strip, bx);
0129 }
0130
0131 return cluster;
0132 }