Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:46 

0001 #include "CLHEP/Units/defs.h"
0002 #include "CLHEP/Units/SystemOfUnits.h"
0003 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0004 #include "SimMuon/GEMDigitizer/interface/ME0SimpleModel.h"
0005 #include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h"
0006 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0007 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 #include "CLHEP/Random/RandFlat.h"
0011 #include "CLHEP/Random/RandPoissonQ.h"
0012 #include "CLHEP/Random/RandGaussQ.h"
0013 #include <cmath>
0014 #include <utility>
0015 #include <map>
0016 
0017 ME0SimpleModel::ME0SimpleModel(const edm::ParameterSet& config)
0018     : ME0DigiModel(config),
0019       averageEfficiency_(config.getParameter<double>("averageEfficiency")),
0020       averageShapingTime_(config.getParameter<double>("averageShapingTime")),
0021       timeResolution_(config.getParameter<double>("timeResolution")),
0022       timeJitter_(config.getParameter<double>("timeJitter")),
0023       averageNoiseRate_(config.getParameter<double>("averageNoiseRate")),
0024       signalPropagationSpeed_(config.getParameter<double>("signalPropagationSpeed")),
0025       bxwidth_(config.getParameter<int>("bxwidth")),
0026       minBunch_(config.getParameter<int>("minBunch")),
0027       maxBunch_(config.getParameter<int>("maxBunch")),
0028       digitizeOnlyMuons_(config.getParameter<bool>("digitizeOnlyMuons")),
0029       doBkgNoise_(config.getParameter<bool>("doBkgNoise")),
0030       doNoiseCLS_(config.getParameter<bool>("doNoiseCLS")),
0031       fixedRollRadius_(config.getParameter<bool>("fixedRollRadius")),
0032       simulateElectronBkg_(config.getParameter<bool>("simulateElectronBkg")),
0033       instLumi_(config.getParameter<double>("instLumi")),
0034       rateFact_(config.getParameter<double>("rateFact")),
0035       referenceInstLumi_(config.getParameter<double>("referenceInstLumi")),
0036       ME0ElecBkgParam0_(config.getParameter<double>("ME0ElecBkgParam0")),
0037       ME0ElecBkgParam1_(config.getParameter<double>("ME0ElecBkgParam1")),
0038       ME0ElecBkgParam2_(config.getParameter<double>("ME0ElecBkgParam2")),
0039       ME0ElecBkgParam3_(config.getParameter<double>("ME0ElecBkgParam3")),
0040       ME0NeuBkgParam0_(config.getParameter<double>("ME0NeuBkgParam0")),
0041       ME0NeuBkgParam1_(config.getParameter<double>("ME0NeuBkgParam1")),
0042       ME0NeuBkgParam2_(config.getParameter<double>("ME0NeuBkgParam2")),
0043       ME0NeuBkgParam3_(config.getParameter<double>("ME0NeuBkgParam3")) {}
0044 
0045 ME0SimpleModel::~ME0SimpleModel() {}
0046 
0047 void ME0SimpleModel::setup() { return; }
0048 
0049 void ME0SimpleModel::simulateSignal(const ME0EtaPartition* roll,
0050                                     const edm::PSimHitContainer& simHits,
0051                                     CLHEP::HepRandomEngine* engine) {
0052   stripDigiSimLinks_.clear();
0053   detectorHitMap_.clear();
0054   stripDigiSimLinks_ = StripDigiSimLinks(roll->id().rawId());
0055 
0056   theME0DigiSimLinks_.clear();
0057   theME0DigiSimLinks_ = ME0DigiSimLinks(roll->id().rawId());
0058   bool digiMuon = false;
0059   bool digiElec = false;
0060   for (const auto& hit : simHits) {
0061     if (std::abs(hit.particleType()) != 13 && digitizeOnlyMuons_)
0062       continue;
0063     double elecEff = 0.;
0064     double partMom = hit.pabs();
0065     double checkMuonEff = CLHEP::RandFlat::shoot(engine, 0., 1.);
0066     double checkElecEff = CLHEP::RandFlat::shoot(engine, 0., 1.);
0067     if (std::abs(hit.particleType()) == 13 && checkMuonEff < averageEfficiency_)
0068       digiMuon = true;
0069     if (std::abs(hit.particleType()) != 13)  //consider all non muon particles with me0 efficiency to electrons
0070     {
0071       if (partMom <= 1.95e-03)
0072         elecEff = 1.7e-05 * std::exp(2.1 * partMom * 1000.);
0073       if (partMom > 1.95e-03 && partMom < 10.e-03)
0074         elecEff =
0075             1.34 * log(7.96e-01 * partMom * 1000. - 5.75e-01) / (1.34 + log(7.96e-01 * partMom * 1000. - 5.75e-01));
0076       if (partMom > 10.e-03)
0077         elecEff = 1.;
0078       if (checkElecEff < elecEff)
0079         digiElec = true;
0080     }
0081     if (!(digiMuon || digiElec))
0082       continue;
0083     const int bx(getSimHitBx(&hit, engine));
0084     const std::vector<std::pair<int, int> >& cluster(simulateClustering(roll, &hit, bx, engine));
0085     for (const auto& digi : cluster) {
0086       detectorHitMap_.emplace(digi, &hit);
0087       strips_.emplace(digi);
0088     }
0089   }
0090 }
0091 
0092 int ME0SimpleModel::getSimHitBx(const PSimHit* simhit, CLHEP::HepRandomEngine* engine) {
0093   int bx = -999;
0094   const LocalPoint& simHitPos(simhit->localPosition());
0095   const float tof(simhit->timeOfFlight());
0096   // random Gaussian time correction due to electronics jitter
0097   float randomJitterTime = CLHEP::RandGaussQ::shoot(engine, 0., timeJitter_);
0098   const ME0DetId& id(simhit->detUnitId());
0099   const ME0EtaPartition* roll(geometry_->etaPartition(id));
0100   if (!roll) {
0101     throw cms::Exception("Geometry") << "ME0SimpleModel::getSimHitBx() - ME0 simhit id does not match any ME0 roll id: "
0102                                      << id << "\n";
0103     return 999;
0104   }
0105   if (roll->id().region() == 0) {
0106     throw cms::Exception("Geometry")
0107         << "ME0SimpleModel::getSimHitBx() - this ME0 id is from barrel, which cannot happen: " << roll->id() << "\n";
0108   }
0109   const int nstrips = roll->nstrips();
0110   float middleStrip = nstrips / 2.;
0111   const LocalPoint& middleOfRoll = roll->centreOfStrip(middleStrip);
0112   const GlobalPoint& globMiddleRol = roll->toGlobal(middleOfRoll);
0113   double muRadius = sqrt(globMiddleRol.x() * globMiddleRol.x() + globMiddleRol.y() * globMiddleRol.y() +
0114                          globMiddleRol.z() * globMiddleRol.z());
0115   double timeCalibrationOffset_ = (muRadius * CLHEP::ns * CLHEP::cm) / (CLHEP::c_light);  //[cm/ns]
0116   const TrapezoidalStripTopology* top(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
0117   const float halfStripLength(0.5 * top->stripLength());
0118   const float distanceFromEdge(halfStripLength - simHitPos.y());
0119 
0120   // signal propagation speed in material in [cm/ns]
0121   double signalPropagationSpeedTrue = signalPropagationSpeed_ * CLHEP::c_light / (CLHEP::ns * CLHEP::cm);
0122   // average time for the signal to propagate from the SimHit to the top of a strip
0123   const float averagePropagationTime(distanceFromEdge / signalPropagationSpeedTrue);
0124   // random Gaussian time correction due to the finite timing resolution of the detector
0125   float randomResolutionTime = CLHEP::RandGaussQ::shoot(engine, 0., timeResolution_);
0126   const float simhitTime(tof + averageShapingTime_ + randomResolutionTime + averagePropagationTime + randomJitterTime);
0127   float referenceTime = 0.;
0128   referenceTime = timeCalibrationOffset_ + halfStripLength / signalPropagationSpeedTrue + averageShapingTime_;
0129   const float timeDifference(simhitTime - referenceTime);
0130   // assign the bunch crossing
0131   bx = static_cast<int>(std::round((timeDifference) / bxwidth_));
0132 
0133   // check time
0134   const bool debug(false);
0135   if (debug) {
0136     LogDebug("ME0SimpleModel") << "checktime "
0137                                << "bx = " << bx << "\tdeltaT = " << timeDifference << "\tsimT =  " << simhitTime
0138                                << "\trefT =  " << referenceTime << "\ttof = " << tof
0139                                << "\tavePropT =  " << averagePropagationTime
0140                                << "\taveRefPropT = " << halfStripLength / signalPropagationSpeedTrue << std::endl;
0141   }
0142   return bx;
0143 }
0144 
0145 void ME0SimpleModel::simulateNoise(const ME0EtaPartition* roll, CLHEP::HepRandomEngine* engine) {
0146   if (!doBkgNoise_)
0147     return;
0148   const ME0DetId me0Id(roll->id());
0149   const int nstrips(roll->nstrips());
0150   double trArea(0.0);
0151   double trStripArea(0.0);
0152 
0153   if (me0Id.region() == 0) {
0154     throw cms::Exception("Geometry")
0155         << "ME0Synchronizer::simulateNoise() - this ME0 id is from barrel, which cannot happen.";
0156   }
0157   const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
0158   const float striplength(top_->stripLength());
0159   trStripArea = (roll->pitch()) * striplength;
0160   trArea = trStripArea * nstrips;
0161   const int nBxing(maxBunch_ - minBunch_ + 1);
0162   const float rollRadius(
0163       fixedRollRadius_
0164           ? top_->radius()
0165           : top_->radius() + CLHEP::RandFlat::shoot(engine, -1. * top_->stripLength() / 2., top_->stripLength() / 2.));
0166 
0167   const float rSqrtR = rollRadius * sqrt(rollRadius);
0168 
0169   //calculate noise from model
0170   double averageNeutralNoiseRatePerRoll = 0.;
0171   double averageNoiseElectronRatePerRoll = 0.;
0172   double averageNoiseRatePerRoll = 0.;
0173   if (me0Id.station() != 1) {
0174     throw cms::Exception("Geometry")
0175         << "ME0SimpleModel::simulateNoise() - this ME0 id is from station 1, which cannot happen: " << roll->id()
0176         << "\n";
0177   } else {
0178     averageNeutralNoiseRatePerRoll = ME0NeuBkgParam0_ * rollRadius * std::exp(ME0NeuBkgParam1_ / rSqrtR) +
0179                                      ME0NeuBkgParam2_ / rSqrtR + ME0NeuBkgParam3_ / (sqrt(rollRadius));
0180     //simulate electron background for ME0
0181     if (simulateElectronBkg_)
0182       averageNoiseElectronRatePerRoll = ME0ElecBkgParam0_ * rSqrtR * std::exp(ME0ElecBkgParam1_ / rSqrtR) +
0183                                         ME0ElecBkgParam2_ / rSqrtR + ME0ElecBkgParam3_ / (sqrt(rollRadius));
0184     averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
0185     averageNoiseRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;
0186   }
0187 
0188   //simulate intrinsic noise
0189   if (simulateIntrinsicNoise_) {
0190     const double aveIntrinsicNoisePerStrip(averageNoiseRate_ * nBxing * bxwidth_ * trStripArea * 1.0e-9);
0191     for (int j = 0; j < nstrips; ++j) {
0192       int randPoissonQ = CLHEP::RandPoissonQ::shoot(engine, aveIntrinsicNoisePerStrip);
0193 
0194       for (int k = 0; k < randPoissonQ; k++) {
0195         const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
0196         strips_.emplace(k + 1, time_hit);
0197       }
0198     }
0199   }  //end simulate intrinsic noise
0200 
0201   //simulate bkg contribution
0202   const double averageNoise(averageNoiseRatePerRoll * nBxing * bxwidth_ * trArea * 1.0e-9);
0203   int randPoissonQ = CLHEP::RandPoissonQ::shoot(engine, averageNoise);
0204 
0205   for (int i = 0; i < randPoissonQ; ++i) {
0206     const int centralStrip(static_cast<int>(CLHEP::RandFlat::shoot(engine, 1, nstrips)));
0207     const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
0208     if (doNoiseCLS_) {
0209       std::vector<std::pair<int, int> > cluster;
0210       cluster.emplace_back(centralStrip, time_hit);
0211       int clusterSize((CLHEP::RandFlat::shoot(engine)) <= 0.53 ? 1 : 2);
0212       if (clusterSize == 2) {
0213         if (CLHEP::RandFlat::shoot(engine) < 0.5) {
0214           if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
0215             cluster.emplace_back(centralStrip - 1, time_hit);
0216         } else {
0217           if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 <= nstrips))
0218             cluster.emplace_back(centralStrip + 1, time_hit);
0219         }
0220       }
0221       for (const auto& digi : cluster) {
0222         strips_.emplace(digi);
0223       }
0224     }  //end doNoiseCLS_
0225     else {
0226       strips_.emplace(centralStrip, time_hit);
0227     }
0228   }
0229   return;
0230 }
0231 
0232 std::vector<std::pair<int, int> > ME0SimpleModel::simulateClustering(const ME0EtaPartition* roll,
0233                                                                      const PSimHit* simHit,
0234                                                                      const int bx,
0235                                                                      CLHEP::HepRandomEngine* engine) {
0236   const StripTopology& topology = roll->specificTopology();
0237   const LocalPoint& hit_position(simHit->localPosition());
0238   const int nstrips(roll->nstrips());
0239   int centralStrip = 0;
0240   if (!(topology.channel(hit_position) + 1 > nstrips))
0241     centralStrip = topology.channel(hit_position) + 1;
0242   else
0243     centralStrip = topology.channel(hit_position);
0244   const GlobalPoint& pointSimHit = roll->toGlobal(hit_position);
0245   const GlobalPoint& pointDigiHit = roll->toGlobal(roll->centreOfStrip(centralStrip));
0246   double deltaX = pointSimHit.x() - pointDigiHit.x();
0247 
0248   // Add central digi to cluster vector
0249   std::vector<std::pair<int, int> > cluster;
0250   cluster.clear();
0251   cluster.emplace_back(centralStrip, bx);
0252 
0253   //simulate cross talk
0254   int clusterSize((CLHEP::RandFlat::shoot(engine)) <= 0.53 ? 1 : 2);
0255   if (clusterSize == 2) {
0256     if (deltaX <= 0) {
0257       if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
0258         cluster.emplace_back(centralStrip - 1, bx);
0259     } else {
0260       if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 <= nstrips))
0261         cluster.emplace_back(centralStrip + 1, bx);
0262     }
0263   }
0264   return cluster;
0265 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!