9 #include "CLHEP/Random/RandFlat.h"
10 #include "CLHEP/Random/RandGaussQ.h"
11 #include "CLHEP/Random/RandPoissonQ.h"
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")),
29 etaproj(
config.getParameter<
bool>(
"useEtaProjectiveGEO")),
30 digitizeOnlyMuons_(
config.getParameter<
bool>(
"digitizeOnlyMuons")),
31 averageEfficiency_(
config.getParameter<double>(
"averageEfficiency")),
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")) {
43 neuBkg.push_back(0.00386257);
45 neuBkg.push_back(16627700);
48 eleBkg.push_back(0.00171409);
51 eleBkg.push_back(-4327.25);
58 CLHEP::HepRandomEngine* engine) {
73 double x = 0.0,
y = 0.0;
80 x = CLHEP::RandGaussQ::shoot(engine,
entry.x(), sigma_u_new);
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;
86 double ex = sigma_u_new;
89 double tof = CLHEP::RandGaussQ::shoot(engine,
hit.timeOfFlight(),
sigma_t);
96 int evtId =
hit.eventId().event();
97 int bx =
hit.eventId().bunchCrossing();
98 int procType =
hit.processType();
100 if (!(evtId == 0 &&
bx == 0 && procType == 0))
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;
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]";
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)
129 if (me0Id.region() == 0) {
130 throw cms::Exception(
"Geometry") <<
"Asking TrapezoidalStripTopology from a ME0 will fail";
136 bottomLength = 2 * bottomLength;
138 topLength = 2 * topLength;
141 double myTanPhi = (topLength - bottomLength) / (height * 2);
142 double rollRadius = top_->radius();
143 trArea = height * (topLength + bottomLength) / 2.0;
147 double initialHeight =
sigma_v;
150 double heightIt = initialHeight;
151 int heightbins = height / heightIt;
154 <<
"[ME0PreRecoDigi :: sNoise][" << roll->
id().
rawId() <<
"] :: roll with id = " << roll->
id();
156 <<
"] :: extracting parameters from the TrapezoidalStripTopology";
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
164 <<
"[ME0PreRecoDigi :: sNoise][" << roll->
id().
rawId() <<
"] :: heightbins = " << heightbins;
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) {
171 heightIt = height - hx * heightIt;
173 double areaIt = heightIt * (bottomIt + topIt) * 1.0 / 2;
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);
181 double myRandY = CLHEP::RandFlat::shoot(engine);
182 double y0_rand = (hx + myRandY) * heightIt;
183 if (hx == heightbins - 1)
184 y0_rand = hx * initialHeight + myRandY * heightIt;
188 double yy_glob = rollRadius + yy_rand;
190 const float rSqrtR = yy_glob *
sqrt(yy_glob);
191 double xMax = topLength / 2.0 - (height / 2.0 - yy_rand) * myTanPhi;
211 double averageElectronRatePerRoll =
221 <<
"[ME0PreRecoDigi :: elebkg][" << roll->
id().
rawId() <<
"]"
222 <<
" evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageElectronRatePerRoll
224 <<
" x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt <<
" [cm^2]) ==> " << std::setw(12)
225 << averageElecRate <<
" [hits]";
227 bool ele_eff = (CLHEP::RandFlat::shoot(engine) < averageElecRate) ?
true :
false;
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;
237 double myRandX = CLHEP::RandFlat::shoot(engine);
238 double xx_rand = 2 *
xMax * (myRandX - 0.5);
239 double ex = sigma_u_new;
243 double myrandBX = CLHEP::RandFlat::shoot(engine);
246 double myrandT = CLHEP::RandFlat::shoot(engine);
249 double myrandP = CLHEP::RandFlat::shoot(engine);
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]";
274 double averageNeutralRatePerRoll =
284 <<
"[ME0PreRecoDigi :: neubkg][" << roll->
id().
rawId() <<
"]"
285 <<
" evaluation of Background Hit Rate at this coord :: " << std::setw(12) << averageNeutralRatePerRoll
287 <<
" x 9 x 25*10^-9 [s] x Area (of strip = " << std::setw(12) << areaIt <<
" [cm^2]) ==> " << std::setw(12)
288 << averageNeutrRate <<
" [hits]";
290 bool neu_eff = (CLHEP::RandFlat::shoot(engine) < averageNeutrRate) ?
true :
false;
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;
300 double myRandX = CLHEP::RandFlat::shoot(engine);
301 double xx_rand = 2 *
xMax * (myRandX - 0.5);
302 double ex = sigma_u_new;
306 double myrandBX = CLHEP::RandFlat::shoot(engine);
309 double myrandT = CLHEP::RandFlat::shoot(engine);
313 double myrandP = CLHEP::RandFlat::shoot(engine);
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]";
340 double rollRadius = top_->radius();
341 double Rmax = rollRadius + height;
342 double Rx = rollRadius +
y;