GEMDigitizer/src/GEMSignalModel.cc

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       bx0filter_(config.getParameter<bool>("bx0filter")),
0024       resolutionX_(config.getParameter<double>("resolutionX")),
0025       cspeed(geant_units::operators::convertMmToCm(CLHEP::c_light)),
0026       // average energy required to remove an electron due to ionization for an Ar/CO2 gas mixture (in the ratio of 70/30) is 28.1 eV
0027       energyMinCut(28.1e-09) {}
0028
0029 GEMSignalModel::~GEMSignalModel() {}
0030
0031 void GEMSignalModel::simulate(const GEMEtaPartition* roll,
0032                               const edm::PSimHitContainer& simHits,
0033                               CLHEP::HepRandomEngine* engine,
0034                               Strips& strips_,
0035                               DetectorHitMap& detectorHitMap_) {
0036   const GEMStripTopology* top(dynamic_cast<const GEMStripTopology*>(&(roll->topology())));
0037   for (const auto& hit : simHits) {
0038     if (hit.energyLoss() < energyMinCut)
0039       continue;
0040     const int bx(getSimHitBx(&hit, engine));
0041     if (bx != 0 and bx0filter_)
0042       continue;
0043     const std::vector<std::pair<int, int> >& cluster(simulateClustering(top, &hit, bx, engine));
0044     for (const auto& digi : cluster) {
0045       detectorHitMap_.emplace(digi, &hit);
0046       strips_.emplace(digi);
0047     }
0048   }
0049 }
0050
0051 int GEMSignalModel::getSimHitBx(const PSimHit* simhit, CLHEP::HepRandomEngine* engine) {
0052   int bx = -999;
0053   const LocalPoint simHitPos(simhit->localPosition());
0054   const float tof(simhit->timeOfFlight());
0055   // random Gaussian time correction due to electronics jitter
0056   float randomJitterTime = CLHEP::RandGaussQ::shoot(engine, 0., timeJitter_);
0057   const GEMDetId id(simhit->detUnitId());
0058   const GEMEtaPartition* roll(geometry_->etaPartition(id));
0059   if (!roll) {
0060     throw cms::Exception("Geometry") << "GEMSignalModel::getSimHitBx() - GEM simhit id does not match any GEM roll id: "
0061                                      << id << "\n";
0062     return 999;
0063   }
0064   if (roll->id().region() == 0) {
0065     throw cms::Exception("Geometry")
0066         << "GEMSignalModel::getSimHitBx() - this GEM id is from barrel, which cannot happen: " << roll->id() << "\n";
0067   }
0068   const int nstrips = roll->nstrips();
0069   float middleStrip = nstrips / 2.;
0070   const LocalPoint& middleOfRoll = roll->centreOfStrip(middleStrip);
0071   const GlobalPoint& globMiddleRol = roll->toGlobal(middleOfRoll);
0072   double muRadius = sqrt(globMiddleRol.x() * globMiddleRol.x() + globMiddleRol.y() * globMiddleRol.y() +
0073                          globMiddleRol.z() * globMiddleRol.z());
0074   double timeCalibrationOffset_ = muRadius / cspeed;  //[ns]
0075
0076   const GEMStripTopology* top(dynamic_cast<const GEMStripTopology*>(&(roll->topology())));
0077   const float halfStripLength(0.5 * top->stripLength());
0078   const float distanceFromEdge(halfStripLength - simHitPos.y());
0079
0080   // signal propagation speed in material in [cm/ns]
0081   double signalPropagationSpeedTrue = signalPropagationSpeed_ * cspeed;
0082
0083   // average time for the signal to propagate from the SimHit to the top of a strip
0084   const float averagePropagationTime(distanceFromEdge / signalPropagationSpeedTrue);
0085   // random Gaussian time correction due to the finite timing resolution of the detector
0086   float randomResolutionTime = CLHEP::RandGaussQ::shoot(engine, 0., timeResolution_);
0087   const float simhitTime(tof + averageShapingTime_ + randomResolutionTime + averagePropagationTime + randomJitterTime);
0088   float referenceTime = 0.;
0089   referenceTime = timeCalibrationOffset_ + halfStripLength / signalPropagationSpeedTrue + averageShapingTime_;
0090   const float timeDifference(simhitTime - referenceTime);
0091   // assign the bunch crossing
0092   bx = static_cast<int>(std::round((timeDifference) / 25.));
0093
0094   // check time
0095   LogDebug("GEMDigiProducer") << "checktime "
0096                               << "bx = " << bx << "\tdeltaT = " << timeDifference << "\tsimT =  " << simhitTime
0097                               << "\trefT =  " << referenceTime << "\ttof = " << tof
0098                               << "\tavePropT =  " << averagePropagationTime
0099                               << "\taveRefPropT = " << halfStripLength / signalPropagationSpeedTrue << "\n";
0100   return bx;
0101 }
0102
0103 std::vector<std::pair<int, int> > GEMSignalModel::simulateClustering(const GEMStripTopology* top,
0104                                                                      const PSimHit* simHit,
0105                                                                      const int bx,
0106                                                                      CLHEP::HepRandomEngine* engine) {
0107   const LocalPoint& hit_entry(simHit->entryPoint());
0108   const LocalPoint& hit_exit(simHit->exitPoint());
0109
0110   LocalPoint start_point, end_point;
0111   if (hit_entry.x() < hit_exit.x()) {
0112     start_point = hit_entry;
0113     end_point = hit_exit;
0114   } else {
0115     start_point = hit_exit;
0116     end_point = hit_entry;
0117   }
0118
0119   // Add Gaussian noise to the points towards outside.
0120   float smeared_start_x = start_point.x() - std::abs(CLHEP::RandGaussQ::shoot(engine, 0, resolutionX_));
0121   float smeared_end_x = end_point.x() + std::abs(CLHEP::RandGaussQ::shoot(engine, 0, resolutionX_));
0122
0123   LocalPoint smeared_start_point(smeared_start_x, start_point.y(), start_point.z());
0124   LocalPoint smeared_end_point(smeared_end_x, end_point.y(), end_point.z());
0125
0126   int cluster_start = top->channel(smeared_start_point);
0127   int cluster_end = top->channel(smeared_end_point);
0128
0129   std::vector<std::pair<int, int> > cluster;
0130   for (int strip = cluster_start; strip <= cluster_end; strip++) {
0131     cluster.emplace_back(strip, bx);
0132   }
0133
0134   return cluster;
0135 }