20 #include "MeasurementCreator.h"
24 #include <PlanarMeasurement.h>
25 #include <ProlateSpacepointMeasurement.h>
26 #include <SpacepointMeasurement.h>
27 #include <WireMeasurement.h>
28 #include <WirePointMeasurement.h>
38 MeasurementCreator::MeasurementCreator() :
50 idealLRResolution_(true),
54 measurementCounter_(0),
61 std::vector<genfit::AbsMeasurement*> MeasurementCreator::create(eMeasurementType type,
double tracklength,
bool& outlier,
int& lr) {
65 std::vector<AbsMeasurement*> retVal;
69 trackModel_->getPosDir(tracklength, point, dir);
72 TVector3 planeNorm(dir);
73 planeNorm.SetTheta(thetaDetPlane_*TMath::Pi()/180);
74 planeNorm.SetPhi(planeNorm.Phi()+phiDetPlane_);
75 static const TVector3 z(0,0,1);
76 static const TVector3 x(1,0,0);
79 TVector3 currentWireDir(wireDir_);
86 if (useSkew_ && (
int)((
double)wireCounter_/(
double)nSuperLayer_)%2 == 1) {
87 TVector3 perp(wireDir_.Cross(dir));
88 if ((
int)((
double)wireCounter_/(
double)nSuperLayer_)%4 == 1){
89 currentWireDir.Rotate(skewAngle_*TMath::Pi()/180, wireDir_.Cross(perp));
91 else currentWireDir.Rotate(-skewAngle_*TMath::Pi()/180, wireDir_.Cross(perp));
93 currentWireDir.SetMag(1.);
97 wirePerp = dir.Cross(currentWireDir);
98 if (gRandom->Uniform(-1,1) >= 0) {
102 wirePerp.SetMag(gRandom->Uniform(minDrift_, maxDrift_));
105 if (outlierProb_ > gRandom->Uniform(1.)) {
107 if(debug_) std::cerr <<
"create outlier" << std::endl;
113 if (debug_) std::cerr <<
"create PixHit" << std::endl;
117 TVectorD hitCoords(2);
119 hitCoords(0) = gRandom->Uniform(-outlierRange_, outlierRange_);
120 hitCoords(1) = gRandom->Uniform(-outlierRange_, outlierRange_);
123 hitCoords(0) = gRandom->Gaus(0,resolution_);
124 hitCoords(1) = gRandom->Gaus(0,resolution_);
127 TMatrixDSym hitCov(2);
128 hitCov(0,0) = resolution_*resolution_;
129 hitCov(1,1) = resolution_*resolution_;
133 retVal.push_back(measurement);
138 if (debug_) std::cerr <<
"create SpacepointHit" << std::endl;
140 TVectorD hitCoords(3);
142 hitCoords(0) = gRandom->Uniform(point.X()-outlierRange_, point.X()+outlierRange_);
143 hitCoords(1) = gRandom->Uniform(point.Y()-outlierRange_, point.Y()+outlierRange_);
144 hitCoords(2) = gRandom->Uniform(point.Z()-outlierRange_, point.Z()+outlierRange_);
147 hitCoords(0) = gRandom->Gaus(point.X(),resolution_);
148 hitCoords(1) = gRandom->Gaus(point.Y(),resolution_);
149 hitCoords(2) = gRandom->Gaus(point.Z(),resolution_);
152 TMatrixDSym hitCov(3);
153 hitCov(0,0) = resolution_*resolution_;
154 hitCov(1,1) = resolution_*resolution_;
155 hitCov(2,2) = resolution_*resolution_;
158 retVal.push_back(measurement);
162 case ProlateSpacepoint: {
163 if (debug_) std::cerr <<
"create ProlateSpacepointHit" << std::endl;
165 TVectorD hitCoords(3);
167 hitCoords(0) = gRandom->Uniform(point.X()-outlierRange_, point.X()+outlierRange_);
168 hitCoords(1) = gRandom->Uniform(point.Y()-outlierRange_, point.Y()+outlierRange_);
169 hitCoords(2) = gRandom->Uniform(point.Z()-outlierRange_, point.Z()+outlierRange_);
172 hitCoords(0) = point.X();
173 hitCoords(1) = point.Y();
174 hitCoords(2) = point.Z();
177 TMatrixDSym hitCov(3);
178 hitCov(0,0) = resolution_*resolution_;
179 hitCov(1,1) = resolution_*resolution_;
180 hitCov(2,2) = resolutionWire_*resolutionWire_;
183 TVector3
xp = currentWireDir.Orthogonal();
185 TVector3 yp = currentWireDir.Cross(xp);
190 rot(0,0) =
xp.X(); rot(0,1) = yp.X(); rot(0,2) = currentWireDir.X();
191 rot(1,0) =
xp.Y(); rot(1,1) = yp.Y(); rot(1,2) = currentWireDir.Y();
192 rot(2,0) =
xp.Z(); rot(2,1) = yp.Z(); rot(2,2) = currentWireDir.Z();
195 TVectorD smearVec(3);
196 smearVec(0) = resolution_;
197 smearVec(1) = resolution_;
198 smearVec(2) = resolutionWire_;
201 hitCoords(0) += gRandom->Gaus(0, smearVec(0));
202 hitCoords(1) += gRandom->Gaus(0, smearVec(1));
203 hitCoords(2) += gRandom->Gaus(0, smearVec(2));
208 hitCov.Similarity(rot);
212 retVal.push_back(measurement);
216 case StripU:
case StripV:
case StripUV : {
217 if (debug_) std::cerr <<
"create StripHit" << std::endl;
220 vU = planeNorm.Cross(z);
221 vV = (planeNorm.Cross(z)).Cross(planeNorm);
224 TVectorD hitCoords(1);
226 hitCoords(0) = gRandom->Uniform(-outlierRange_, outlierRange_);
228 hitCoords(0) = gRandom->Gaus(0,resolution_);
230 TMatrixDSym hitCov(1);
231 hitCov(0,0) = resolution_*resolution_;
237 retVal.push_back(measurement);
240 if (type == StripUV) {
242 hitCoords(0) = gRandom->Uniform(-outlierRange_, outlierRange_);
244 hitCoords(0) = gRandom->Gaus(0,resolution_);
246 hitCov(0,0) = resolution_*resolution_;
251 retVal.push_back(measurement);
257 if (debug_) std::cerr <<
"create WireHit" << std::endl;
260 wirePerp.SetMag(gRandom->Uniform(outlierRange_));
263 TVectorD hitCoords(7);
264 hitCoords(0) = (point-wirePerp-currentWireDir).X();
265 hitCoords(1) = (point-wirePerp-currentWireDir).Y();
266 hitCoords(2) = (point-wirePerp-currentWireDir).Z();
268 hitCoords(3) = (point-wirePerp+currentWireDir).X();
269 hitCoords(4) = (point-wirePerp+currentWireDir).Y();
270 hitCoords(5) = (point-wirePerp+currentWireDir).Z();
273 hitCoords(6) = gRandom->Uniform(outlierRange_);
275 hitCoords(6) = gRandom->Gaus(wirePerp.Mag(), resolution_);
277 TMatrixDSym hitCov(7);
278 hitCov(6,6) = resolution_*resolution_;
282 if (idealLRResolution_){
286 retVal.push_back(measurement);
291 if (debug_) std::cerr <<
"create WirePointHit" << std::endl;
294 wirePerp.SetMag(gRandom->Uniform(outlierRange_));
297 TVectorD hitCoords(8);
298 hitCoords(0) = (point-wirePerp-currentWireDir).X();
299 hitCoords(1) = (point-wirePerp-currentWireDir).Y();
300 hitCoords(2) = (point-wirePerp-currentWireDir).Z();
302 hitCoords(3) = (point-wirePerp+currentWireDir).X();
303 hitCoords(4) = (point-wirePerp+currentWireDir).Y();
304 hitCoords(5) = (point-wirePerp+currentWireDir).Z();
307 hitCoords(6) = gRandom->Uniform(outlierRange_);
308 hitCoords(7) = gRandom->Uniform(currentWireDir.Mag()-outlierRange_, currentWireDir.Mag()+outlierRange_);
311 hitCoords(6) = gRandom->Gaus(wirePerp.Mag(), resolution_);
312 hitCoords(7) = gRandom->Gaus(currentWireDir.Mag(), resolutionWire_);
316 TMatrixDSym hitCov(8);
317 hitCov(6,6) = resolution_*resolution_;
318 hitCov(7,7) = resolutionWire_*resolutionWire_;
321 if (idealLRResolution_){
325 retVal.push_back(measurement);
330 std::cerr <<
"measurement type not defined!" << std::endl;
339 void MeasurementCreator::reset() {
341 measurementCounter_ = 0;