CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ME0SimpleModel.cc
Go to the documentation of this file.
1 #include "CLHEP/Units/defs.h"
2 #include "CLHEP/Units/SystemOfUnits.h"
3 #include "CLHEP/Units/GlobalPhysicalConstants.h"
9 #include "CLHEP/Random/RandFlat.h"
10 #include "CLHEP/Random/RandPoissonQ.h"
11 #include "CLHEP/Random/RandGaussQ.h"
12 #include <cmath>
13 #include <utility>
14 #include <map>
15 
16 
18 ME0DigiModel(config)
19 , averageEfficiency_(config.getParameter<double> ("averageEfficiency"))
20 , averageShapingTime_(config.getParameter<double> ("averageShapingTime"))
21 , timeResolution_(config.getParameter<double> ("timeResolution"))
22 , timeJitter_(config.getParameter<double> ("timeJitter"))
23 , averageNoiseRate_(config.getParameter<double> ("averageNoiseRate"))
24 , signalPropagationSpeed_(config.getParameter<double> ("signalPropagationSpeed"))
25 , bxwidth_(config.getParameter<int> ("bxwidth"))
26 , minBunch_(config.getParameter<int> ("minBunch"))
27 , maxBunch_(config.getParameter<int> ("maxBunch"))
28 , digitizeOnlyMuons_(config.getParameter<bool> ("digitizeOnlyMuons"))
29 , doBkgNoise_(config.getParameter<bool> ("doBkgNoise"))
30 , doNoiseCLS_(config.getParameter<bool> ("doNoiseCLS"))
31 , fixedRollRadius_(config.getParameter<bool> ("fixedRollRadius"))
32 , simulateElectronBkg_(config.getParameter<bool> ("simulateElectronBkg"))
33 , instLumi_(config.getParameter<double> ("instLumi"))
34 , rateFact_(config.getParameter<double> ("rateFact"))
35 , referenceInstLumi_(config.getParameter<double> ("referenceInstLumi"))
36 , ME0ElecBkgParam0_(config.getParameter<double> ("ME0ElecBkgParam0"))
37 , ME0ElecBkgParam1_(config.getParameter<double> ("ME0ElecBkgParam1"))
38 , ME0ElecBkgParam2_(config.getParameter<double> ("ME0ElecBkgParam2"))
39 , ME0ElecBkgParam3_(config.getParameter<double> ("ME0ElecBkgParam3"))
40 , ME0NeuBkgParam0_(config.getParameter<double> ("ME0NeuBkgParam0"))
41 , ME0NeuBkgParam1_(config.getParameter<double> ("ME0NeuBkgParam1"))
42 , ME0NeuBkgParam2_(config.getParameter<double> ("ME0NeuBkgParam2"))
43 , ME0NeuBkgParam3_(config.getParameter<double> ("ME0NeuBkgParam3"))
44 {
45 }
46 
48 {
49 }
50 
52 {
53  return;
54 }
55 
56 void ME0SimpleModel::simulateSignal(const ME0EtaPartition* roll, const edm::PSimHitContainer& simHits, CLHEP::HepRandomEngine* engine)
57 {
59  detectorHitMap_.clear();
61 
64  bool digiMuon = false;
65  bool digiElec = false;
66  for (const auto& hit : simHits)
67  {
68  if (std::abs(hit.particleType()) != 13 && digitizeOnlyMuons_)
69  continue;
70  double elecEff = 0.;
71  double partMom = hit.pabs();
72  double checkMuonEff = CLHEP::RandFlat::shoot(engine, 0., 1.);
73  double checkElecEff = CLHEP::RandFlat::shoot(engine, 0., 1.);
74  if (std::abs(hit.particleType()) == 13 && checkMuonEff < averageEfficiency_)
75  digiMuon = true;
76  if (std::abs(hit.particleType()) != 13) //consider all non muon particles with me0 efficiency to electrons
77  {
78  if (partMom <= 1.95e-03)
79  elecEff = 1.7e-05 * std::exp(2.1 * partMom * 1000.);
80  if (partMom > 1.95e-03 && partMom < 10.e-03)
81  elecEff = 1.34 * log(7.96e-01 * partMom * 1000. - 5.75e-01)
82  / (1.34 + log(7.96e-01 * partMom * 1000. - 5.75e-01));
83  if (partMom > 10.e-03)
84  elecEff = 1.;
85  if (checkElecEff < elecEff)
86  digiElec = true;
87  }
88  if (!(digiMuon || digiElec))
89  continue;
90  const int bx(getSimHitBx(&hit, engine));
91  const std::vector<std::pair<int, int> >& cluster(simulateClustering(roll, &hit, bx, engine));
92  for (const auto & digi : cluster)
93  {
94  detectorHitMap_.emplace(digi,&hit);
95  strips_.emplace(digi);
96  }
97  }
98 }
99 
100 int ME0SimpleModel::getSimHitBx(const PSimHit* simhit, CLHEP::HepRandomEngine* engine)
101 {
102  int bx = -999;
103  const LocalPoint& simHitPos(simhit->localPosition());
104  const float tof(simhit->timeOfFlight());
105  // random Gaussian time correction due to electronics jitter
106  float randomJitterTime = CLHEP::RandGaussQ::shoot(engine, 0., timeJitter_);
107  const ME0DetId& id(simhit->detUnitId());
108  const ME0EtaPartition* roll(geometry_->etaPartition(id));
109  if (!roll)
110  {
111  throw cms::Exception("Geometry")<< "ME0SimpleModel::getSimHitBx() - ME0 simhit id does not match any ME0 roll id: " << id << "\n";
112  return 999;
113  }
114  if (roll->id().region() == 0)
115  {
116  throw cms::Exception("Geometry") << "ME0SimpleModel::getSimHitBx() - this ME0 id is from barrel, which cannot happen: " << roll->id() << "\n";
117  }
118  const int nstrips = roll->nstrips();
119  float middleStrip = nstrips/2.;
120  const LocalPoint& middleOfRoll = roll->centreOfStrip(middleStrip);
121  const GlobalPoint& globMiddleRol = roll->toGlobal(middleOfRoll);
122  double muRadius = sqrt(globMiddleRol.x()*globMiddleRol.x() + globMiddleRol.y()*globMiddleRol.y() +globMiddleRol.z()*globMiddleRol.z());
123  double timeCalibrationOffset_ = (muRadius*CLHEP::ns*CLHEP::cm)/(CLHEP::c_light); //[cm/ns]
124  const TrapezoidalStripTopology* top(dynamic_cast<const TrapezoidalStripTopology*> (&(roll->topology())));
125  const float halfStripLength(0.5 * top->stripLength());
126  const float distanceFromEdge(halfStripLength - simHitPos.y());
127 
128  // signal propagation speed in material in [cm/ns]
129  double signalPropagationSpeedTrue = signalPropagationSpeed_ * CLHEP::c_light/(CLHEP::ns*CLHEP::cm);
130  // average time for the signal to propagate from the SimHit to the top of a strip
131  const float averagePropagationTime(distanceFromEdge / signalPropagationSpeedTrue);
132  // random Gaussian time correction due to the finite timing resolution of the detector
133  float randomResolutionTime = CLHEP::RandGaussQ::shoot(engine, 0., timeResolution_);
134  const float simhitTime(tof + averageShapingTime_ + randomResolutionTime + averagePropagationTime + randomJitterTime);
135  float referenceTime = 0.;
136  referenceTime = timeCalibrationOffset_ + halfStripLength / signalPropagationSpeedTrue + averageShapingTime_;
137  const float timeDifference(simhitTime - referenceTime);
138  // assign the bunch crossing
139  bx = static_cast<int> (std::round((timeDifference) / bxwidth_));
140 
141  // check time
142  const bool debug(false);
143  if (debug)
144  {
145  LogDebug("ME0SimpleModel")<< "checktime " << "bx = " << bx << "\tdeltaT = " << timeDifference << "\tsimT = " << simhitTime
146  << "\trefT = " << referenceTime << "\ttof = " << tof << "\tavePropT = " << averagePropagationTime
147  << "\taveRefPropT = " << halfStripLength / signalPropagationSpeedTrue << std::endl;
148  }
149  return bx;
150 }
151 
152 void ME0SimpleModel::simulateNoise(const ME0EtaPartition* roll, CLHEP::HepRandomEngine* engine)
153 {
154  if (!doBkgNoise_)
155  return;
156  const ME0DetId me0Id(roll->id());
157  const int nstrips(roll->nstrips());
158  double trArea(0.0);
159  double trStripArea(0.0);
160 
161  if (me0Id.region() == 0)
162  {
163  throw cms::Exception("Geometry")
164  << "ME0Synchronizer::simulateNoise() - this ME0 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  const float rollRadius(fixedRollRadius_ ? top_->radius() :
172  top_->radius() + CLHEP::RandFlat::shoot(engine, -1.*top_->stripLength()/2., top_->stripLength()/2.));
173 
174  const float rSqrtR = rollRadius * sqrt(rollRadius);
175 
176  //calculate noise from model
177  double averageNeutralNoiseRatePerRoll = 0.;
178  double averageNoiseElectronRatePerRoll = 0.;
179  double averageNoiseRatePerRoll = 0.;
180  if (me0Id.station() != 1)
181  {
182  throw cms::Exception("Geometry") << "ME0SimpleModel::simulateNoise() - this ME0 id is from station 1, which cannot happen: " <<roll->id() << "\n";
183  }
184  else
185  {
186  averageNeutralNoiseRatePerRoll = ME0NeuBkgParam0_ * rollRadius* std::exp(ME0NeuBkgParam1_/rSqrtR) + ME0NeuBkgParam2_/rSqrtR + ME0NeuBkgParam3_/(sqrt(rollRadius));
187  //simulate electron background for ME0
189  averageNoiseElectronRatePerRoll = ME0ElecBkgParam0_ * rSqrtR* std::exp(ME0ElecBkgParam1_/rSqrtR) + ME0ElecBkgParam2_/rSqrtR + ME0ElecBkgParam3_/(sqrt(rollRadius));
190  averageNoiseRatePerRoll = averageNeutralNoiseRatePerRoll + averageNoiseElectronRatePerRoll;
191  averageNoiseRatePerRoll *= instLumi_*rateFact_*1.0/referenceInstLumi_;
192  }
193 
194  //simulate intrinsic noise
196  {
197  const double aveIntrinsicNoisePerStrip(averageNoiseRate_ * nBxing * bxwidth_ * trStripArea * 1.0e-9);
198  for(int j = 0; j < nstrips; ++j)
199  {
200  int randPoissonQ = CLHEP::RandPoissonQ::shoot(engine, aveIntrinsicNoisePerStrip);
201 
202  for (int k = 0; k < randPoissonQ; k++ )
203  {
204  const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
205  strips_.emplace(k+1,time_hit);
206  }
207  }
208  }//end simulate intrinsic noise
209 
210  //simulate bkg contribution
211  const double averageNoise(averageNoiseRatePerRoll * nBxing * bxwidth_ * trArea * 1.0e-9);
212  int randPoissonQ = CLHEP::RandPoissonQ::shoot(engine, averageNoise);
213 
214  for (int i = 0; i < randPoissonQ; ++i)
215  {
216  const int centralStrip(static_cast<int> (CLHEP::RandFlat::shoot(engine, 1, nstrips)));
217  const int time_hit(static_cast<int>(CLHEP::RandFlat::shoot(engine, nBxing)) + minBunch_);
218  if (doNoiseCLS_)
219  {
220  std::vector < std::pair<int, int> > cluster;
221  cluster.emplace_back(centralStrip, time_hit);
222  int clusterSize((CLHEP::RandFlat::shoot(engine)) <= 0.53 ? 1 : 2);
223  if (clusterSize == 2)
224  {
225  if(CLHEP::RandFlat::shoot(engine) < 0.5)
226  {
227  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
228  cluster.emplace_back(centralStrip - 1, time_hit);
229  }
230  else
231  {
232  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 <= nstrips))
233  cluster.emplace_back(centralStrip + 1, time_hit);
234  }
235  }
236  for (const auto& digi : cluster)
237  {
238  strips_.emplace(digi);
239  }
240  } //end doNoiseCLS_
241  else
242  {
243  strips_.emplace(centralStrip, time_hit);
244  }
245  }
246  return;
247 }
248 
249 std::vector<std::pair<int, int> > ME0SimpleModel::simulateClustering(const ME0EtaPartition* roll,
250  const PSimHit* simHit, const int bx, CLHEP::HepRandomEngine* engine)
251 {
252  const StripTopology& topology = roll->specificTopology();
253  const LocalPoint& hit_position(simHit->localPosition());
254  const int nstrips(roll->nstrips());
255  int centralStrip = 0;
256  if (!(topology.channel(hit_position) + 1 > nstrips))
257  centralStrip = topology.channel(hit_position) + 1;
258  else
259  centralStrip = topology.channel(hit_position);
260  const GlobalPoint& pointSimHit = roll->toGlobal(hit_position);
261  const GlobalPoint& pointDigiHit = roll->toGlobal(roll->centreOfStrip(centralStrip));
262  double deltaX = pointSimHit.x() - pointDigiHit.x();
263 
264  // Add central digi to cluster vector
265  std::vector < std::pair<int, int> > cluster;
266  cluster.clear();
267  cluster.emplace_back(centralStrip, bx);
268 
269  //simulate cross talk
270  int clusterSize((CLHEP::RandFlat::shoot(engine)) <= 0.53 ? 1 : 2);
271  if (clusterSize == 2)
272  {
273  if (deltaX <= 0)
274  {
275  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip - 1 > 0))
276  cluster.emplace_back(centralStrip - 1, bx);
277  }
278  else
279  {
280  if (CLHEP::RandFlat::shoot(engine) < averageEfficiency_ && (centralStrip + 1 <= nstrips))
281  cluster.emplace_back(centralStrip + 1, bx);
282  }
283  }
284  return cluster;
285 }
#define LogDebug(id)
const StripTopology & specificTopology() const
double averageShapingTime_
void setup() override
ME0SimpleModel(const edm::ParameterSet &)
edm::DetSet< ME0DigiSimLink > ME0DigiSimLinks
Definition: ME0DigiModel.h:37
CaloTopology const * topology(0)
double ME0NeuBkgParam2_
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:54
double averageEfficiency_
bool simulateElectronBkg_
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:50
T y() const
Definition: PV3DBase.h:63
StripDigiSimLinks stripDigiSimLinks_
Definition: ME0DigiModel.h:79
std::set< std::pair< int, int > > strips_
Definition: ME0DigiModel.h:65
LocalPoint centreOfStrip(int strip) const
Definition: config.py:1
const ME0EtaPartition * etaPartition(ME0DetId id) const
Return a etaPartition given its id.
Definition: ME0Geometry.cc:64
int getSimHitBx(const PSimHit *, CLHEP::HepRandomEngine *)
const ME0Geometry * geometry_
Definition: ME0DigiModel.h:63
float timeOfFlight() const
Definition: PSimHit.h:73
Local3DPoint localPosition() const
Definition: PSimHit.h:52
double averageNoiseRate_
T sqrt(T t)
Definition: SSEVec.h:18
ME0DetId id() const
T z() const
Definition: PV3DBase.h:64
float pitch() const
double ME0ElecBkgParam0_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double referenceInstLumi_
double ME0NeuBkgParam0_
virtual int channel(const LocalPoint &p) const =0
int nstrips() const
Return the chamber this roll belongs to.
int k[5][pyjets_maxn]
double ME0ElecBkgParam2_
ME0DigiSimLinks theME0DigiSimLinks_
Definition: ME0DigiModel.h:80
#define debug
Definition: HDRShower.cc:19
double ME0NeuBkgParam3_
double signalPropagationSpeed_
double ME0ElecBkgParam3_
void simulateNoise(const ME0EtaPartition *, CLHEP::HepRandomEngine *) override
std::vector< std::pair< int, int > > simulateClustering(const ME0EtaPartition *, const PSimHit *, const int, CLHEP::HepRandomEngine *) override
~ME0SimpleModel() override
bool simulateIntrinsicNoise_
DetectorHitMap detectorHitMap_
Definition: ME0DigiModel.h:78
void clear()
Definition: DetSet.h:71
edm::DetSet< StripDigiSimLink > StripDigiSimLinks
Definition: ME0DigiModel.h:36
std::vector< PSimHit > PSimHitContainer
void simulateSignal(const ME0EtaPartition *, const edm::PSimHitContainer &, CLHEP::HepRandomEngine *) override
double ME0ElecBkgParam1_
T x() const
Definition: PV3DBase.h:62
const Topology & topology() const override
double ME0NeuBkgParam1_
unsigned int detUnitId() const
Definition: PSimHit.h:97
double timeResolution_