ME0SimpleModel.cc

Go to the documentation of this file.

#include "CLHEP/Units/defs.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Units/GlobalPhysicalConstants.h"
#include "SimMuon/GEMDigitizer/interface/ME0SimpleModel.h"
#include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h"
#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
#include "Geometry/GEMGeometry/interface/ME0Geometry.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandPoissonQ.h"
#include "CLHEP/Random/RandGaussQ.h"
#include <cmath>
#include <utility>
#include <map>

ME0SimpleModel::ME0SimpleModel(const edm::ParameterSet& config)
    : ME0DigiModel(config),
      averageEfficiency_(config.getParameter<double>("averageEfficiency")),
      averageShapingTime_(config.getParameter<double>("averageShapingTime")),
      timeResolution_(config.getParameter<double>("timeResolution")),
      timeJitter_(config.getParameter<double>("timeJitter")),
      averageNoiseRate_(config.getParameter<double>("averageNoiseRate")),
      signalPropagationSpeed_(config.getParameter<double>("signalPropagationSpeed")),
      bxwidth_(config.getParameter<int>("bxwidth")),
      minBunch_(config.getParameter<int>("minBunch")),
      maxBunch_(config.getParameter<int>("maxBunch")),
      digitizeOnlyMuons_(config.getParameter<bool>("digitizeOnlyMuons")),
      doBkgNoise_(config.getParameter<bool>("doBkgNoise")),
      doNoiseCLS_(config.getParameter<bool>("doNoiseCLS")),
      fixedRollRadius_(config.getParameter<bool>("fixedRollRadius")),
      simulateElectronBkg_(config.getParameter<bool>("simulateElectronBkg")),
      instLumi_(config.getParameter<double>("instLumi")),
      rateFact_(config.getParameter<double>("rateFact")),
      referenceInstLumi_(config.getParameter<double>("referenceInstLumi")),
      ME0ElecBkgParam0_(config.getParameter<double>("ME0ElecBkgParam0")),
      ME0ElecBkgParam1_(config.getParameter<double>("ME0ElecBkgParam1")),
      ME0ElecBkgParam2_(config.getParameter<double>("ME0ElecBkgParam2")),
      ME0ElecBkgParam3_(config.getParameter<double>("ME0ElecBkgParam3")),
      ME0NeuBkgParam0_(config.getParameter<double>("ME0NeuBkgParam0")),
      ME0NeuBkgParam1_(config.getParameter<double>("ME0NeuBkgParam1")),
      ME0NeuBkgParam2_(config.getParameter<double>("ME0NeuBkgParam2")),
      ME0NeuBkgParam3_(config.getParameter<double>("ME0NeuBkgParam3")) {}

ME0SimpleModel::~ME0SimpleModel() {}

void ME0SimpleModel::setup() { return; }

void ME0SimpleModel::simulateSignal(const ME0EtaPartition* roll,
                                    const edm::PSimHitContainer& simHits,
                                    CLHEP::HepRandomEngine* engine) {
  stripDigiSimLinks_.clear();
  detectorHitMap_.clear();
  stripDigiSimLinks_ = StripDigiSimLinks(roll->id().rawId());

  theME0DigiSimLinks_.clear();
  theME0DigiSimLinks_ = ME0DigiSimLinks(roll->id().rawId());
  bool digiMuon = false;
  bool digiElec = false;
  for (const auto& hit : simHits) {
    if (std::abs(hit.particleType()) != 13 && digitizeOnlyMuons_)
      continue;
    double elecEff = 0.;
    double partMom = hit.pabs();
    double checkMuonEff = CLHEP::RandFlat::shoot(engine, 0., 1.);
    double checkElecEff = CLHEP::RandFlat::shoot(engine, 0., 1.);
    if (std::abs(hit.particleType()) == 13 && checkMuonEff < averageEfficiency_)
      digiMuon = true;
    if (std::abs(hit.particleType()) != 13)  //consider all non muon particles with me0 efficiency to electrons
    {
      if (partMom <= 1.95e-03)
        elecEff = 1.7e-05 * std::exp(2.1 * partMom * 1000.);
      if (partMom > 1.95e-03 && partMom < 10.e-03)
        elecEff =
            1.34 * log(7.96e-01 * partMom * 1000. - 5.75e-01) / (1.34 + log(7.96e-01 * partMom * 1000. - 5.75e-01));
      if (partMom > 10.e-03)
        elecEff = 1.;
      if (checkElecEff < elecEff)
        digiElec = true;
    }
    if (!(digiMuon || digiElec))
      continue;
    const int bx(getSimHitBx(&hit, engine));
    const std::vector<std::pair<int, int> >& cluster(simulateClustering(roll, &hit, bx, engine));
    for (const auto& digi : cluster) {
      detectorHitMap_.emplace(digi, &hit);
      strips_.emplace(digi);
    }
  }
}

int ME0SimpleModel::getSimHitBx(const PSimHit* simhit, CLHEP::HepRandomEngine* engine) {
  int bx = -999;
  const LocalPoint& simHitPos(simhit->localPosition());
  const float tof(simhit->timeOfFlight());
  // random Gaussian time correction due to electronics jitter
  float randomJitterTime = CLHEP::RandGaussQ::shoot(engine, 0., timeJitter_);
  const ME0DetId& id(simhit->detUnitId());
  const ME0EtaPartition* roll(geometry_->etaPartition(id));
  if (!roll) {
    throw cms::Exception("Geometry") << "ME0SimpleModel::getSimHitBx() - ME0 simhit id does not match any ME0 roll id: "
                                     << id << "\n";
    return 999;
  }
  if (roll->id().region() == 0) {
    throw cms::Exception("Geometry")
        << "ME0SimpleModel::getSimHitBx() - this ME0 id is from barrel, which cannot happen: " << roll->id() << "\n";
  }
  const int nstrips = roll->nstrips();
  float middleStrip = nstrips / 2.;
  const LocalPoint& middleOfRoll = roll->centreOfStrip(middleStrip);
  const GlobalPoint& globMiddleRol = roll->toGlobal(middleOfRoll);
  double muRadius = sqrt(globMiddleRol.x() * globMiddleRol.x() + globMiddleRol.y() * globMiddleRol.y() +
                         globMiddleRol.z() * globMiddleRol.z());
  double timeCalibrationOffset_ = (muRadius * CLHEP::ns * CLHEP::cm) / (CLHEP::c_light);  //[cm/ns]
  const TrapezoidalStripTopology* top(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
  const float halfStripLength(0.5 * top->stripLength());
  const float distanceFromEdge(halfStripLength - simHitPos.y());

  // signal propagation speed in material in [cm/ns]
  double signalPropagationSpeedTrue = signalPropagationSpeed_ * CLHEP::c_light / (CLHEP::ns * CLHEP::cm);
  // average time for the signal to propagate from the SimHit to the top of a strip
  const float averagePropagationTime(distanceFromEdge / signalPropagationSpeedTrue);
  // random Gaussian time correction due to the finite timing resolution of the detector
  float randomResolutionTime = CLHEP::RandGaussQ::shoot(engine, 0., timeResolution_);
  const float simhitTime(tof + averageShapingTime_ + randomResolutionTime + averagePropagationTime + randomJitterTime);
  float referenceTime = 0.;
  referenceTime = timeCalibrationOffset_ + halfStripLength / signalPropagationSpeedTrue + averageShapingTime_;
  const float timeDifference(simhitTime - referenceTime);
  // assign the bunch crossing
  bx = static_cast<int>(std::round((timeDifference) / bxwidth_));

  // check time
  const bool debug(false);
  if (debug) {
    LogDebug("ME0SimpleModel") << "checktime "
                               << "bx = " << bx << "\tdeltaT = " << timeDifference << "\tsimT =  " << simhitTime
                               << "\trefT =  " << referenceTime << "\ttof = " << tof
                               << "\tavePropT =  " << averagePropagationTime
                               << "\taveRefPropT = " << halfStripLength / signalPropagationSpeedTrue << std::endl;
  }
  return bx;
}

void ME0SimpleModel::simulateNoise(const ME0EtaPartition* roll, CLHEP::HepRandomEngine* engine) {
  if (!doBkgNoise_)
    return;
  const ME0DetId me0Id(roll->id());
  const int nstrips(roll->nstrips());
  double trArea(0.0);
  double trStripArea(0.0);

  if (me0Id.region() == 0) {
    throw cms::Exception("Geometry")
        << "ME0Synchronizer::simulateNoise() - this ME0 id is from barrel, which cannot happen.";
  }
  const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
  const float striplength(top_->stripLength());
  trStripArea = (roll->pitch()) * striplength;
  trArea = trStripArea * nstrips;
  const int nBxing(maxBunch_ - minBunch_ + 1);
  const float rollRadius(
      fixedRollRadius_
          ? top_->radius()
          : top_->radius() + CLHEP::RandFlat::shoot(engine, -1. * top_->stripLength() / 2., top_->stripLength() / 2.));

  const float rSqrtR = rollRadius * sqrt(rollRadius);

  //calculate noise from model
  double averageNeutralNoiseRatePerRoll = 0.;
  double averageNoiseElectronRatePerRoll = 0.;
  double averageNoiseRatePerRoll = 0.;
  if (me0Id.station() != 1) {
    throw cms::Exception("Geometry")
        << "ME0SimpleModel::simulateNoise() - this ME0 id is from station 1, which cannot happen: " << roll->id()
        << "\n";
  } else {
    averageNeutralNoiseRatePerRoll = ME0NeuBkgParam0_ * rollRadius * std::exp(ME0NeuBkgParam1_ / rSqrtR) +
                                     ME0NeuBkgParam2_ / rSqrtR + ME0NeuBkgParam3_ / (sqrt(rollRadius));
    //simulate electron background for ME0
    if (simulateElectronBkg_)
      averageNoiseElectronRatePerRoll = ME0ElecBkgParam0_ * rSqrtR * std::exp(ME0ElecBkgParam1_ / rSqrtR) +
                                        ME0ElecBkgParam2_ / rSqrtR + ME0ElecBkgParam3_ / (sqrt(rollRadius));
    averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
    averageNoiseRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;
  }

  //simulate intrinsic noise
  if (simulateIntrinsicNoise_) {
    const double aveIntrinsicNoisePerStrip(averageNoiseRate_ * nBxing * bxwidth_ * trStripArea * 1.0e-9);
    for (int j = 0; j < nstrips; ++j) {
      int randPoissonQ = CLHEP::RandPoissonQ::shoot(engine, aveIntrinsicNoisePerStrip);

      for (int k = 0; k < randPoissonQ; k++) {
        const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
        strips_.emplace(k + 1, time_hit);
      }
    }
  }  //end simulate intrinsic noise

  //simulate bkg contribution
  const double averageNoise(averageNoiseRatePerRoll * nBxing * bxwidth_ * trArea * 1.0e-9);
  int randPoissonQ = CLHEP::RandPoissonQ::shoot(engine, averageNoise);

  for (int i = 0; i < randPoissonQ; ++i) {
    const int centralStrip(static_cast<int>(CLHEP::RandFlat::shoot(engine, 1, nstrips)));
    const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
    if (doNoiseCLS_) {
      std::vector<std::pair<int, int> > cluster;
      cluster.emplace_back(centralStrip, time_hit);
      int clusterSize((CLHEP::RandFlat::shoot(engine)) <= 0.53 ? 1 : 2);
      if (clusterSize == 2) {
        if (CLHEP::RandFlat::shoot(engine) < 0.5) {
          if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
            cluster.emplace_back(centralStrip - 1, time_hit);
        } else {
          if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 <= nstrips))
            cluster.emplace_back(centralStrip + 1, time_hit);
        }
      }
      for (const auto& digi : cluster) {
        strips_.emplace(digi);
      }
    }  //end doNoiseCLS_
    else {
      strips_.emplace(centralStrip, time_hit);
    }
  }
  return;
}

std::vector<std::pair<int, int> > ME0SimpleModel::simulateClustering(const ME0EtaPartition* roll,
                                                                     const PSimHit* simHit,
                                                                     const int bx,
                                                                     CLHEP::HepRandomEngine* engine) {
  const StripTopology& topology = roll->specificTopology();
  const LocalPoint& hit_position(simHit->localPosition());
  const int nstrips(roll->nstrips());
  int centralStrip = 0;
  if (!(topology.channel(hit_position) + 1 > nstrips))
    centralStrip = topology.channel(hit_position) + 1;
  else
    centralStrip = topology.channel(hit_position);
  const GlobalPoint& pointSimHit = roll->toGlobal(hit_position);
  const GlobalPoint& pointDigiHit = roll->toGlobal(roll->centreOfStrip(centralStrip));
  double deltaX = pointSimHit.x() - pointDigiHit.x();

  // Add central digi to cluster vector
  std::vector<std::pair<int, int> > cluster;
  cluster.clear();
  cluster.emplace_back(centralStrip, bx);

  //simulate cross talk
  int clusterSize((CLHEP::RandFlat::shoot(engine)) <= 0.53 ? 1 : 2);
  if (clusterSize == 2) {
    if (deltaX <= 0) {
      if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
        cluster.emplace_back(centralStrip - 1, bx);
    } else {
      if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 <= nstrips))
        cluster.emplace_back(centralStrip + 1, bx);
    }
  }
  return cluster;
}