GEMDigitizer/src/GEMBkgModel.cc

0001 #include "SimMuon/GEMDigitizer/interface/GEMBkgModel.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/Utilities/interface/Exception.h"
0006 #include "CLHEP/Random/RandFlat.h"
0007 #include "CLHEP/Random/RandPoissonQ.h"
0008 #include "CLHEP/Random/RandGaussQ.h"
0009 #include <cmath>
0010 #include <utility>
0011 #include <map>
0012
0013 GEMBkgModel::GEMBkgModel(const edm::ParameterSet& config)
0014     : GEMDigiModel(config),
0015       clusterSizeCut(0.53),
0016       averageEfficiency_(config.getParameter<double>("averageEfficiency")),
0017       minBunch_(config.getParameter<int>("minBunch")),
0018       maxBunch_(config.getParameter<int>("maxBunch")),
0019       simulateNoiseCLS_(config.getParameter<bool>("simulateNoiseCLS")),
0020       fixedRollRadius_(config.getParameter<bool>("fixedRollRadius")),
0021       simulateElectronBkg_(config.getParameter<bool>("simulateElectronBkg")),
0022       instLumi_(config.getParameter<double>("instLumi")),
0023       rateFact_(config.getParameter<double>("rateFact")),
0024       bxWidth_(config.getParameter<double>("bxWidth")),
0025       referenceInstLumi_(config.getParameter<double>("referenceInstLumi")),
0026       GE11ElecBkgParam0_(config.getParameter<double>("GE11ElecBkgParam0")),
0027       GE11ElecBkgParam1_(config.getParameter<double>("GE11ElecBkgParam1")),
0028       GE11ElecBkgParam2_(config.getParameter<double>("GE11ElecBkgParam2")),
0029       GE21ElecBkgParam0_(config.getParameter<double>("GE21ElecBkgParam0")),
0030       GE21ElecBkgParam1_(config.getParameter<double>("GE21ElecBkgParam1")),
0031       GE21ElecBkgParam2_(config.getParameter<double>("GE21ElecBkgParam2")),
0032       GE11ModNeuBkgParam0_(config.getParameter<double>("GE11ModNeuBkgParam0")),
0033       GE11ModNeuBkgParam1_(config.getParameter<double>("GE11ModNeuBkgParam1")),
0034       GE11ModNeuBkgParam2_(config.getParameter<double>("GE11ModNeuBkgParam2")),
0035       GE21ModNeuBkgParam0_(config.getParameter<double>("GE11ModNeuBkgParam0")),
0036       GE21ModNeuBkgParam1_(config.getParameter<double>("GE11ModNeuBkgParam1")),
0037       GE21ModNeuBkgParam2_(config.getParameter<double>("GE11ModNeuBkgParam2")) {}
0038
0039 GEMBkgModel::~GEMBkgModel() {}
0040
0041 void GEMBkgModel::simulate(const GEMEtaPartition* roll,
0042                            const edm::PSimHitContainer&,
0043                            CLHEP::HepRandomEngine* engine,
0044                            Strips& strips_,
0045                            DetectorHitMap& detectorHitMap_) {
0046   const GEMDetId& gemId(roll->id());
0047   const int nstrips(roll->nstrips());
0048   double trArea(0.0);
0049   double trStripArea(0.0);
0050   if (gemId.region() == 0) {
0051     throw cms::Exception("Geometry") << "GEMBkgModel::simulate() - this GEM id is from barrel, which cannot happen.";
0052   }
0053   const GEMStripTopology* top_(dynamic_cast<const GEMStripTopology*>(&(roll->topology())));
0054   const float striplength(top_->stripLength());
0055   trStripArea = (roll->pitch()) * striplength;
0056   trArea = trStripArea * nstrips;
0057   const int nBxing(maxBunch_ - minBunch_ + 1);
0058   const float rollRadius(
0059       fixedRollRadius_
0060           ? top_->radius()
0061           : top_->radius() + CLHEP::RandFlat::shoot(engine, -1. * top_->stripLength() / 2., top_->stripLength() / 2.));
0062
0063   //calculate noise from model
0064   double averageNeutralNoiseRatePerRoll = 0.;
0065   double averageNoiseElectronRatePerRoll = 0.;
0066   double averageNoiseRatePerRoll = 0.;
0067   if (gemId.station() == 1) {
0068     //simulate neutral background for GE1/1
0069     averageNeutralNoiseRatePerRoll =
0070         (GE11ModNeuBkgParam0_ + GE11ModNeuBkgParam1_ * rollRadius +
0071          GE11ModNeuBkgParam2_ * rollRadius * rollRadius);  //simulate electron background for GE1/1
0072     if (simulateElectronBkg_)
0073       averageNoiseElectronRatePerRoll =
0074           (GE11ElecBkgParam0_ + GE11ElecBkgParam1_ * rollRadius + GE11ElecBkgParam2_ * rollRadius * rollRadius);
0075
0076     // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
0077     averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
0078     averageNoiseRatePerRoll *= instLumi_ * rateFact_ / referenceInstLumi_;
0079   } else if (gemId.station() == 2) {
0080     //simulate neutral background for GE2/1
0081     averageNeutralNoiseRatePerRoll =
0082         (GE21ModNeuBkgParam0_ + GE21ModNeuBkgParam1_ * rollRadius + GE21ModNeuBkgParam2_ * rollRadius * rollRadius);
0083     //simulate electron background for GE2/1
0084     if (simulateElectronBkg_)
0085       averageNoiseElectronRatePerRoll =
0086           (GE21ElecBkgParam0_ + GE21ElecBkgParam1_ * rollRadius + GE21ElecBkgParam2_ * rollRadius * rollRadius);
0087
0088     // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
0089     averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
0090     averageNoiseRatePerRoll *= instLumi_ * rateFact_ / referenceInstLumi_;
0091   }
0092
0093   //simulate bkg contribution
0094   const double averageNoise(averageNoiseRatePerRoll * nBxing * trArea * bxWidth_);
0095   CLHEP::RandPoissonQ randPoissonQ(*engine, averageNoise);
0096   const int n_hits(randPoissonQ.fire());
0097   for (int i = 0; i < n_hits; ++i) {
0098     const int centralStrip(static_cast<int>(CLHEP::RandFlat::shoot(engine, 1, nstrips)));
0099     const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
0100     if (simulateNoiseCLS_) {
0101       std::vector<std::pair<int, int> > cluster_;
0102       cluster_.clear();
0103       cluster_.emplace_back(centralStrip, time_hit);
0104       int clusterSize((CLHEP::RandFlat::shoot(engine)) <= clusterSizeCut ? 1 : 2);
0105       if (clusterSize == 2) {
0106         if (CLHEP::RandFlat::shoot(engine) < 0.5) {
0107           if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
0108             cluster_.emplace_back(centralStrip - 1, time_hit);
0109         } else {
0110           if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 <= nstrips))
0111             cluster_.emplace_back(centralStrip + 1, time_hit);
0112         }
0113       }
0114       for (const auto& digi : cluster_) {
0115         strips_.emplace(digi);
0116       }
0117     }  //end simulateNoiseCLS_
0118     else {
0119       strips_.emplace(centralStrip, time_hit);
0120     }
0121   }
0122   return;
0123 }