9 #include "CLHEP/Random/RandFlat.h"
10 #include "CLHEP/Random/RandPoissonQ.h"
11 #include "CLHEP/Random/RandGaussQ.h"
23 const double COSMIC_PAR(37.62);
28 , averageEfficiency_(config.getParameter<double> (
"averageEfficiency"))
29 , averageShapingTime_(config.getParameter<double> (
"averageShapingTime"))
30 , timeResolution_(config.getParameter<double> (
"timeResolution"))
31 , timeJitter_(config.getParameter<double> (
"timeJitter"))
32 , averageNoiseRate_(config.getParameter<double> (
"averageNoiseRate"))
34 , clsParametrization_(config.getParameter<std::vector<double>>(
"clsParametrization"))
35 , signalPropagationSpeed_(config.getParameter<double> (
"signalPropagationSpeed"))
36 , cosmics_(config.getParameter<bool> (
"cosmics"))
37 , bxwidth_(config.getParameter<int> (
"bxwidth"))
38 , minBunch_(config.getParameter<int> (
"minBunch"))
39 , maxBunch_(config.getParameter<int> (
"maxBunch"))
40 , digitizeOnlyMuons_(config.getParameter<bool> (
"digitizeOnlyMuons"))
41 , doBkgNoise_(config.getParameter<bool> (
"doBkgNoise"))
42 , doNoiseCLS_(config.getParameter<bool> (
"doNoiseCLS"))
43 , fixedRollRadius_(config.getParameter<bool> (
"fixedRollRadius"))
44 , scaleLumi_(config.getParameter<double> (
"scaleLumi"))
45 , simulateElectronBkg_(config.getParameter<bool> (
"simulateElectronBkg"))
46 , constNeuGE11_(config.getParameter<double> (
"constNeuGE11"))
47 , slopeNeuGE11_(config.getParameter<double> (
"slopeNeuGE11"))
48 , GE21NeuBkgParams_(config.getParameter<std::vector<double>>(
"GE21NeuBkgParams"))
49 , GE11ElecBkgParams_(config.getParameter<std::vector<double>>(
"GE11ElecBkgParams"))
50 , GE21ElecBkgParams_(config.getParameter<std::vector<double>>(
"GE21ElecBkgParams"))
71 for (edm::PSimHitContainer::const_iterator
hit = simHits.begin();
hit != simHits.end(); ++
hit)
80 for (
auto & digi : cluster)
94 float randomJitterTime = CLHEP::RandGaussQ::shoot(engine, 0.,
timeJitter_);
101 throw cms::Exception(
"Geometry")<<
"GEMSimpleModel::getSimHitBx() - GEM simhit id does not match any GEM roll id: " <<
id <<
"\n";
105 if (roll->id().region() == 0)
107 throw cms::Exception(
"Geometry") <<
"GEMSimpleModel::getSimHitBx() - this GEM id is from barrel, which cannot happen: " << roll->id() <<
"\n";
111 const double cspeed = 299792458;
112 const int nstrips = roll->nstrips();
113 float middleStrip = nstrips/2.;
114 LocalPoint middleOfRoll = roll->centreOfStrip(middleStrip);
115 GlobalPoint globMiddleRol = roll->toGlobal(middleOfRoll);
116 double muRadius =
sqrt(globMiddleRol.
x()*globMiddleRol.
x() + globMiddleRol.
y()*globMiddleRol.
y() +globMiddleRol.
z()*globMiddleRol.
z());
117 double timeCalibrationOffset_ = (muRadius *1
e+9)/(cspeed*1
e+2);
119 const TrapezoidalStripTopology*
top(dynamic_cast<const TrapezoidalStripTopology*> (&(roll->topology())));
120 const float halfStripLength(0.5 *
top->stripLength());
121 const float distanceFromEdge(halfStripLength - simHitPos.y());
127 const float averagePropagationTime(distanceFromEdge / signalPropagationSpeedTrue);
129 float randomResolutionTime = CLHEP::RandGaussQ::shoot(engine, 0.,
timeResolution_);
131 const float simhitTime(tof +
averageShapingTime_ + randomResolutionTime + averagePropagationTime + randomJitterTime);
133 float referenceTime = 0.;
134 referenceTime = timeCalibrationOffset_ + halfStripLength / signalPropagationSpeedTrue +
averageShapingTime_;
135 const float timeDifference(
cosmics_ ? (simhitTime - referenceTime) / COSMIC_PAR : simhitTime - referenceTime);
138 bx =
static_cast<int> (std::round((timeDifference) /
bxwidth_));
141 const bool debug(
false);
144 std::cout <<
"checktime " <<
"bx = " << bx <<
"\tdeltaT = " << timeDifference <<
"\tsimT = " << simhitTime
145 <<
"\trefT = " << referenceTime <<
"\ttof = " << tof <<
"\tavePropT = " << averagePropagationTime
146 <<
"\taveRefPropT = " << halfStripLength / signalPropagationSpeedTrue << std::endl;
157 const int nstrips(roll->
nstrips());
159 double trStripArea(0.0);
161 if (gemId.region() == 0)
164 <<
"GEMSynchronizer::simulateNoise() - this GEM id is from barrel, which cannot happen.";
166 const TrapezoidalStripTopology* top_(dynamic_cast<const TrapezoidalStripTopology*> (&(roll->
topology())));
167 const float striplength(top_->stripLength());
168 trStripArea = (roll->
pitch()) * striplength;
169 trArea = trStripArea * nstrips;
173 top_->radius() + CLHEP::RandFlat::shoot(engine, -1.*top_->stripLength()/2., top_->stripLength()/2.));
176 double averageNeutralNoiseRatePerRoll = 0.;
177 double averageNoiseElectronRatePerRoll = 0.;
178 double averageNoiseRatePerRoll = 0.;
180 if(gemId.station() == 1)
193 averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
196 if(gemId.station() == 2 || gemId.station() == 3)
215 +
GE21ElecBkgParams_[6]*rollRadius*rollRadius*rollRadius*rollRadius*rollRadius*rollRadius;
217 averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
224 for(
int j = 0;
j < nstrips; ++
j)
226 CLHEP::RandPoissonQ randPoissonQ(*engine, aveIntrinsicNoisePerStrip);
227 const int n_intrHits(randPoissonQ.fire());
229 for (
int k = 0;
k < n_intrHits;
k++ )
231 const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) +
minBunch_);
232 std::pair<int, int> digi(
k+1,time_hit);
239 const double averageNoise(averageNoiseRatePerRoll * nBxing *
bxwidth_ * trArea * 1.0
e-9 *
scaleLumi_);
240 CLHEP::RandPoissonQ randPoissonQ(*engine, averageNoise);
241 const int n_hits(randPoissonQ.fire());
243 for (
int i = 0;
i < n_hits; ++
i)
245 const int centralStrip(static_cast<int> (CLHEP::RandFlat::shoot(engine, 1, nstrips)));
246 const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) +
minBunch_);
250 std::vector<std::pair<int, int> > cluster_;
252 cluster_.push_back(std::pair<int, int>(centralStrip, time_hit));
255 double randForCls = CLHEP::RandFlat::shoot(engine);
277 if (clusterSize % 2 != 0)
279 int clsR = (clusterSize - 1) / 2;
280 for (
int i = 1;
i <= clsR; ++
i)
283 cluster_.push_back(std::pair<int, int>(centralStrip -
i, time_hit));
285 cluster_.push_back(std::pair<int, int>(centralStrip +
i, time_hit));
289 if (clusterSize % 2 == 0)
291 int clsR = (clusterSize - 2) / 2;
292 if(CLHEP::RandFlat::shoot(engine) < 0.5)
295 cluster_.push_back(std::pair<int, int>(centralStrip - 1, time_hit));
296 for (
int i = 1;
i <= clsR; ++
i)
299 cluster_.push_back(std::pair<int, int>(centralStrip - 1 -
i, time_hit));
301 cluster_.push_back(std::pair<int, int>(centralStrip +
i, time_hit));
307 if (CLHEP::RandFlat::shoot(engine) <
averageEfficiency_ && (centralStrip + 1 <= nstrips))
308 cluster_.push_back(std::pair<int, int>(centralStrip + 1, time_hit));
309 for (
int i = 1;
i <= clsR; ++
i)
311 if (CLHEP::RandFlat::shoot(engine) <
averageEfficiency_ && (centralStrip + 1 +
i <= nstrips))
312 cluster_.push_back(std::pair<int, int>(centralStrip + 1 +
i, time_hit));
314 cluster_.push_back(std::pair<int, int>(centralStrip -
i, time_hit));
318 for(
auto & digi : cluster_)
325 std::pair<int, int> digi(centralStrip, time_hit);
333 const PSimHit* simHit,
const int bx, CLHEP::HepRandomEngine* engine)
339 const int nstrips(roll->
nstrips());
341 int centralStrip = 0;
342 if (!(topology.
channel(hit_position) + 1 > nstrips))
343 centralStrip = topology.
channel(hit_position) + 1;
345 centralStrip = topology.
channel(hit_position);
349 double deltaphi = pointSimHit.
phi() - pointDigiHit.
phi();
352 std::vector<std::pair<int, int> > cluster_;
354 cluster_.push_back(std::pair<int, int>(centralStrip, bx));
358 double randForCls = CLHEP::RandFlat::shoot(engine);
383 if (clusterSize % 2 != 0)
385 int clsR = (clusterSize - 1) / 2;
386 for (
int i = 1;
i <= clsR; ++
i)
389 cluster_.push_back(std::pair<int, int>(centralStrip -
i, bx));
391 cluster_.push_back(std::pair<int, int>(centralStrip +
i, bx));
395 if (clusterSize % 2 == 0)
397 int clsR = (clusterSize - 2) / 2;
401 cluster_.push_back(std::pair<int, int>(centralStrip - 1, bx));
402 for (
int i = 1;
i <= clsR; ++
i)
405 cluster_.push_back(std::pair<int, int>(centralStrip - 1 -
i, bx));
407 cluster_.push_back(std::pair<int, int>(centralStrip +
i, bx));
412 if (CLHEP::RandFlat::shoot(engine) <
averageEfficiency_ && (centralStrip + 1 <= nstrips))
413 cluster_.push_back(std::pair<int, int>(centralStrip + 1, bx));
414 for (
int i = 1;
i <= clsR; ++
i)
416 if (CLHEP::RandFlat::shoot(engine) <
averageEfficiency_ && (centralStrip + 1 +
i <= nstrips))
417 cluster_.push_back(std::pair<int, int>(centralStrip + 1 +
i, bx));
419 cluster_.push_back(std::pair<int, int>(centralStrip -
i, bx));
std::vector< double > clsParametrization_
double signalPropagationSpeed_
CaloTopology const * topology(0)
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Geom::Phi< T > phi() const
const GEMEtaPartition * etaPartition(GEMDetId id) const
Return a etaPartition given its id.
uint32_t rawId() const
get the raw id
float timeOfFlight() const
bool simulateIntrinsicNoise_
std::vector< std::pair< int, int > > simulateClustering(const GEMEtaPartition *, const PSimHit *, const int, CLHEP::HepRandomEngine *engine) override
Local3DPoint localPosition() const
double averageEfficiency_
std::vector< double > GE11ElecBkgParams_
Abs< T >::type abs(const T &t)
virtual int channel(const LocalPoint &p) const =0
GEMSimpleModel(const edm::ParameterSet &)
float pabs() const
fast and more accurate access to momentumAtEntry().mag()
Container::value_type value_type
int getSimHitBx(const PSimHit *, CLHEP::HepRandomEngine *engine)
void simulateNoise(const GEMEtaPartition *, CLHEP::HepRandomEngine *engine) override
bool simulateElectronBkg_
void simulateSignal(const GEMEtaPartition *, const edm::PSimHitContainer &, CLHEP::HepRandomEngine *engine) override
const GEMGeometry * geometry_
double averageShapingTime_
std::set< std::pair< int, int > > strips_
std::vector< double > GE21ElecBkgParams_
DetectorHitMap detectorHitMap_
std::vector< PSimHit > PSimHitContainer
StripDigiSimLinks stripDigiSimLinks_
std::vector< double > GE21NeuBkgParams_
edm::DetSet< StripDigiSimLink > StripDigiSimLinks
unsigned int detUnitId() const