CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GEMSimpleModel.cc
Go to the documentation of this file.
2 
6 
8 
9 #include "CLHEP/Random/RandFlat.h"
10 #include "CLHEP/Random/RandPoissonQ.h"
11 #include "CLHEP/Random/RandGaussQ.h"
12 
13 #include <cmath>
14 #include <utility>
15 #include <map>
16 
17 #include "TMath.h" /* exp */
18 
19 
20 namespace
21 {
22  // "magic" parameter for cosmics
23  const double COSMIC_PAR(37.62);
24 }
25 
27 GEMDigiModel(config)
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"))
33 //, averageClusterSize_(config.getParameter<double> ("averageClusterSize"))
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"))
51 
52 {
53 
54 }
55 
57 {
58 }
59 
61 {
62  return;
63 }
64 
65 void GEMSimpleModel::simulateSignal(const GEMEtaPartition* roll, const edm::PSimHitContainer& simHits, CLHEP::HepRandomEngine* engine)
66 {
68  detectorHitMap_.clear();
70 
71  for (edm::PSimHitContainer::const_iterator hit = simHits.begin(); hit != simHits.end(); ++hit)
72  {
73  if (std::abs(hit->particleType()) != 13 && digitizeOnlyMuons_)
74  continue;
75 
76  // Check GEM efficiency
77  if (CLHEP::RandFlat::shoot(engine) > averageEfficiency_) continue;
78  const int bx(getSimHitBx(&(*hit), engine));
79  const std::vector<std::pair<int, int> > cluster(simulateClustering(roll, &(*hit), bx, engine));
80  for (auto & digi : cluster)
81  {
83  strips_.insert(digi);
84  }
85  }
86 }
87 
88 int GEMSimpleModel::getSimHitBx(const PSimHit* simhit, CLHEP::HepRandomEngine* engine)
89 {
90  int bx = -999;
91  const LocalPoint simHitPos(simhit->localPosition());
92  const float tof(simhit->timeOfFlight());
93  // random Gaussian time correction due to electronics jitter
94  float randomJitterTime = CLHEP::RandGaussQ::shoot(engine, 0., timeJitter_);
95 
96  const GEMDetId id(simhit->detUnitId());
97  const GEMEtaPartition* roll(geometry_->etaPartition(id));
98 
99  if (!roll)
100  {
101  throw cms::Exception("Geometry")<< "GEMSimpleModel::getSimHitBx() - GEM simhit id does not match any GEM roll id: " << id << "\n";
102  return 999;
103  }
104 
105  if (roll->id().region() == 0)
106  {
107  throw cms::Exception("Geometry") << "GEMSimpleModel::getSimHitBx() - this GEM id is from barrel, which cannot happen: " << roll->id() << "\n";
108  }
109 
110  // signal propagation speed in vacuum in [m/s]
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 *1e+9)/(cspeed*1e+2); //[ns]
118 
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());
122 
123  // signal propagation speed in material in [cm/ns]
124  double signalPropagationSpeedTrue = signalPropagationSpeed_ * cspeed * 1e+2 * 1e-9;
125 
126  // average time for the signal to propagate from the SimHit to the top of a strip
127  const float averagePropagationTime(distanceFromEdge / signalPropagationSpeedTrue);
128  // random Gaussian time correction due to the finite timing resolution of the detector
129  float randomResolutionTime = CLHEP::RandGaussQ::shoot(engine, 0., timeResolution_);
130 
131  const float simhitTime(tof + averageShapingTime_ + randomResolutionTime + averagePropagationTime + randomJitterTime);
132 
133  float referenceTime = 0.;
134  referenceTime = timeCalibrationOffset_ + halfStripLength / signalPropagationSpeedTrue + averageShapingTime_;
135  const float timeDifference(cosmics_ ? (simhitTime - referenceTime) / COSMIC_PAR : simhitTime - referenceTime);
136 
137  // assign the bunch crossing
138  bx = static_cast<int> (std::round((timeDifference) / bxwidth_));
139 
140  // check time
141  const bool debug(false);
142  if (debug)
143  {
144  std::cout << "checktime " << "bx = " << bx << "\tdeltaT = " << timeDifference << "\tsimT = " << simhitTime
145  << "\trefT = " << referenceTime << "\ttof = " << tof << "\tavePropT = " << averagePropagationTime
146  << "\taveRefPropT = " << halfStripLength / signalPropagationSpeedTrue << std::endl;
147  }
148  return bx;
149 }
150 
151 void GEMSimpleModel::simulateNoise(const GEMEtaPartition* roll, CLHEP::HepRandomEngine* engine)
152 {
153  if (!doBkgNoise_)
154  return;
155 
156  const GEMDetId gemId(roll->id());
157  const int nstrips(roll->nstrips());
158  double trArea(0.0);
159  double trStripArea(0.0);
160 
161  if (gemId.region() == 0)
162  {
163  throw cms::Exception("Geometry")
164  << "GEMSynchronizer::simulateNoise() - this GEM id is from barrel, which cannot happen.";
165  }
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;
170  const int nBxing(maxBunch_ - minBunch_ + 1);
171 
172  const float rollRadius(fixedRollRadius_ ? top_->radius() :
173  top_->radius() + CLHEP::RandFlat::shoot(engine, -1.*top_->stripLength()/2., top_->stripLength()/2.));
174 
175  //calculate noise from model
176  double averageNeutralNoiseRatePerRoll = 0.;
177  double averageNoiseElectronRatePerRoll = 0.;
178  double averageNoiseRatePerRoll = 0.;
179 
180  if(gemId.station() == 1)
181  {
182 //simulate neutral background for GE1/1
183  averageNeutralNoiseRatePerRoll = constNeuGE11_ * TMath::Exp(slopeNeuGE11_*rollRadius);
184 
185 //simulate eletron background for GE1/1
186 //the product is faster than Power or pow:
188  averageNoiseElectronRatePerRoll = GE11ElecBkgParams_[0]
189  + GE11ElecBkgParams_[1]*rollRadius
190  + GE11ElecBkgParams_[2]*rollRadius*rollRadius
191  + GE11ElecBkgParams_[3]*rollRadius*rollRadius*rollRadius;
192 
193  averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
194  }
195 
196  if(gemId.station() == 2 || gemId.station() == 3)
197  {
198 //simulate neutral background for GE2/1
199  averageNeutralNoiseRatePerRoll = GE21NeuBkgParams_[0]
200  + GE21NeuBkgParams_[1]*rollRadius
201  + GE21NeuBkgParams_[2]*rollRadius*rollRadius
202  + GE21NeuBkgParams_[3]*rollRadius*rollRadius*rollRadius
203  + GE21NeuBkgParams_[4]*rollRadius*rollRadius*rollRadius*rollRadius
204  + GE21NeuBkgParams_[5]*rollRadius*rollRadius*rollRadius*rollRadius*rollRadius;
205 
206 
207 //simulate eletron background for GE2/1
209  averageNoiseElectronRatePerRoll = GE21ElecBkgParams_[0]
210  + GE21ElecBkgParams_[1]*rollRadius
211  + GE21ElecBkgParams_[2]*rollRadius*rollRadius
212  + GE21ElecBkgParams_[3]*rollRadius*rollRadius*rollRadius
213  + GE21ElecBkgParams_[4]*rollRadius*rollRadius*rollRadius*rollRadius
214  + GE21ElecBkgParams_[5]*rollRadius*rollRadius*rollRadius*rollRadius*rollRadius
215  + GE21ElecBkgParams_[6]*rollRadius*rollRadius*rollRadius*rollRadius*rollRadius*rollRadius;
216 
217  averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
218  }
219 
220  //simulate intrinsic noise
222  {
223  const double aveIntrinsicNoisePerStrip(averageNoiseRate_ * nBxing * bxwidth_ * trStripArea * 1.0e-9);
224  for(int j = 0; j < nstrips; ++j)
225  {
226  CLHEP::RandPoissonQ randPoissonQ(*engine, aveIntrinsicNoisePerStrip);
227  const int n_intrHits(randPoissonQ.fire());
228 
229  for (int k = 0; k < n_intrHits; k++ )
230  {
231  const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
232  std::pair<int, int> digi(k+1,time_hit);
233  strips_.insert(digi);
234  }
235  }
236  }//end simulate intrinsic noise
237 
238  //simulate bkg contribution
239  const double averageNoise(averageNoiseRatePerRoll * nBxing * bxwidth_ * trArea * 1.0e-9 * scaleLumi_);
240  CLHEP::RandPoissonQ randPoissonQ(*engine, averageNoise);
241  const int n_hits(randPoissonQ.fire());
242 
243  for (int i = 0; i < n_hits; ++i)
244  {
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_);
247 
248  if (doNoiseCLS_)
249  {
250  std::vector<std::pair<int, int> > cluster_;
251  cluster_.clear();
252  cluster_.push_back(std::pair<int, int>(centralStrip, time_hit));
253 
254  int clusterSize = 0;
255  double randForCls = CLHEP::RandFlat::shoot(engine);
256 
257  if(randForCls <= clsParametrization_[0] && randForCls >= 0.)
258  clusterSize = 1;
259  else if(randForCls <= clsParametrization_[1] && randForCls > clsParametrization_[0])
260  clusterSize = 2;
261  else if(randForCls <= clsParametrization_[2] && randForCls > clsParametrization_[1])
262  clusterSize = 3;
263  else if(randForCls <= clsParametrization_[3] && randForCls > clsParametrization_[2])
264  clusterSize = 4;
265  else if(randForCls <= clsParametrization_[4] && randForCls > clsParametrization_[3])
266  clusterSize = 5;
267  else if(randForCls <= clsParametrization_[5] && randForCls > clsParametrization_[4])
268  clusterSize = 6;
269  else if(randForCls <= clsParametrization_[6] && randForCls > clsParametrization_[5])
270  clusterSize = 7;
271  else if(randForCls <= clsParametrization_[7] && randForCls > clsParametrization_[6])
272  clusterSize = 8;
273  else if(randForCls <= clsParametrization_[8] && randForCls > clsParametrization_[7])
274  clusterSize = 9;
275 
276  //odd cls
277  if (clusterSize % 2 != 0)
278  {
279  int clsR = (clusterSize - 1) / 2;
280  for (int i = 1; i <= clsR; ++i)
281  {
282  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - i > 0))
283  cluster_.push_back(std::pair<int, int>(centralStrip - i, time_hit));
284  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + i <= nstrips))
285  cluster_.push_back(std::pair<int, int>(centralStrip + i, time_hit));
286  }
287  }
288  //even cls
289  if (clusterSize % 2 == 0)
290  {
291  int clsR = (clusterSize - 2) / 2;
292  if(CLHEP::RandFlat::shoot(engine) < 0.5)
293  {
294  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
295  cluster_.push_back(std::pair<int, int>(centralStrip - 1, time_hit));
296  for (int i = 1; i <= clsR; ++i)
297  {
298  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 - i > 0))
299  cluster_.push_back(std::pair<int, int>(centralStrip - 1 - i, time_hit));
300  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + i <= nstrips))
301  cluster_.push_back(std::pair<int, int>(centralStrip + i, time_hit));
302  }
303  }
304 
305  else
306  {
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)
310  {
311  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 + i <= nstrips))
312  cluster_.push_back(std::pair<int, int>(centralStrip + 1 + i, time_hit));
313  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - i < 0))
314  cluster_.push_back(std::pair<int, int>(centralStrip - i, time_hit));
315  }
316  }
317  }
318  for(auto & digi : cluster_)
319  {
320  strips_.insert(digi);
321  }
322  }//end doNoiseCLS_
323  else
324  {
325  std::pair<int, int> digi(centralStrip, time_hit);
326  strips_.insert(digi);
327  }
328  }
329  return;
330 }
331 
332 std::vector<std::pair<int, int> > GEMSimpleModel::simulateClustering(const GEMEtaPartition* roll,
333  const PSimHit* simHit, const int bx, CLHEP::HepRandomEngine* engine)
334 {
335  // const Topology& topology(roll->specs()->topology());
336  const StripTopology& topology = roll->specificTopology();
337  // const LocalPoint& entry(simHit->entryPoint());
338  const LocalPoint& hit_position(simHit->localPosition());
339  const int nstrips(roll->nstrips());
340 
341  int centralStrip = 0;
342  if (!(topology.channel(hit_position) + 1 > nstrips))
343  centralStrip = topology.channel(hit_position) + 1;
344  else
345  centralStrip = topology.channel(hit_position);
346 
347  GlobalPoint pointSimHit = roll->toGlobal(hit_position);
348  GlobalPoint pointDigiHit = roll->toGlobal(roll->centreOfStrip(centralStrip));
349  double deltaphi = pointSimHit.phi() - pointDigiHit.phi();
350 
351  // Add central digi to cluster vector
352  std::vector<std::pair<int, int> > cluster_;
353  cluster_.clear();
354  cluster_.push_back(std::pair<int, int>(centralStrip, bx));
355 
356  // get the cluster size
357  int clusterSize = 0;
358  double randForCls = CLHEP::RandFlat::shoot(engine);
359 
360  if(randForCls <= clsParametrization_[0] && randForCls >= 0.)
361  clusterSize = 1;
362  else if(randForCls <= clsParametrization_[1] && randForCls > clsParametrization_[0])
363  clusterSize = 2;
364  else if(randForCls <= clsParametrization_[2] && randForCls > clsParametrization_[1])
365  clusterSize = 3;
366  else if(randForCls <= clsParametrization_[3] && randForCls > clsParametrization_[2])
367  clusterSize = 4;
368  else if(randForCls <= clsParametrization_[4] && randForCls > clsParametrization_[3])
369  clusterSize = 5;
370  else if(randForCls <= clsParametrization_[5] && randForCls > clsParametrization_[4])
371  clusterSize = 6;
372  else if(randForCls <= clsParametrization_[6] && randForCls > clsParametrization_[5])
373  clusterSize = 7;
374  else if(randForCls <= clsParametrization_[7] && randForCls > clsParametrization_[6])
375  clusterSize = 8;
376  else if(randForCls <= clsParametrization_[8] && randForCls > clsParametrization_[7])
377  clusterSize = 9;
378 
379  if (abs(simHit->particleType()) != 13 && fabs(simHit->pabs()) < minPabsNoiseCLS_)
380  return cluster_;
381 
382  //odd cls
383  if (clusterSize % 2 != 0)
384  {
385  int clsR = (clusterSize - 1) / 2;
386  for (int i = 1; i <= clsR; ++i)
387  {
388  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - i > 0))
389  cluster_.push_back(std::pair<int, int>(centralStrip - i, bx));
390  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + i <= nstrips))
391  cluster_.push_back(std::pair<int, int>(centralStrip + i, bx));
392  }
393  }
394  //even cls
395  if (clusterSize % 2 == 0)
396  {
397  int clsR = (clusterSize - 2) / 2;
398  if (deltaphi <= 0)
399  {
400  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
401  cluster_.push_back(std::pair<int, int>(centralStrip - 1, bx));
402  for (int i = 1; i <= clsR; ++i)
403  {
404  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 - i > 0))
405  cluster_.push_back(std::pair<int, int>(centralStrip - 1 - i, bx));
406  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + i <= nstrips))
407  cluster_.push_back(std::pair<int, int>(centralStrip + i, bx));
408  }
409  }
410  else
411  {
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)
415  {
416  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 + i <= nstrips))
417  cluster_.push_back(std::pair<int, int>(centralStrip + 1 + i, bx));
418  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - i < 0))
419  cluster_.push_back(std::pair<int, int>(centralStrip - i, bx));
420  }
421  }
422  }
423  return cluster_;
424 
425 }
426 
427 
double minPabsNoiseCLS_
std::vector< double > clsParametrization_
int i
Definition: DBlmapReader.cc:9
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.
Definition: GeomDet.h:52
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
double averageNoiseRate_
T y() const
Definition: PV3DBase.h:63
int getSimHitBx(const PSimHit *, CLHEP::HepRandomEngine *)
const Topology & topology() const
const GEMEtaPartition * etaPartition(GEMDetId id) const
Return a GEMEtaPartition given its id.
Definition: GEMGeometry.cc:103
GEMDetId id() const
LocalPoint centreOfStrip(int strip) const
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
int nstrips() const
number of readout strips in partition
float timeOfFlight() const
Definition: PSimHit.h:69
bool simulateIntrinsicNoise_
Local3DPoint localPosition() const
Definition: PSimHit.h:44
T sqrt(T t)
Definition: SSEVec.h:48
double averageEfficiency_
T z() const
Definition: PV3DBase.h:64
std::vector< double > GE11ElecBkgParams_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
std::vector< std::pair< int, int > > simulateClustering(const GEMEtaPartition *, const PSimHit *, const int, CLHEP::HepRandomEngine *) override
virtual int channel(const LocalPoint &p) const =0
const StripTopology & specificTopology() const
GEMSimpleModel(const edm::ParameterSet &)
double constNeuGE11_
float pabs() const
fast and more accurate access to momentumAtEntry().mag()
Definition: PSimHit.h:63
Container::value_type value_type
float pitch() const
bool simulateElectronBkg_
#define debug
Definition: HDRShower.cc:19
tuple simHits
Definition: trackerHits.py:16
void simulateNoise(const GEMEtaPartition *, CLHEP::HepRandomEngine *) override
const GEMGeometry * geometry_
Definition: GEMDigiModel.h:60
double averageShapingTime_
int particleType() const
Definition: PSimHit.h:85
std::set< std::pair< int, int > > strips_
Definition: GEMDigiModel.h:62
std::vector< double > GE21ElecBkgParams_
void clear()
Definition: DetSet.h:69
DetectorHitMap detectorHitMap_
Definition: GEMDigiModel.h:74
double timeResolution_
tuple cout
Definition: gather_cfg.py:121
std::vector< PSimHit > PSimHitContainer
StripDigiSimLinks stripDigiSimLinks_
Definition: GEMDigiModel.h:75
std::vector< double > GE21NeuBkgParams_
edm::DetSet< StripDigiSimLink > StripDigiSimLinks
Definition: GEMDigiModel.h:35
T x() const
Definition: PV3DBase.h:62
double slopeNeuGE11_
void simulateSignal(const GEMEtaPartition *, const edm::PSimHitContainer &, CLHEP::HepRandomEngine *) override
unsigned int detUnitId() const
Definition: PSimHit.h:93