|
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
|
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerGeometryBuilder/interface/ProxyStripTopology.h"
#include "Geometry/CommonTopologies/interface/TkRadialStripTopology.h"
#include "TFile.h"
#include "TProfile.h"
class ValidateRadial : public edm::one::EDAnalyzer<> {
public:
ValidateRadial(const edm::ParameterSet&);
~ValidateRadial() override;
void beginJob() override {}
void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
void endJob() override {}
private:
std::vector<const TkRadialStripTopology*> get_list_of_radial_topologies(const edm::Event&, const edm::EventSetup&);
void test_topology(const TkRadialStripTopology*, unsigned);
bool pass_frame_change_test(const TkRadialStripTopology* t, float strip, float stripErr2, bool);
bool EQUAL(const double a, const double b) { return fabs(a - b) < epsilon_; }
const double epsilon_;
TFile* file_;
const bool printOut_;
const bool posOnly_;
const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tokTrackerGeometry_;
mutable float maxerrU = 0.;
mutable float maxerrUV = 0.;
};
ValidateRadial::ValidateRadial(const edm::ParameterSet& cfg)
: epsilon_(cfg.getParameter<double>("Epsilon")),
file_(new TFile(cfg.getParameter<std::string>("FileName").c_str(), "RECREATE")),
printOut_(cfg.getParameter<bool>("PrintOut")),
posOnly_(cfg.getParameter<bool>("PosOnly")),
tokTrackerGeometry_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()) {
std::cout << "I'm ALIVE" << std::endl;
}
void ValidateRadial::analyze(const edm::Event& e, const edm::EventSetup& es) {
std::vector<const TkRadialStripTopology*> topologies = get_list_of_radial_topologies(e, es);
for (unsigned i = 0; i < topologies.size(); i++) {
test_topology(topologies[i], i);
}
file_->Close();
}
std::vector<const TkRadialStripTopology*> ValidateRadial::get_list_of_radial_topologies(const edm::Event& e,
const edm::EventSetup& es) {
std::vector<const TkRadialStripTopology*> topos;
auto theTrackerGeometry = es.getHandle(tokTrackerGeometry_);
const uint32_t radial_detids[] = {402666125, //TID r1
402668833, //TID r2
402673476, //TID r3
470066725, //TEC r1
470390853, //TEC r2
470114664, //TEC r3
470131344, //TEC r4
470079661, //TEC r5
470049476, //TEC r6
470045428}; //TEC r7
for (unsigned int radial_detid : radial_detids) {
auto g = dynamic_cast<const StripGeomDetUnit*>(theTrackerGeometry->idToDet(radial_detid));
if (!g)
std::cout << "no geom for " << radial_detid << std::endl;
auto const topol = &g->specificTopology();
const TkRadialStripTopology* rt = nullptr;
auto const proxyT = dynamic_cast<const ProxyStripTopology*>(topol);
if (proxyT)
rt = dynamic_cast<const TkRadialStripTopology*>(&(proxyT->specificTopology()));
else
rt = dynamic_cast<const TkRadialStripTopology*>(topol);
if (!rt)
std::cout << "no radial topology for " << radial_detid << std::endl;
else
topos.emplace_back(rt);
}
return topos;
}
#include "Geometry/CommonTopologies/interface/CSCRadialStripTopology.h"
void compare(const TkRadialStripTopology& t,
const CSCRadialStripTopology& ot,
const float strip,
const float stripErr2) {
const LocalPoint lp = t.localPosition(strip);
const LocalError le = t.localError(strip, stripErr2);
const MeasurementPoint mp = t.measurementPosition(lp);
const MeasurementError me = t.measurementError(lp, le);
const float newstrip = t.strip(lp);
const LocalPoint olp = ot.localPosition(strip);
const LocalError ole = ot.localError(strip, stripErr2);
const MeasurementPoint omp = ot.measurementPosition(lp);
const MeasurementError ome = ot.measurementError(lp, le);
const float onewstrip = ot.strip(lp);
if (fabs(lp.x() - olp.x()) > 0.001)
std::cout << "FAILED " << lp.x() << " " << olp.x() << std::endl;
if (fabs(lp.y() - olp.y()) > 0.001)
std::cout << "FAILED " << lp.y() << " " << olp.y() << std::endl;
if (fabs(newstrip - onewstrip) > 0.001)
std::cout << "FAILED " << newstrip << " " << onewstrip << std::endl;
if (fabs(le.xx() - ole.xx()) > 0.001)
std::cout << "FAILED " << le.xx() << " " << ole.xx() << std::endl;
if (fabs(le.yy() - ole.yy()) > 0.001)
std::cout << "FAILED " << le.yy() << " " << ole.xy() << std::endl;
if (fabs(le.xy() - ole.xy()) > 0.001)
std::cout << "FAILED " << le.xy() << " " << ole.yy() << std::endl;
if (fabs(mp.x() - omp.x()) > 0.001)
std::cout << "FAILED " << mp.x() << " " << omp.x() << std::endl;
if (fabs(me.uu() - ome.uu()) > 0.001)
std::cout << "FAILED " << me.uu() << " " << ome.uu() << std::endl;
}
void ValidateRadial::test_topology(const TkRadialStripTopology* t, unsigned i) {
CSCRadialStripTopology oldt(t->nstrips(),
t->angularWidth(),
t->detHeight(),
t->centreToIntersection(),
t->yAxisOrientation(),
t->yCentreOfStripPlane());
std::cout << "\nTK\n" << *t << std::endl;
std::cout << "\nCSC\n" << oldt << std::endl;
TProfile prof(("se2limit1" + std::to_string(i)).c_str(),
"Precision Limit of recoverable strip error (1st order);strip;strip error",
t->nstrips() / 8,
0,
t->nstrips());
TProfile prof2(("se2limit2" + std::to_string(i)).c_str(),
"Precision Limit of recoverable strip error (2nd order);strip;strip error",
t->nstrips() / 8,
0,
t->nstrips());
for (float strip = 0; strip < t->nstrips(); strip += 0.5) {
for (float stripErr2 = 0.03; stripErr2 > 1e-10; stripErr2 /= 1.1) {
compare(*t, oldt, strip, stripErr2);
if (!pass_frame_change_test(t, strip, stripErr2, false)) {
prof.Fill(strip, sqrt(stripErr2));
break;
}
}
for (float stripErr2 = 0.03; stripErr2 > 1e-10; stripErr2 /= 1.1)
if (!pass_frame_change_test(t, strip, stripErr2, true)) {
prof2.Fill(strip, sqrt(stripErr2));
break;
}
}
prof.Write();
prof2.Write();
}
ValidateRadial::~ValidateRadial() {
std::cout << "ValidateRadial max UU, max UV " << maxerrU << " " << maxerrUV << std::endl;
}
bool ValidateRadial::pass_frame_change_test(const TkRadialStripTopology* t,
const float strip,
const float stripErr2,
const bool secondOrder) {
const LocalPoint lp = t->localPosition(strip);
const LocalError le = t->localError(strip, stripErr2);
const MeasurementPoint mp = t->measurementPosition(lp);
const MeasurementError me = t->measurementError(lp, le);
const float newstrip = t->strip(lp);
const LocalPoint newlp = t->localPosition(mp);
const LocalError newle = t->localError(mp, me);
const MeasurementPoint newmp = t->measurementPosition(newlp);
const MeasurementError newme = t->measurementError(newlp, newle);
maxerrU = std::max(maxerrU, std::abs(me.uu() - stripErr2) / stripErr2);
maxerrUV = std::max(maxerrUV, std::abs(me.uv()));
const bool pass1p = fabs(strip - newstrip) < 0.001 && fabs(strip - mp.x()) < 0.001 && fabs(mp.x() - newstrip) < 0.001;
const bool pass1e =
(fabs(me.uu() - stripErr2) / stripErr2 < epsilon_) && (fabs(me.uv()) < epsilon_) && me.uu() > 0 && me.vv() > 0;
const bool pass2p = fabs(strip - newmp.x()) < 0.001;
const bool pass2e = (fabs(newme.uu() - stripErr2) / stripErr2 < epsilon_) && (fabs(newme.uv()) < epsilon_) &&
newme.uu() > 0 && newme.vv() > 0;
const bool passp = (secondOrder ? pass2p : pass1p);
const bool passe = (secondOrder ? pass2e : pass1e);
if (printOut_ && ((!passp) || ((!posOnly_) && (!passe))))
std::cout << "FAILED "
<< "(" << strip << ", " << newstrip << ", " << mp.x() << ")\t"
<< "(" << stripErr2 << ", " << me.uu() << ")\t\t" << (me.uv()) << std::endl;
return passp & passe;
}
DEFINE_FWK_MODULE(ValidateRadial);
|