|
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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
|
#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <iostream>
#include <cmath>
#include <algorithm>
//#define EDM_ML_DEBUG
TrapezoidalStripTopology::TrapezoidalStripTopology(int ns, float p, float l, float r0)
: theNumberOfStrips(ns), thePitch(p), theDistToBeam(r0), theDetHeight(l) {
theOffset = -theNumberOfStrips / 2. * thePitch;
theYAxOr = 1;
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("CommonTopologies") << "Constructing TrapezoidalStripTopology with nstrips = " << ns
<< " pitch = " << p << " length = " << l << " r0 =" << r0;
#endif
}
TrapezoidalStripTopology::TrapezoidalStripTopology(int ns, float p, float l, float r0, int yAx)
: theNumberOfStrips(ns), thePitch(p), theDistToBeam(r0), theDetHeight(l), theYAxOr(yAx) {
theOffset = -theNumberOfStrips / 2. * thePitch;
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("CommonTopologies") << "Constructing TrapezoidalStripTopology with nstrips = " << ns
<< " pitch = " << p << " length = " << l << " r0 =" << r0
<< " yAxOrientation =" << yAx;
#endif
}
LocalPoint TrapezoidalStripTopology::localPosition(float strip) const {
return LocalPoint(strip * thePitch + theOffset, 0.0);
}
LocalPoint TrapezoidalStripTopology::localPosition(const MeasurementPoint& mp) const {
float y = mp.y() * theDetHeight;
float x = (mp.x() * thePitch + theOffset) * (theYAxOr * y + theDistToBeam) / theDistToBeam;
return LocalPoint(x, y);
}
LocalError TrapezoidalStripTopology::localError(float strip, float stripErr2) const {
float lt, lc2, ls2, lslc;
float localL2, localP2;
float sl2, sp2;
// angle from strip to local frame (see CMS TN / 95-170)
lt = -(strip * thePitch + theOffset) * theYAxOr / theDistToBeam;
lc2 = 1.f / (1. + lt * lt);
lslc = lt * lc2;
ls2 = 1.f - lc2;
localL2 = theDetHeight * theDetHeight / lc2;
localP2 = thePitch * thePitch * lc2;
sl2 = localL2 / 12.;
sp2 = stripErr2 * localP2;
return LocalError(lc2 * sp2 + ls2 * sl2, lslc * (sp2 - sl2), ls2 * sp2 + lc2 * sl2);
}
LocalError TrapezoidalStripTopology::localError(const MeasurementPoint& mp, const MeasurementError& merr) const {
float lt, lc2, ls2, lslc;
float localL, localP;
float sl2, sp2, spl;
// angle from strip to local frame (see CMS TN / 95-170)
lt = -(mp.x() * thePitch + theOffset) * theYAxOr / theDistToBeam;
lc2 = 1. / (1. + lt * lt);
lslc = lt * lc2;
ls2 = 1.f - lc2;
localL = theDetHeight / std::sqrt(lc2);
localP = localPitch(localPosition(mp));
sp2 = merr.uu() * localP * localP;
sl2 = merr.vv() * localL * localL;
spl = merr.uv() * localP * localL;
return LocalError(lc2 * sp2 + ls2 * sl2 - 2 * lslc * spl,
lslc * (sp2 - sl2) + (lc2 - ls2) * spl,
ls2 * sp2 + lc2 * sl2 + 2 * lslc * spl);
}
float TrapezoidalStripTopology::strip(const LocalPoint& lp) const {
float aStrip = ((lp.x() * theDistToBeam / (theYAxOr * lp.y() + theDistToBeam)) - theOffset) / thePitch;
if (aStrip < 0)
aStrip = 0;
else if (aStrip > theNumberOfStrips)
aStrip = theNumberOfStrips;
return aStrip;
}
MeasurementPoint TrapezoidalStripTopology::measurementPosition(const LocalPoint& lp) const {
return MeasurementPoint(((lp.x() * theDistToBeam / (theYAxOr * lp.y() + theDistToBeam)) - theOffset) / thePitch,
lp.y() / theDetHeight);
}
MeasurementError TrapezoidalStripTopology::measurementError(const LocalPoint& lp, const LocalError& lerr) const {
float lt, lc2, ls2, lslc;
float localL, localP;
float sl2, sp2, spl;
lt = -lp.x() / (theYAxOr * lp.y() + theDistToBeam) * theYAxOr;
lc2 = 1. / (1. + lt * lt);
lslc = lt * lc2;
ls2 = 1. - lc2;
localL = theDetHeight / std::sqrt(lc2);
localP = localPitch(lp);
sp2 = lc2 * lerr.xx() + ls2 * lerr.yy() + 2 * lslc * lerr.xy();
sl2 = ls2 * lerr.xx() + lc2 * lerr.yy() - 2 * lslc * lerr.xy();
spl = lslc * (lerr.yy() - lerr.xx()) + (lc2 - ls2) * lerr.xy();
return MeasurementError(sp2 / (localP * localP), spl / (localP * localL), sl2 / (localL * localL));
}
int TrapezoidalStripTopology::channel(const LocalPoint& lp) const {
return std::min(int(strip(lp)), theNumberOfStrips - 1);
}
float TrapezoidalStripTopology::pitch() const { return thePitch; }
float TrapezoidalStripTopology::localPitch(const LocalPoint& lp) const {
float x = lp.x();
float y = theYAxOr * lp.y() + theDistToBeam;
return thePitch * y / (theDistToBeam * std::sqrt(1.f + x * x / (y * y)));
}
float TrapezoidalStripTopology::stripAngle(float strip) const {
return std::atan(-(strip * thePitch + theOffset) * theYAxOr / theDistToBeam);
}
int TrapezoidalStripTopology::nstrips() const { return theNumberOfStrips; }
float TrapezoidalStripTopology::shiftOffset(float pitch_fraction) {
theOffset += thePitch * pitch_fraction;
return theOffset;
}
float TrapezoidalStripTopology::localStripLength(const LocalPoint& lp) const {
float ltan = -lp.x() / (theYAxOr * lp.y() + theDistToBeam) * theYAxOr;
float localL = theDetHeight * std::sqrt(1.f + ltan * ltan);
// float lcos2 = 1.f/(1.f+ltan*ltan);
// float localL = theDetHeight / std::sqrt(lcos2);
return localL;
}
|