CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
ME0PreRecoGaussianModel.cc
Go to the documentation of this file.
1 #include "SimMuon/GEMDigitizer/interface/ME0PreRecoGaussianModel.h"
2 #include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h"
3 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
4 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
5 
6 #include "FWCore/ServiceRegistry/interface/Service.h"
7 #include "FWCore/MessageLogger/interface/MessageLogger.h"
8 
9 #include "CLHEP/Random/RandFlat.h"
10 #include "CLHEP/Random/RandGaussQ.h"
11 #include "CLHEP/Random/RandPoissonQ.h"
12 
13 #include <cmath>
14 #include <utility>
15 #include <map>
16 
17 const int bxwidth = 25; // [ns]
18 
19 ME0PreRecoGaussianModel::ME0PreRecoGaussianModel(const edm::ParameterSet& config)
20  : ME0DigiPreRecoModel(config),
21  sigma_t(config.getParameter<double>("timeResolution")),
22  sigma_u(config.getParameter<double>("phiResolution")),
23  sigma_v(config.getParameter<double>("etaResolution")),
24  error_u(config.getParameter<double>("phiError")),
25  error_v(config.getParameter<double>("etaError")),
26  gaussianSmearing_(config.getParameter<bool>("gaussianSmearing")),
27  constPhiSmearing_(config.getParameter<bool>("constantPhiSpatialResolution")),
28  corr(config.getParameter<bool>("useCorrelation")),
29  etaproj(config.getParameter<bool>("useEtaProjectiveGEO")),
30  digitizeOnlyMuons_(config.getParameter<bool>("digitizeOnlyMuons")),
31  averageEfficiency_(config.getParameter<double>("averageEfficiency")),
32  // simulateIntrinsicNoise_(config.getParameter<bool>("simulateIntrinsicNoise")),
33  // averageNoiseRate_(config.getParameter<double>("averageNoiseRate")),
34  simulateElectronBkg_(config.getParameter<bool>("simulateElectronBkg")),
35  simulateNeutralBkg_(config.getParameter<bool>("simulateNeutralBkg")),
36  minBunch_(config.getParameter<int>("minBunch")),
37  maxBunch_(config.getParameter<int>("maxBunch")),
38  instLumi_(config.getParameter<double>("instLumi")),
39  rateFact_(config.getParameter<double>("rateFact")),
40  referenceInstLumi_(config.getParameter<double>("referenceInstLumi")) {
41  // polynomial parametrisation of neutral (n+g) and electron background
42  // This is the background for an Instantaneous Luminosity of L = 5E34 cm^-2 s^-1
43  neuBkg.push_back(0.00386257);
44  neuBkg.push_back(6344.65);
45  neuBkg.push_back(16627700);
46  neuBkg.push_back(-102098);
47 
48  eleBkg.push_back(0.00171409);
49  eleBkg.push_back(4900.56);
50  eleBkg.push_back(710909);
51  eleBkg.push_back(-4327.25);
52 }
53 
54 ME0PreRecoGaussianModel::~ME0PreRecoGaussianModel() {}
55 
56 void ME0PreRecoGaussianModel::simulateSignal(const ME0EtaPartition* roll,
57  const edm::PSimHitContainer& simHits,
58  CLHEP::HepRandomEngine* engine) {
59  for (const auto& hit : simHits) {
60  // Digitize only Muons?
61  if (std::abs(hit.particleType()) != 13 && digitizeOnlyMuons_)
62  continue;
63  // Digitize only in [minBunch,maxBunch] window
64  // window is: [(2n+1)*bxw/2, (2n+3)*bxw/2], n = [minBunch, maxBunch]
65  if (hit.timeOfFlight() < (2 * minBunch_ + 1) * bxwidth * 1.0 / 2 ||
66  hit.timeOfFlight() > (2 * maxBunch_ + 3) * bxwidth * 1.0 / 2)
67  continue;
68  // is GEM efficient?
69  if (CLHEP::RandFlat::shoot(engine) > averageEfficiency_)
70  continue;
71  // create digi
72  auto entry = hit.entryPoint();
73  double x = 0.0, y = 0.0;
74 
75  double sigma_u_new = sigma_u;
76  if (constPhiSmearing_)
77  sigma_u_new = correctSigmaU(roll, entry.y());
78 
79  if (gaussianSmearing_) { // Gaussian Smearing
80  x = CLHEP::RandGaussQ::shoot(engine, entry.x(), sigma_u_new);
81  y = CLHEP::RandGaussQ::shoot(engine, entry.y(), sigma_v);
82  } else { // Uniform Smearing ... use the sigmas as boundaries
83  x = entry.x() + (CLHEP::RandFlat::shoot(engine) - 0.5) * sigma_u_new;
84  y = entry.y() + (CLHEP::RandFlat::shoot(engine) - 0.5) * sigma_v;
85  }
86  double ex = sigma_u_new;
87  double ey = sigma_v;
88  double corr = 0.;
89  double tof = CLHEP::RandGaussQ::shoot(engine, hit.timeOfFlight(), sigma_t);
90  int pdgid = hit.particleType();
91  if (ex == 0)
92  ex = error_u; //errors cannot be zero
93  if (ey == 0)
94  ey = error_v;
95 
96  int evtId = hit.eventId().event();
97  int bx = hit.eventId().bunchCrossing();
98  int procType = hit.processType();
99  int res = 1;
100  if (!(evtId == 0 && bx == 0 && procType == 0))
101  res = 2;
102 
103  digi_.emplace(x, y, ex, ey, corr, tof, pdgid, res);
104 
105  edm::LogVerbatim("ME0PreRecoGaussianModel")
106  << "[ME0PreRecoDigi :: simulateSignal] :: simhit in " << roll->id() << " at loc x = " << std::setw(8)
107  << entry.x() << " [cm]"
108  << " loc y = " << std::setw(8) << entry.y() << " [cm] time = " << std::setw(8) << hit.timeOfFlight()
109  << " [ns] pdgid = " << std::showpos << std::setw(4) << pdgid;
110  edm::LogVerbatim("ME0PreRecoGaussianModel")
111  << "[ME0PreRecoDigi :: simulateSignal] :: digi in " << roll->id() << " at loc x = " << std::setw(8) << x
112  << " [cm] loc y = " << std::setw(8) << y << " [cm]"
113  << " time = " << std::setw(8) << tof << " [ns]";
114  edm::LogVerbatim("ME0PreRecoGaussianModel")
115  << "[ME0PreRecoDigi :: simulateSignal] :: digi in " << roll->id() << " with DX = " << std::setw(8)
116  << (entry.x() - x) << " [cm]"
117  << " DY = " << std::setw(8) << (entry.y() - y) << " [cm] DT = " << std::setw(8) << (hit.timeOfFlight() - tof)
118  << " [ns]";
119  }
120 }
121 
122 void ME0PreRecoGaussianModel::simulateNoise(const ME0EtaPartition* roll, CLHEP::HepRandomEngine* engine) {
123  double trArea(0.0);
124  const ME0DetId me0Id(roll->id());
125 
126  // Extract detailed information from the Strip Topology:
127  // base_bottom, base_top, height, strips, pads
128  // note that (0,0) is in the middle of the roll ==> all param are at all half length
129  if (me0Id.region() == 0) {
130  throw cms::Exception("Geometry") << "Asking TrapezoidalStripTopology from a ME0 will fail";
131  } // not sure we really need this
132  const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
133 
134  auto& parameters(roll->specs()->parameters());
135  double bottomLength(parameters[0]);
136  bottomLength = 2 * bottomLength; // bottom is largest length, so furtest away from beamline
137  double topLength(parameters[1]);
138  topLength = 2 * topLength; // top is shortest length, so closest to beamline
139  double height(parameters[2]);
140  height = 2 * height;
141  double myTanPhi = (topLength - bottomLength) / (height * 2);
142  double rollRadius = top_->radius();
143  trArea = height * (topLength + bottomLength) / 2.0;
144 
145  // Divide the detector area in different strips
146  // take smearing in y-coord as height for each strip
147  double initialHeight = sigma_v;
148  if (sigma_v < 1.0)
149  initialHeight = 1.0;
150  double heightIt = initialHeight;
151  int heightbins = height / heightIt; // round down
152 
153  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
154  << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: roll with id = " << roll->id();
155  edm::LogVerbatim("ME0PreRecoGaussianModelNoise") << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId()
156  << "] :: extracting parameters from the TrapezoidalStripTopology";
157  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
158  << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: bottom = " << bottomLength
159  << " [cm] top = " << topLength << " [cm] height = " << height << " [cm]"
160  << " area = " << trArea << " [cm^2] Rmid = " << rollRadius
161  << " [cm] => Rmin = " << rollRadius - height * 1.0 / 2.0 << " [cm] Rmax = " << rollRadius + height * 1.0 / 2.0
162  << " [cm]";
163  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
164  << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: heightbins = " << heightbins;
165 
166  for (int hx = 0; hx < heightbins; ++hx) {
167  double bottomIt = bottomLength + hx * 2 * tan(10. / 180 * 3.14) * heightIt;
168  double topIt = bottomLength + (hx + 1) * 2 * tan(10. / 180 * 3.14) * heightIt;
169  if (hx == heightbins - 1) {
170  topIt = topLength; // last bin ... make strip a bit larger to cover entire roll
171  heightIt = height - hx * heightIt;
172  }
173  double areaIt = heightIt * (bottomIt + topIt) * 1.0 / 2;
174 
175  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
176  << "[ME0PreRecoDigi :: sNoise][" << roll->id().rawId() << "] :: height = " << std::setw(12) << heightIt
177  << " [cm] bottom = " << std::setw(12) << bottomIt << " [cm]"
178  << " top = " << std::setw(12) << topIt << " [cm] area = " << std::setw(12) << areaIt
179  << " [cm^2] || sin(10) = " << sin(10. / 180 * 3.14);
180 
181  double myRandY = CLHEP::RandFlat::shoot(engine);
182  double y0_rand = (hx + myRandY) * heightIt; // Y coord, measured from the bottom of the roll
183  if (hx == heightbins - 1)
184  y0_rand = hx * initialHeight + myRandY * heightIt;
185  double yy_rand =
186  (y0_rand -
187  height * 1.0 / 2); // Y coord, measured from the middle of the roll, which is the Y coord in Local Coords
188  double yy_glob = rollRadius + yy_rand; // R coord in Global Coords
189  // max length in x for given y coordinate (cfr trapezoidal eta partition)
190  const float rSqrtR = yy_glob * sqrt(yy_glob);
191  double xMax = topLength / 2.0 - (height / 2.0 - yy_rand) * myTanPhi;
192 
193  double sigma_u_new = sigma_u;
194  if (constPhiSmearing_)
195  sigma_u_new = correctSigmaU(roll, yy_rand);
196 
197  // 1) Intrinsic Noise ... Not implemented right now
198  // ------------------------------------------------
199  // if (simulateIntrinsicNoise_)
200  // {
201  // }
202 
203  // 2) Background Noise
204  // ----------------------------
205 
206  // 2a) electron background
207  // -----------------------
208  if (simulateElectronBkg_) {
209  // Extract / Calculate the Average Electron Rate
210  // for the given global Y coord from Parametrization
211  double averageElectronRatePerRoll =
212  eleBkg[0] * rSqrtR * std::exp(eleBkg[1] / rSqrtR) + eleBkg[2] / rSqrtR + eleBkg[3] / (sqrt(yy_glob));
213  // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
214  averageElectronRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;
215 
216  // Rate [Hz/cm^2] * Nbx * 25*10^-9 [s] * Area [cm] = # hits in this roll in this bx
217  const double averageElecRate(averageElectronRatePerRoll * (maxBunch_ - minBunch_ + 1) * (bxwidth * 1.0e-9) *
218  areaIt);
219 
220  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
221  << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "]" /* "] :: BX = "<<std::showpos<<bx*/
222  << " evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageElectronRatePerRoll
223  << " [Hz/cm^2]"
224  << " x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt << " [cm^2]) ==> " << std::setw(12)
225  << averageElecRate << " [hits]";
226 
227  bool ele_eff = (CLHEP::RandFlat::shoot(engine) < averageElecRate) ? true : false;
228 
229  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
230  << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "] :: myRandY = " << std::setw(12) << myRandY
231  << " => local y = " << std::setw(12) << yy_rand << " [cm]"
232  << " => global y (global R) = " << std::setw(12) << yy_glob << " [cm] || Probability = " << std::setw(12)
233  << averageElecRate << " => efficient? " << ele_eff << std::endl;
234 
235  if (ele_eff) {
236  //calculate xx_rand at a given yy_rand
237  double myRandX = CLHEP::RandFlat::shoot(engine);
238  double xx_rand = 2 * xMax * (myRandX - 0.5);
239  double ex = sigma_u_new;
240  double ey = sigma_v;
241  double corr = 0.;
242  // extract random BX
243  double myrandBX = CLHEP::RandFlat::shoot(engine);
244  int bx = int((maxBunch_ - minBunch_ + 1) * myrandBX) + minBunch_;
245  // extract random time in this BX
246  double myrandT = CLHEP::RandFlat::shoot(engine);
247  double minBXtime = (bx - 0.5) * bxwidth; // double maxBXtime = (bx+0.5)*bxwidth;
248  double time = myrandT * bxwidth + minBXtime;
249  double myrandP = CLHEP::RandFlat::shoot(engine);
250  int pdgid = 0;
251  if (myrandP <= 0.5)
252  pdgid = -11; // electron
253  else
254  pdgid = 11; // positron
255  if (ex == 0)
256  ex = error_u; //errors cannot be zero
257  if (ey == 0)
258  ey = error_v;
259 
260  digi_.emplace(xx_rand, yy_rand, ex, ey, corr, time, pdgid, 0);
261 
262  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
263  << "[ME0PreRecoDigi :: elebkg][" << roll->id().rawId() << "] =====> electron hit in " << roll->id()
264  << " pdgid = " << pdgid << " bx = " << bx << " ==> digitized"
265  << " at loc x = " << xx_rand << " loc y = " << yy_rand << " time = " << time << " [ns]";
266  }
267  } // end if electron bkg
268 
269  // 2b) neutral (n+g) background
270  // ----------------------------
271  if (simulateNeutralBkg_) {
272  // Extract / Calculate the Average Neutral Rate
273  // for the given global Y coord from Parametrization
274  double averageNeutralRatePerRoll =
275  neuBkg[0] * yy_glob * std::exp(neuBkg[1] / rSqrtR) + neuBkg[2] / rSqrtR + neuBkg[3] / (sqrt(yy_glob));
276  // Scale up/down for desired instantaneous lumi (reference is 5E34, double from config is in units of 1E34)
277  averageNeutralRatePerRoll *= instLumi_ * rateFact_ * 1.0 / referenceInstLumi_;
278 
279  // Rate [Hz/cm^2] * Nbx * 25*10^-9 [s] * Area [cm] = # hits in this roll
280  const double averageNeutrRate(averageNeutralRatePerRoll * (maxBunch_ - minBunch_ + 1) * (bxwidth * 1.0e-9) *
281  areaIt);
282 
283  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
284  << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "]" /* "] :: BX = "<<std::showpos<<bx*/
285  << " evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageNeutralRatePerRoll
286  << " [Hz/cm^2]"
287  << " x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt << " [cm^2]) ==> " << std::setw(12)
288  << averageNeutrRate << " [hits]";
289 
290  bool neu_eff = (CLHEP::RandFlat::shoot(engine) < averageNeutrRate) ? true : false;
291 
292  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
293  << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "] :: myRandY = " << std::setw(12) << myRandY
294  << " => local y = " << std::setw(12) << yy_rand << " [cm]"
295  << " => global y (global R) = " << std::setw(12) << yy_glob << " [cm] || Probability = " << std::setw(12)
296  << averageNeutrRate << " => efficient? " << neu_eff << std::endl;
297 
298  if (neu_eff) {
299  //calculate xx_rand at a given yy_rand
300  double myRandX = CLHEP::RandFlat::shoot(engine);
301  double xx_rand = 2 * xMax * (myRandX - 0.5);
302  double ex = sigma_u_new;
303  double ey = sigma_v;
304  double corr = 0.;
305  // extract random BX
306  double myrandBX = CLHEP::RandFlat::shoot(engine);
307  int bx = int((maxBunch_ - minBunch_ + 1) * myrandBX) + minBunch_;
308  // extract random time in this BX
309  double myrandT = CLHEP::RandFlat::shoot(engine);
310  double minBXtime = (bx - 0.5) * bxwidth;
311  double time = myrandT * bxwidth + minBXtime;
312  int pdgid = 0;
313  double myrandP = CLHEP::RandFlat::shoot(engine);
314  if (myrandP <= 0.08)
315  pdgid = 2112; // neutrons: GEM sensitivity for neutrons: 0.08%
316  else
317  pdgid =
318  22; // photons: GEM sensitivity for photons: 1.04% ==> neutron fraction = (0.08 / 1.04) = 0.077 = 0.08
319  if (ex == 0)
320  ex = error_u; //errors cannot be zero
321  if (ey == 0)
322  ey = error_v;
323 
324  digi_.emplace(xx_rand, yy_rand, ex, ey, corr, time, pdgid, 0);
325 
326  edm::LogVerbatim("ME0PreRecoGaussianModelNoise")
327  << "[ME0PreRecoDigi :: neubkg][" << roll->id().rawId() << "] ======> neutral hit in " << roll->id()
328  << " pdgid = " << pdgid << " bx = " << bx << " ==> digitized"
329  << " at loc x = " << xx_rand << " loc y = " << yy_rand << " time = " << time << " [ns]";
330  }
331 
332  } // end if neutral bkg
333  } // end loop over strips (= pseudo rolls)
334 }
335 
336 double ME0PreRecoGaussianModel::correctSigmaU(const ME0EtaPartition* roll, double y) {
337  const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology())));
338  auto& parameters(roll->specs()->parameters());
339  double height(parameters[2]); // height = height from Center of Roll
340  double rollRadius = top_->radius(); // rollRadius = Radius at Center of Roll
341  double Rmax = rollRadius + height; // MaxRadius = Radius at top of Roll
342  double Rx = rollRadius + y; // y in [-height,+height]
343  double sigma_u_new = Rx / Rmax * sigma_u;
344  return sigma_u_new;
345 }
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
ME0EtaPartition::specs
const ME0EtaPartitionSpecs * specs() const
Definition: ME0EtaPartition.h:17
ME0PreRecoGaussianModel::neuBkg
std::vector< double > neuBkg
Definition: ME0PreRecoGaussianModel.h:53
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
ME0EtaPartitionSpecs::parameters
const ME0Specs & parameters() const
Definition: ME0EtaPartitionSpecs.cc:47
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
ME0DigiPreRecoModel::digi_
std::set< ME0DigiPreReco > digi_
Definition: ME0DigiPreRecoModel.h:48
ME0PreRecoGaussianModel::correctSigmaU
double correctSigmaU(const ME0EtaPartition *, double)
Definition: ME0PreRecoGaussianModel.cc:336
bxwidth
const int bxwidth
Definition: ME0PreRecoGaussianModel.cc:17
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
ME0EtaPartition::id
ME0DetId id() const
Definition: ME0EtaPartition.h:18
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ME0PreRecoGaussianModel::simulateSignal
void simulateSignal(const ME0EtaPartition *, const edm::PSimHitContainer &, CLHEP::HepRandomEngine *) override
Definition: ME0PreRecoGaussianModel.cc:56
trackerHits.simHits
tuple simHits
Definition: trackerHits.py:16
ME0PreRecoGaussianModel::simulateNoise
void simulateNoise(const ME0EtaPartition *, CLHEP::HepRandomEngine *) override
Definition: ME0PreRecoGaussianModel.cc:122
ME0PreRecoGaussianModel::ME0PreRecoGaussianModel
ME0PreRecoGaussianModel(const edm::ParameterSet &)
Definition: ME0PreRecoGaussianModel.cc:19
submitPVResolutionJobs.config
tuple config
parse the configuration file
Definition: submitPVResolutionJobs.py:281
ME0EtaPartition::topology
const Topology & topology() const override
Definition: ME0EtaPartition.cc:14
mps_splice.entry
list entry
Definition: mps_splice.py:68
Gflash::Rmax
const double Rmax[kNumberCalorimeter]
Definition: GflashNameSpace.h:31
edm::PSimHitContainer
std::vector< PSimHit > PSimHitContainer
Definition: PSimHitContainer.h:11
ME0PreRecoGaussianModel::eleBkg
std::vector< double > eleBkg
Definition: ME0PreRecoGaussianModel.h:53
Generated for CMSSW Reference Manual by doxygen 1.8.5