/
githubmirror
/
cmssw
Обзор
Документация
Войти
/
githubmirror
/
cmssw
Код
Запросы
0
Пакеты
0
Релизы
0
Аналитика
Безопасность
master
SimMuon/GEMDigitizer/src/ME0PreRecoGaussianModel.cc
345 строк
16 KB
Shahzad Malik Muzaffar
[SIMULATION-UPGRADE] Changes suggested by new llvm18 clang-format
05 сен 2024, 17:11
05 сен 2024, 17:11
07d786d
Код
Авторство
О чём код?
#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; }