ME0PreRecoGaussianModel.cc

Go to the documentation of this file.

#include "SimMuon/GEMDigitizer/interface/ME0PreRecoGaussianModel.h"
#include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h"
#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
#include "Geometry/GEMGeometry/interface/ME0Geometry.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandGaussQ.h"
#include "CLHEP/Random/RandPoissonQ.h"

#include <cmath>
#include <utility>
#include <map>

const int bxwidth = 25;  // [ns]

ME0PreRecoGaussianModel::ME0PreRecoGaussianModel(const edm::ParameterSet& config)
    : ME0DigiPreRecoModel(config),
      sigma_t(config.getParameter<double>("timeResolution")),
      sigma_u(config.getParameter<double>("phiResolution")),
      sigma_v(config.getParameter<double>("etaResolution")),
      error_u(config.getParameter<double>("phiError")),
      error_v(config.getParameter<double>("etaError")),
      gaussianSmearing_(config.getParameter<bool>("gaussianSmearing")),
      constPhiSmearing_(config.getParameter<bool>("constantPhiSpatialResolution")),
      corr(config.getParameter<bool>("useCorrelation")),
      etaproj(config.getParameter<bool>("useEtaProjectiveGEO")),
      digitizeOnlyMuons_(config.getParameter<bool>("digitizeOnlyMuons")),
      averageEfficiency_(config.getParameter<double>("averageEfficiency")),
      // simulateIntrinsicNoise_(config.getParameter<bool>("simulateIntrinsicNoise")),
      // averageNoiseRate_(config.getParameter<double>("averageNoiseRate")),
      simulateElectronBkg_(config.getParameter<bool>("simulateElectronBkg")),
      simulateNeutralBkg_(config.getParameter<bool>("simulateNeutralBkg")),
      minBunch_(config.getParameter<int>("minBunch")),
      maxBunch_(config.getParameter<int>("maxBunch")),
      instLumi_(config.getParameter<double>("instLumi")),
      rateFact_(config.getParameter<double>("rateFact")),
      referenceInstLumi_(config.getParameter<double>("referenceInstLumi")) {
  // polynomial parametrisation of neutral (n+g) and electron background
  // This is the background for an Instantaneous Luminosity of L = 5E34 cm^-2 s^-1
  neuBkg.push_back(0.00386257);
  neuBkg.push_back(6344.65);
  neuBkg.push_back(16627700);
  neuBkg.push_back(-102098);

  eleBkg.push_back(0.00171409);
  eleBkg.push_back(4900.56);
  eleBkg.push_back(710909);
  eleBkg.push_back(-4327.25);
}

ME0PreRecoGaussianModel::~ME0PreRecoGaussianModel() {}

void ME0PreRecoGaussianModel::simulateSignal(const ME0EtaPartition* roll,
                                             const edm::PSimHitContainer& simHits,
                                             CLHEP::HepRandomEngine* engine) {
  for (const auto& hit : simHits) {
    // Digitize only Muons?
    if (std::abs(hit.particleType()) != 13 && digitizeOnlyMuons_)
      continue;
    // Digitize only in [minBunch,maxBunch] window
    // window is: [(2n+1)*bxw/2, (2n+3)*bxw/2], n = [minBunch, maxBunch]
    if (hit.timeOfFlight() < (2 * minBunch_ + 1) * bxwidth * 1.0 / 2 ||
        hit.timeOfFlight() > (2 * maxBunch_ + 3) * bxwidth * 1.0 / 2)
      continue;
    // is GEM efficient?
    if (CLHEP::RandFlat::shoot(engine) > averageEfficiency_)
      continue;
    // create digi
    auto entry = hit.entryPoint();
    double x = 0.0, y = 0.0;

    double sigma_u_new = sigma_u;
    if (constPhiSmearing_)
      sigma_u_new = correctSigmaU(roll, entry.y());

    if (gaussianSmearing_) {  // Gaussian Smearing
      x = CLHEP::RandGaussQ::shoot(engine, entry.x(), sigma_u_new);
      y = CLHEP::RandGaussQ::shoot(engine, entry.y(), sigma_v);
    } else {  // Uniform Smearing ... use the sigmas as boundaries
      x = entry.x() + (CLHEP::RandFlat::shoot(engine) - 0.5) * sigma_u_new;
      y = entry.y() + (CLHEP::RandFlat::shoot(engine) - 0.5) * sigma_v;
    }
    double ex = sigma_u_new;
    double ey = sigma_v;
    double corr = 0.;
    double tof = CLHEP::RandGaussQ::shoot(engine, hit.timeOfFlight(), sigma_t);
    int pdgid = hit.particleType();
    if (ex == 0)
      ex = error_u;  //errors cannot be zero
    if (ey == 0)
      ey = error_v;

    int evtId = hit.eventId().event();
    int bx = hit.eventId().bunchCrossing();
    int procType = hit.processType();
    int res = 1;
    if (!(evtId == 0 && bx == 0 && procType == 0))
      res = 2;

    digi_.emplace(x, y, ex, ey, corr, tof, pdgid, res);

    edm::LogVerbatim("ME0PreRecoGaussianModel")
        << "[ME0PreRecoDigi :: simulateSignal] :: simhit in " << roll->id() << " at loc x = " << std::setw(8)
        << entry.x() << " [cm]"
        << " loc y = " << std::setw(8) << entry.y() << " [cm] time = " << std::setw(8) << hit.timeOfFlight()
        << " [ns] pdgid = " << std::showpos << std::setw(4) << pdgid;
    edm::LogVerbatim("ME0PreRecoGaussianModel")
        << "[ME0PreRecoDigi :: simulateSignal] :: digi   in " << roll->id() << " at loc x = " << std::setw(8) << x
        << " [cm] loc y = " << std::setw(8) << y << " [cm]"
        << " time = " << std::setw(8) << tof << " [ns]";
    edm::LogVerbatim("ME0PreRecoGaussianModel")
        << "[ME0PreRecoDigi :: simulateSignal] :: digi   in " << roll->id() << " with DX = " << std::setw(8)
        << (entry.x() - x) << " [cm]"
        << " DY = " << std::setw(8) << (entry.y() - y) << " [cm] DT = " << std::setw(8) << (hit.timeOfFlight() - tof)
        << " [ns]";
  }
}

void ME0PreRecoGaussianModel::simulateNoise(const ME0EtaPartition* roll, CLHEP::HepRandomEngine* engine) {
  double trArea(0.0);
  const ME0DetId me0Id(roll->id());

  // Extract detailed information from the Strip Topology:
  // base_bottom, base_top, height, strips, pads
  // note that (0,0) is in the middle of the roll ==> all param are at all half length
  if (me0Id.region() == 0) {
    throw cms::Exception("Geometry") << "Asking TrapezoidalStripTopology from a ME0 will fail";
  }  // not sure we really need this
  const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));

  auto& parameters(roll->specs()->parameters());
  double bottomLength(parameters[0]);
  bottomLength = 2 * bottomLength;  // bottom is largest length, so furtest away from beamline
  double topLength(parameters[1]);
  topLength = 2 * topLength;  // top is shortest length, so closest to beamline
  double height(parameters[2]);
  height = 2 * height;
  double myTanPhi = (topLength - bottomLength) / (height * 2);
  double rollRadius = top_->radius();
  trArea = height * (topLength + bottomLength) / 2.0;

  // Divide the detector area in different strips
  // take smearing in y-coord as height for each strip
  double initialHeight = sigma_v;
  if (sigma_v < 1.0)
    initialHeight = 1.0;
  double heightIt = initialHeight;
  int heightbins = height / heightIt;  // round down

  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
      << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: roll with id = " << roll->id();
  edm::LogVerbatim("ME0PreRecoGaussianModelNoise") << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId()
                                                   << "] :: extracting parameters from the TrapezoidalStripTopology";
  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
      << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: bottom = " << bottomLength
      << " [cm] top  = " << topLength << " [cm] height = " << height << " [cm]"
      << " area  = " << trArea << " [cm^2] Rmid = " << rollRadius
      << " [cm] => Rmin = " << rollRadius - height * 1.0 / 2.0 << " [cm] Rmax = " << rollRadius + height * 1.0 / 2.0
      << " [cm]";
  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
      << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: heightbins = " << heightbins;

  for (int hx = 0; hx < heightbins; ++hx) {
    double bottomIt = bottomLength + hx * 2 * tan(10. / 180 * 3.14) * heightIt;
    double topIt = bottomLength + (hx + 1) * 2 * tan(10. / 180 * 3.14) * heightIt;
    if (hx == heightbins - 1) {
      topIt = topLength;  // last bin ... make strip a bit larger to cover entire roll
      heightIt = height - hx * heightIt;
    }
    double areaIt = heightIt * (bottomIt + topIt) * 1.0 / 2;

    edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
        << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: height = " << std::setw(12) << heightIt
        << " [cm] bottom = " << std::setw(12) << bottomIt << " [cm]"
        << " top = " << std::setw(12) << topIt << " [cm] area = " << std::setw(12) << areaIt
        << " [cm^2] || sin(10) = " << sin(10. / 180 * 3.14);

    double myRandY = CLHEP::RandFlat::shoot(engine);
    double y0_rand = (hx + myRandY) * heightIt;  // Y coord, measured from the bottom of the roll
    if (hx == heightbins - 1)
      y0_rand = hx * initialHeight + myRandY * heightIt;
    double yy_rand =
        (y0_rand -
         height * 1.0 / 2);  // Y coord, measured from the middle of the roll, which is the Y coord in Local Coords
    double yy_glob = rollRadius + yy_rand;  // R coord in Global Coords
    // max length in x for given y coordinate (cfr trapezoidal eta partition)
    const float rSqrtR = yy_glob * sqrt(yy_glob);
    double xMax = topLength / 2.0 - (height / 2.0 - yy_rand) * myTanPhi;

    double sigma_u_new = sigma_u;
    if (constPhiSmearing_)
      sigma_u_new = correctSigmaU(roll, yy_rand);

    // 1) Intrinsic Noise ... Not implemented right now
    // ------------------------------------------------
    // if (simulateIntrinsicNoise_)
    // {
    // }

    // 2) Background Noise
    // ----------------------------

    // 2a) electron background
    // -----------------------
    if (simulateElectronBkg_) {
      // Extract / Calculate the Average Electron Rate
      // for the given global Y coord from Parametrization
      double averageElectronRatePerRoll =
          eleBkg[0] * rSqrtR * std::exp(eleBkg[1] / rSqrtR) + eleBkg[2] / rSqrtR + eleBkg[3] / (sqrt(yy_glob));
      // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
      averageElectronRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;

      // Rate [Hz/cm^2] * Nbx * 25*10^-9 [s] * Area [cm] = # hits in this roll in this bx
      const double averageElecRate(averageElectronRatePerRoll * (maxBunch_ - minBunch_ + 1) * (bxwidth * 1.0e-9) *
                                   areaIt);

      edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
          << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "]" /* "] :: BX = "<<std::showpos<<bx*/
          << " evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageElectronRatePerRoll
          << " [Hz/cm^2]"
          << " x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt << " [cm^2]) ==> " << std::setw(12)
          << averageElecRate << " [hits]";

      bool ele_eff = (CLHEP::RandFlat::shoot(engine) < averageElecRate) ? true : false;

      edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
          << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "] :: myRandY = " << std::setw(12) << myRandY
          << " => local y = " << std::setw(12) << yy_rand << " [cm]"
          << " => global y (global R) = " << std::setw(12) << yy_glob << " [cm] || Probability = " << std::setw(12)
          << averageElecRate << " => efficient? " << ele_eff << std::endl;

      if (ele_eff) {
        //calculate xx_rand at a given yy_rand
        double myRandX = CLHEP::RandFlat::shoot(engine);
        double xx_rand = 2 * xMax * (myRandX - 0.5);
        double ex = sigma_u_new;
        double ey = sigma_v;
        double corr = 0.;
        // extract random BX
        double myrandBX = CLHEP::RandFlat::shoot(engine);
        int bx = int((maxBunch_ - minBunch_ + 1) * myrandBX) + minBunch_;
        // extract random time in this BX
        double myrandT = CLHEP::RandFlat::shoot(engine);
        double minBXtime = (bx - 0.5) * bxwidth;  // double maxBXtime = (bx+0.5)*bxwidth;
        double time = myrandT * bxwidth + minBXtime;
        double myrandP = CLHEP::RandFlat::shoot(engine);
        int pdgid = 0;
        if (myrandP <= 0.5)
          pdgid = -11;  // electron
        else
          pdgid = 11;  // positron
        if (ex == 0)
          ex = error_u;  //errors cannot be zero
        if (ey == 0)
          ey = error_v;

        digi_.emplace(xx_rand, yy_rand, ex, ey, corr, time, pdgid, 0);

        edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
            << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "] =====> electron hit in " << roll->id()
            << " pdgid = " << pdgid << " bx = " << bx << " ==> digitized"
            << " at loc x = " << xx_rand << " loc y = " << yy_rand << " time = " << time << " [ns]";
      }
    }  // end if electron bkg

    // 2b) neutral (n+g) background
    // ----------------------------
    if (simulateNeutralBkg_) {
      // Extract / Calculate the Average Neutral Rate
      // for the given global Y coord from Parametrization
      double averageNeutralRatePerRoll =
          neuBkg[0] * yy_glob * std::exp(neuBkg[1] / rSqrtR) + neuBkg[2] / rSqrtR + neuBkg[3] / (sqrt(yy_glob));
      // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
      averageNeutralRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;

      // Rate [Hz/cm^2] * Nbx * 25*10^-9 [s] * Area [cm] = # hits in this roll
      const double averageNeutrRate(averageNeutralRatePerRoll * (maxBunch_ - minBunch_ + 1) * (bxwidth * 1.0e-9) *
                                    areaIt);

      edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
          << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "]" /* "] :: BX = "<<std::showpos<<bx*/
          << " evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageNeutralRatePerRoll
          << " [Hz/cm^2]"
          << " x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt << " [cm^2]) ==> " << std::setw(12)
          << averageNeutrRate << " [hits]";

      bool neu_eff = (CLHEP::RandFlat::shoot(engine) < averageNeutrRate) ? true : false;

      edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
          << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "] :: myRandY = " << std::setw(12) << myRandY
          << " => local y = " << std::setw(12) << yy_rand << " [cm]"
          << " => global y (global R) = " << std::setw(12) << yy_glob << " [cm] || Probability = " << std::setw(12)
          << averageNeutrRate << " => efficient? " << neu_eff << std::endl;

      if (neu_eff) {
        //calculate xx_rand at a given yy_rand
        double myRandX = CLHEP::RandFlat::shoot(engine);
        double xx_rand = 2 * xMax * (myRandX - 0.5);
        double ex = sigma_u_new;
        double ey = sigma_v;
        double corr = 0.;
        // extract random BX
        double myrandBX = CLHEP::RandFlat::shoot(engine);
        int bx = int((maxBunch_ - minBunch_ + 1) * myrandBX) + minBunch_;
        // extract random time in this BX
        double myrandT = CLHEP::RandFlat::shoot(engine);
        double minBXtime = (bx - 0.5) * bxwidth;
        double time = myrandT * bxwidth + minBXtime;
        int pdgid = 0;
        double myrandP = CLHEP::RandFlat::shoot(engine);
        if (myrandP <= 0.08)
          pdgid = 2112;  // neutrons: GEM sensitivity for neutrons: 0.08%
        else
          pdgid =
              22;  // photons:  GEM sensitivity for photons:  1.04% ==> neutron fraction = (0.08 / 1.04) = 0.077 = 0.08
        if (ex == 0)
          ex = error_u;  //errors cannot be zero
        if (ey == 0)
          ey = error_v;

        digi_.emplace(xx_rand, yy_rand, ex, ey, corr, time, pdgid, 0);

        edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
            << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "] ======> neutral hit in " << roll->id()
            << " pdgid = " << pdgid << " bx = " << bx << " ==> digitized"
            << " at loc x = " << xx_rand << " loc y = " << yy_rand << " time = " << time << " [ns]";
      }

    }  // end if neutral bkg
  }    // end loop over strips (= pseudo rolls)
}

double ME0PreRecoGaussianModel::correctSigmaU(const ME0EtaPartition* roll, double y) {
  const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
  auto& parameters(roll->specs()->parameters());
  double height(parameters[2]);        // height     = height from Center of Roll
  double rollRadius = top_->radius();  // rollRadius = Radius at Center of Roll
  double Rmax = rollRadius + height;   // MaxRadius  = Radius at top of Roll
  double Rx = rollRadius + y;          // y in [-height,+height]
  double sigma_u_new = Rx / Rmax * sigma_u;
  return sigma_u_new;
}