Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:25:24 

0001 #include "SimMuon/GEMDigitizer/interface/ME0PreRecoGaussianModel.h"
0002 #include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h"
0003 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0004 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0005 
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 
0009 #include "CLHEP/Random/RandFlat.h"
0010 #include "CLHEP/Random/RandGaussQ.h"
0011 #include "CLHEP/Random/RandPoissonQ.h"
0012 
0013 #include <cmath>
0014 #include <utility>
0015 #include <map>
0016 
0017 const int bxwidth = 25;  // [ns]
0018 
0019 ME0PreRecoGaussianModel::ME0PreRecoGaussianModel(const edm::ParameterSet& config)
0020     : ME0DigiPreRecoModel(config),
0021       sigma_t(config.getParameter<double>("timeResolution")),
0022       sigma_u(config.getParameter<double>("phiResolution")),
0023       sigma_v(config.getParameter<double>("etaResolution")),
0024       error_u(config.getParameter<double>("phiError")),
0025       error_v(config.getParameter<double>("etaError")),
0026       gaussianSmearing_(config.getParameter<bool>("gaussianSmearing")),
0027       constPhiSmearing_(config.getParameter<bool>("constantPhiSpatialResolution")),
0028       corr(config.getParameter<bool>("useCorrelation")),
0029       etaproj(config.getParameter<bool>("useEtaProjectiveGEO")),
0030       digitizeOnlyMuons_(config.getParameter<bool>("digitizeOnlyMuons")),
0031       averageEfficiency_(config.getParameter<double>("averageEfficiency")),
0032       // simulateIntrinsicNoise_(config.getParameter<bool>("simulateIntrinsicNoise")),
0033       // averageNoiseRate_(config.getParameter<double>("averageNoiseRate")),
0034       simulateElectronBkg_(config.getParameter<bool>("simulateElectronBkg")),
0035       simulateNeutralBkg_(config.getParameter<bool>("simulateNeutralBkg")),
0036       minBunch_(config.getParameter<int>("minBunch")),
0037       maxBunch_(config.getParameter<int>("maxBunch")),
0038       instLumi_(config.getParameter<double>("instLumi")),
0039       rateFact_(config.getParameter<double>("rateFact")),
0040       referenceInstLumi_(config.getParameter<double>("referenceInstLumi")) {
0041   // polynomial parametrisation of neutral (n+g) and electron background
0042   // This is the background for an Instantaneous Luminosity of L = 5E34 cm^-2 s^-1
0043   neuBkg.push_back(0.00386257);
0044   neuBkg.push_back(6344.65);
0045   neuBkg.push_back(16627700);
0046   neuBkg.push_back(-102098);
0047 
0048   eleBkg.push_back(0.00171409);
0049   eleBkg.push_back(4900.56);
0050   eleBkg.push_back(710909);
0051   eleBkg.push_back(-4327.25);
0052 }
0053 
0054 ME0PreRecoGaussianModel::~ME0PreRecoGaussianModel() {}
0055 
0056 void ME0PreRecoGaussianModel::simulateSignal(const ME0EtaPartition* roll,
0057                                              const edm::PSimHitContainer& simHits,
0058                                              CLHEP::HepRandomEngine* engine) {
0059   for (const auto& hit : simHits) {
0060     // Digitize only Muons?
0061     if (std::abs(hit.particleType()) != 13 && digitizeOnlyMuons_)
0062       continue;
0063     // Digitize only in [minBunch,maxBunch] window
0064     // window is: [(2n+1)*bxw/2, (2n+3)*bxw/2], n = [minBunch, maxBunch]
0065     if (hit.timeOfFlight() < (2 * minBunch_ + 1) * bxwidth * 1.0 / 2 ||
0066         hit.timeOfFlight() > (2 * maxBunch_ + 3) * bxwidth * 1.0 / 2)
0067       continue;
0068     // is GEM efficient?
0069     if (CLHEP::RandFlat::shoot(engine) > averageEfficiency_)
0070       continue;
0071     // create digi
0072     auto entry = hit.entryPoint();
0073     double x = 0.0, y = 0.0;
0074 
0075     double sigma_u_new = sigma_u;
0076     if (constPhiSmearing_)
0077       sigma_u_new = correctSigmaU(roll, entry.y());
0078 
0079     if (gaussianSmearing_) {  // Gaussian Smearing
0080       x = CLHEP::RandGaussQ::shoot(engine, entry.x(), sigma_u_new);
0081       y = CLHEP::RandGaussQ::shoot(engine, entry.y(), sigma_v);
0082     } else {  // Uniform Smearing ... use the sigmas as boundaries
0083       x = entry.x() + (CLHEP::RandFlat::shoot(engine) - 0.5) * sigma_u_new;
0084       y = entry.y() + (CLHEP::RandFlat::shoot(engine) - 0.5) * sigma_v;
0085     }
0086     double ex = sigma_u_new;
0087     double ey = sigma_v;
0088     double corr = 0.;
0089     double tof = CLHEP::RandGaussQ::shoot(engine, hit.timeOfFlight(), sigma_t);
0090     int pdgid = hit.particleType();
0091     if (ex == 0)
0092       ex = error_u;  //errors cannot be zero
0093     if (ey == 0)
0094       ey = error_v;
0095 
0096     int evtId = hit.eventId().event();
0097     int bx = hit.eventId().bunchCrossing();
0098     int procType = hit.processType();
0099     int res = 1;
0100     if (!(evtId == 0 && bx == 0 && procType == 0))
0101       res = 2;
0102 
0103     digi_.emplace(x, y, ex, ey, corr, tof, pdgid, res);
0104 
0105     edm::LogVerbatim("ME0PreRecoGaussianModel")
0106         << "[ME0PreRecoDigi :: simulateSignal] :: simhit in " << roll->id() << " at loc x = " << std::setw(8)
0107         << entry.x() << " [cm]"
0108         << " loc y = " << std::setw(8) << entry.y() << " [cm] time = " << std::setw(8) << hit.timeOfFlight()
0109         << " [ns] pdgid = " << std::showpos << std::setw(4) << pdgid;
0110     edm::LogVerbatim("ME0PreRecoGaussianModel")
0111         << "[ME0PreRecoDigi :: simulateSignal] :: digi   in " << roll->id() << " at loc x = " << std::setw(8) << x
0112         << " [cm] loc y = " << std::setw(8) << y << " [cm]"
0113         << " time = " << std::setw(8) << tof << " [ns]";
0114     edm::LogVerbatim("ME0PreRecoGaussianModel")
0115         << "[ME0PreRecoDigi :: simulateSignal] :: digi   in " << roll->id() << " with DX = " << std::setw(8)
0116         << (entry.x() - x) << " [cm]"
0117         << " DY = " << std::setw(8) << (entry.y() - y) << " [cm] DT = " << std::setw(8) << (hit.timeOfFlight() - tof)
0118         << " [ns]";
0119   }
0120 }
0121 
0122 void ME0PreRecoGaussianModel::simulateNoise(const ME0EtaPartition* roll, CLHEP::HepRandomEngine* engine) {
0123   double trArea(0.0);
0124   const ME0DetId me0Id(roll->id());
0125 
0126   // Extract detailed information from the Strip Topology:
0127   // base_bottom, base_top, height, strips, pads
0128   // note that (0,0) is in the middle of the roll ==> all param are at all half length
0129   if (me0Id.region() == 0) {
0130     throw cms::Exception("Geometry") << "Asking TrapezoidalStripTopology from a ME0 will fail";
0131   }  // not sure we really need this
0132   const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
0133 
0134   auto& parameters(roll->specs()->parameters());
0135   double bottomLength(parameters[0]);
0136   bottomLength = 2 * bottomLength;  // bottom is largest length, so furtest away from beamline
0137   double topLength(parameters[1]);
0138   topLength = 2 * topLength;  // top is shortest length, so closest to beamline
0139   double height(parameters[2]);
0140   height = 2 * height;
0141   double myTanPhi = (topLength - bottomLength) / (height * 2);
0142   double rollRadius = top_->radius();
0143   trArea = height * (topLength + bottomLength) / 2.0;
0144 
0145   // Divide the detector area in different strips
0146   // take smearing in y-coord as height for each strip
0147   double initialHeight = sigma_v;
0148   if (sigma_v < 1.0)
0149     initialHeight = 1.0;
0150   double heightIt = initialHeight;
0151   int heightbins = height / heightIt;  // round down
0152 
0153   edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0154       << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: roll with id = " << roll->id();
0155   edm::LogVerbatim("ME0PreRecoGaussianModelNoise") << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId()
0156                                                    << "] :: extracting parameters from the TrapezoidalStripTopology";
0157   edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0158       << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: bottom = " << bottomLength
0159       << " [cm] top  = " << topLength << " [cm] height = " << height << " [cm]"
0160       << " area  = " << trArea << " [cm^2] Rmid = " << rollRadius
0161       << " [cm] => Rmin = " << rollRadius - height * 1.0 / 2.0 << " [cm] Rmax = " << rollRadius + height * 1.0 / 2.0
0162       << " [cm]";
0163   edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0164       << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: heightbins = " << heightbins;
0165 
0166   for (int hx = 0; hx < heightbins; ++hx) {
0167     double bottomIt = bottomLength + hx * 2 * tan(10. / 180 * 3.14) * heightIt;
0168     double topIt = bottomLength + (hx + 1) * 2 * tan(10. / 180 * 3.14) * heightIt;
0169     if (hx == heightbins - 1) {
0170       topIt = topLength;  // last bin ... make strip a bit larger to cover entire roll
0171       heightIt = height - hx * heightIt;
0172     }
0173     double areaIt = heightIt * (bottomIt + topIt) * 1.0 / 2;
0174 
0175     edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0176         << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: height = " << std::setw(12) << heightIt
0177         << " [cm] bottom = " << std::setw(12) << bottomIt << " [cm]"
0178         << " top = " << std::setw(12) << topIt << " [cm] area = " << std::setw(12) << areaIt
0179         << " [cm^2] || sin(10) = " << sin(10. / 180 * 3.14);
0180 
0181     double myRandY = CLHEP::RandFlat::shoot(engine);
0182     double y0_rand = (hx + myRandY) * heightIt;  // Y coord, measured from the bottom of the roll
0183     if (hx == heightbins - 1)
0184       y0_rand = hx * initialHeight + myRandY * heightIt;
0185     double yy_rand =
0186         (y0_rand -
0187          height * 1.0 / 2);  // Y coord, measured from the middle of the roll, which is the Y coord in Local Coords
0188     double yy_glob = rollRadius + yy_rand;  // R coord in Global Coords
0189     // max length in x for given y coordinate (cfr trapezoidal eta partition)
0190     const float rSqrtR = yy_glob * sqrt(yy_glob);
0191     double xMax = topLength / 2.0 - (height / 2.0 - yy_rand) * myTanPhi;
0192 
0193     double sigma_u_new = sigma_u;
0194     if (constPhiSmearing_)
0195       sigma_u_new = correctSigmaU(roll, yy_rand);
0196 
0197     // 1) Intrinsic Noise ... Not implemented right now
0198     // ------------------------------------------------
0199     // if (simulateIntrinsicNoise_)
0200     // {
0201     // }
0202 
0203     // 2) Background Noise
0204     // ----------------------------
0205 
0206     // 2a) electron background
0207     // -----------------------
0208     if (simulateElectronBkg_) {
0209       // Extract / Calculate the Average Electron Rate
0210       // for the given global Y coord from Parametrization
0211       double averageElectronRatePerRoll =
0212           eleBkg[0] * rSqrtR * std::exp(eleBkg[1] / rSqrtR) + eleBkg[2] / rSqrtR + eleBkg[3] / (sqrt(yy_glob));
0213       // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
0214       averageElectronRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;
0215 
0216       // Rate [Hz/cm^2] * Nbx * 25*10^-9 [s] * Area [cm] = # hits in this roll in this bx
0217       const double averageElecRate(averageElectronRatePerRoll * (maxBunch_ - minBunch_ + 1) * (bxwidth * 1.0e-9) *
0218                                    areaIt);
0219 
0220       edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0221           << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "]" /* "] :: BX = "<<std::showpos<<bx*/
0222           << " evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageElectronRatePerRoll
0223           << " [Hz/cm^2]"
0224           << " x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt << " [cm^2]) ==> " << std::setw(12)
0225           << averageElecRate << " [hits]";
0226 
0227       bool ele_eff = (CLHEP::RandFlat::shoot(engine) < averageElecRate) ? true : false;
0228 
0229       edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0230           << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "] :: myRandY = " << std::setw(12) << myRandY
0231           << " => local y = " << std::setw(12) << yy_rand << " [cm]"
0232           << " => global y (global R) = " << std::setw(12) << yy_glob << " [cm] || Probability = " << std::setw(12)
0233           << averageElecRate << " => efficient? " << ele_eff << std::endl;
0234 
0235       if (ele_eff) {
0236         //calculate xx_rand at a given yy_rand
0237         double myRandX = CLHEP::RandFlat::shoot(engine);
0238         double xx_rand = 2 * xMax * (myRandX - 0.5);
0239         double ex = sigma_u_new;
0240         double ey = sigma_v;
0241         double corr = 0.;
0242         // extract random BX
0243         double myrandBX = CLHEP::RandFlat::shoot(engine);
0244         int bx = int((maxBunch_ - minBunch_ + 1) * myrandBX) + minBunch_;
0245         // extract random time in this BX
0246         double myrandT = CLHEP::RandFlat::shoot(engine);
0247         double minBXtime = (bx - 0.5) * bxwidth;  // double maxBXtime = (bx+0.5)*bxwidth;
0248         double time = myrandT * bxwidth + minBXtime;
0249         double myrandP = CLHEP::RandFlat::shoot(engine);
0250         int pdgid = 0;
0251         if (myrandP <= 0.5)
0252           pdgid = -11;  // electron
0253         else
0254           pdgid = 11;  // positron
0255         if (ex == 0)
0256           ex = error_u;  //errors cannot be zero
0257         if (ey == 0)
0258           ey = error_v;
0259 
0260         digi_.emplace(xx_rand, yy_rand, ex, ey, corr, time, pdgid, 0);
0261 
0262         edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0263             << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "] =====> electron hit in " << roll->id()
0264             << " pdgid = " << pdgid << " bx = " << bx << " ==> digitized"
0265             << " at loc x = " << xx_rand << " loc y = " << yy_rand << " time = " << time << " [ns]";
0266       }
0267     }  // end if electron bkg
0268 
0269     // 2b) neutral (n+g) background
0270     // ----------------------------
0271     if (simulateNeutralBkg_) {
0272       // Extract / Calculate the Average Neutral Rate
0273       // for the given global Y coord from Parametrization
0274       double averageNeutralRatePerRoll =
0275           neuBkg[0] * yy_glob * std::exp(neuBkg[1] / rSqrtR) + neuBkg[2] / rSqrtR + neuBkg[3] / (sqrt(yy_glob));
0276       // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
0277       averageNeutralRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;
0278 
0279       // Rate [Hz/cm^2] * Nbx * 25*10^-9 [s] * Area [cm] = # hits in this roll
0280       const double averageNeutrRate(averageNeutralRatePerRoll * (maxBunch_ - minBunch_ + 1) * (bxwidth * 1.0e-9) *
0281                                     areaIt);
0282 
0283       edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0284           << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "]" /* "] :: BX = "<<std::showpos<<bx*/
0285           << " evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageNeutralRatePerRoll
0286           << " [Hz/cm^2]"
0287           << " x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt << " [cm^2]) ==> " << std::setw(12)
0288           << averageNeutrRate << " [hits]";
0289 
0290       bool neu_eff = (CLHEP::RandFlat::shoot(engine) < averageNeutrRate) ? true : false;
0291 
0292       edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0293           << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "] :: myRandY = " << std::setw(12) << myRandY
0294           << " => local y = " << std::setw(12) << yy_rand << " [cm]"
0295           << " => global y (global R) = " << std::setw(12) << yy_glob << " [cm] || Probability = " << std::setw(12)
0296           << averageNeutrRate << " => efficient? " << neu_eff << std::endl;
0297 
0298       if (neu_eff) {
0299         //calculate xx_rand at a given yy_rand
0300         double myRandX = CLHEP::RandFlat::shoot(engine);
0301         double xx_rand = 2 * xMax * (myRandX - 0.5);
0302         double ex = sigma_u_new;
0303         double ey = sigma_v;
0304         double corr = 0.;
0305         // extract random BX
0306         double myrandBX = CLHEP::RandFlat::shoot(engine);
0307         int bx = int((maxBunch_ - minBunch_ + 1) * myrandBX) + minBunch_;
0308         // extract random time in this BX
0309         double myrandT = CLHEP::RandFlat::shoot(engine);
0310         double minBXtime = (bx - 0.5) * bxwidth;
0311         double time = myrandT * bxwidth + minBXtime;
0312         int pdgid = 0;
0313         double myrandP = CLHEP::RandFlat::shoot(engine);
0314         if (myrandP <= 0.08)
0315           pdgid = 2112;  // neutrons: GEM sensitivity for neutrons: 0.08%
0316         else
0317           pdgid =
0318               22;  // photons:  GEM sensitivity for photons:  1.04% ==> neutron fraction = (0.08 / 1.04) = 0.077 = 0.08
0319         if (ex == 0)
0320           ex = error_u;  //errors cannot be zero
0321         if (ey == 0)
0322           ey = error_v;
0323 
0324         digi_.emplace(xx_rand, yy_rand, ex, ey, corr, time, pdgid, 0);
0325 
0326         edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
0327             << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "] ======> neutral hit in " << roll->id()
0328             << " pdgid = " << pdgid << " bx = " << bx << " ==> digitized"
0329             << " at loc x = " << xx_rand << " loc y = " << yy_rand << " time = " << time << " [ns]";
0330       }
0331 
0332     }  // end if neutral bkg
0333   }    // end loop over strips (= pseudo rolls)
0334 }
0335 
0336 double ME0PreRecoGaussianModel::correctSigmaU(const ME0EtaPartition* roll, double y) {
0337   const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
0338   auto& parameters(roll->specs()->parameters());
0339   double height(parameters[2]);        // height     = height from Center of Roll
0340   double rollRadius = top_->radius();  // rollRadius = Radius at Center of Roll
0341   double Rmax = rollRadius + height;   // MaxRadius  = Radius at top of Roll
0342   double Rx = rollRadius + y;          // y in [-height,+height]
0343   double sigma_u_new = Rx / Rmax * sigma_u;
0344   return sigma_u_new;
0345 }
[ 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!