1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
#include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h"
#include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithmFactory.h"
#include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitProduct.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
#include "DataFormats/TrackerRecHit2D/interface/FastSingleTrackerRecHit.h"
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include <string>
#include <iostream>
/* This plugin performs simple Gaussian smearing (GS) of the sim hit position
* per module. If y resolution < 0 (default) the module length/sqrt(12) is taken as uncertainty.
*/
class TrackingRecHitStripGSPlugin : public TrackingRecHitAlgorithm {
private:
double _resolutionX;
double _resolutionX2;
double _resolutionY;
double _resolutionY2;
constexpr static double INV12 = 1.0 / 12.0;
public:
TrackingRecHitStripGSPlugin(const std::string& name,
const edm::ParameterSet& config,
edm::ConsumesCollector& consumesCollector)
: TrackingRecHitAlgorithm(name, config, consumesCollector),
_resolutionX(0.001),
_resolutionX2(_resolutionX * _resolutionX),
_resolutionY(-1),
_resolutionY2(_resolutionY * _resolutionY)
{
if (config.exists("resolutionX")) {
_resolutionX = config.getParameter<double>("resolutionX");
_resolutionX2 = _resolutionX * _resolutionX;
}
if (config.exists("resolutionY")) {
_resolutionY = config.getParameter<double>("resolutionY");
_resolutionY2 = _resolutionY * _resolutionY;
}
}
TrackingRecHitProductPtr process(TrackingRecHitProductPtr product) const override {
for (const std::pair<unsigned int, const PSimHit*>& simHitIdPair : product->getSimHitIdPairs()) {
const PSimHit* simHit = simHitIdPair.second;
const Local3DPoint& simHitPosition = simHit->localPosition();
const GeomDet* geomDet = this->getTrackerGeometry().idToDetUnit(product->getDetId());
const Plane& plane = geomDet->surface();
const Bounds& bounds = plane.bounds();
const double boundY = bounds.length();
Local3DPoint recHitPosition;
unsigned int retry = 0;
do {
recHitPosition = Local3DPoint(simHitPosition.x() + this->getRandomEngine().gaussShoot(0.0, _resolutionX),
0.0,
0.0 //fix y & z coordinates to center of module
);
// If we tried to generate thePosition, and it's out of the bounds
// for 10 times, then take and return the simHit's location.
++retry;
if (retry > 10) {
recHitPosition = Local3DPoint(simHitPosition.x(), 0.0, 0.0 //fix y & z coordinates to center of module
);
break;
}
} while (not bounds.inside(recHitPosition));
LocalError error(
//xx (variance)
_resolutionX2,
//xy (covariance)
0.0,
//take here the provided y resolution or (lenght/sqrt(12))^2
_resolutionY < 0.0 ? boundY * boundY * INV12 : _resolutionY2);
FastSingleTrackerRecHit recHit(recHitPosition, error, *geomDet, fastTrackerRecHitType::siStrip1D);
product->addRecHit(recHit, {simHitIdPair});
}
return product;
}
};
DEFINE_EDM_PLUGIN(TrackingRecHitAlgorithmFactory, TrackingRecHitStripGSPlugin, "TrackingRecHitStripGSPlugin");
|