Belle II Software  release-05-01-25
MeasurementCreator.cc
1 /* Copyright 2008-2010, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "MeasurementCreator.h"
21 
22 #include <iostream>
23 
24 #include <PlanarMeasurement.h>
25 #include <ProlateSpacepointMeasurement.h>
26 #include <SpacepointMeasurement.h>
27 #include <WireMeasurement.h>
28 #include <WirePointMeasurement.h>
29 
30 #include <TRandom.h>
31 #include <TMath.h>
32 
33 #include <assert.h>
34 #include <math.h>
35 
36 namespace genfit {
37 
38 MeasurementCreator::MeasurementCreator() :
39  trackModel_(nullptr),
40  resolution_(0.01),
41  resolutionWire_(0.1),
42  outlierProb_(0),
43  outlierRange_(2),
44  thetaDetPlane_(90),
45  phiDetPlane_(0),
46  wireCounter_(0),
47  wireDir_(0.,0.,1.),
48  minDrift_(0),
49  maxDrift_(2),
50  idealLRResolution_(true),
51  useSkew_(false),
52  skewAngle_(5),
53  nSuperLayer_(5),
54  measurementCounter_(0),
55  debug_(false)
56 {
57  ;
58 }
59 
60 
61 std::vector<genfit::AbsMeasurement*> MeasurementCreator::create(eMeasurementType type, double tracklength, bool& outlier, int& lr) {
62 
63  outlier = false;
64  lr = 0;
65  std::vector<AbsMeasurement*> retVal;
66  genfit::AbsMeasurement* measurement;
67 
68  TVector3 point, dir;
69  trackModel_->getPosDir(tracklength, point, dir);
70 
71 
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);
77 
78 
79  TVector3 currentWireDir(wireDir_);
80  TVector3 wirePerp;
81 
82  if (type == Wire ||
83  type == WirePoint){
84 
85  // skew layers
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));
90  }
91  else currentWireDir.Rotate(-skewAngle_*TMath::Pi()/180, wireDir_.Cross(perp));
92  }
93  currentWireDir.SetMag(1.);
94 
95  // left/right
96  lr = 1;
97  wirePerp = dir.Cross(currentWireDir);
98  if (gRandom->Uniform(-1,1) >= 0) {
99  wirePerp *= -1.;
100  lr = -1;
101  }
102  wirePerp.SetMag(gRandom->Uniform(minDrift_, maxDrift_));
103  }
104 
105  if (outlierProb_ > gRandom->Uniform(1.)) {
106  outlier = true;
107  if(debug_) std::cerr << "create outlier" << std::endl;
108  }
109 
110 
111  switch(type){
112  case Pixel: {
113  if (debug_) std::cerr << "create PixHit" << std::endl;
114 
115  genfit::SharedPlanePtr plane(new genfit::DetPlane(point, planeNorm.Cross(z), (planeNorm.Cross(z)).Cross(planeNorm)));
116 
117  TVectorD hitCoords(2);
118  if (outlier) {
119  hitCoords(0) = gRandom->Uniform(-outlierRange_, outlierRange_);
120  hitCoords(1) = gRandom->Uniform(-outlierRange_, outlierRange_);
121  }
122  else {
123  hitCoords(0) = gRandom->Gaus(0,resolution_);
124  hitCoords(1) = gRandom->Gaus(0,resolution_);
125  }
126 
127  TMatrixDSym hitCov(2);
128  hitCov(0,0) = resolution_*resolution_;
129  hitCov(1,1) = resolution_*resolution_;
130 
131  measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
132  static_cast<genfit::PlanarMeasurement*>(measurement)->setPlane(plane, measurementCounter_);
133  retVal.push_back(measurement);
134  }
135  break;
136 
137  case Spacepoint: {
138  if (debug_) std::cerr << "create SpacepointHit" << std::endl;
139 
140  TVectorD hitCoords(3);
141  if (outlier) {
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_);
145  }
146  else {
147  hitCoords(0) = gRandom->Gaus(point.X(),resolution_);
148  hitCoords(1) = gRandom->Gaus(point.Y(),resolution_);
149  hitCoords(2) = gRandom->Gaus(point.Z(),resolution_);
150  }
151 
152  TMatrixDSym hitCov(3);
153  hitCov(0,0) = resolution_*resolution_;
154  hitCov(1,1) = resolution_*resolution_;
155  hitCov(2,2) = resolution_*resolution_;
156 
157  measurement = new genfit::SpacepointMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
158  retVal.push_back(measurement);
159  }
160  break;
161 
162  case ProlateSpacepoint: {
163  if (debug_) std::cerr << "create ProlateSpacepointHit" << std::endl;
164 
165  TVectorD hitCoords(3);
166  if (outlier) {
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_);
170  }
171  else {
172  hitCoords(0) = point.X();
173  hitCoords(1) = point.Y();
174  hitCoords(2) = point.Z();
175  }
176 
177  TMatrixDSym hitCov(3);
178  hitCov(0,0) = resolution_*resolution_;
179  hitCov(1,1) = resolution_*resolution_;
180  hitCov(2,2) = resolutionWire_*resolutionWire_;
181 
182  // rotation matrix
183  TVector3 xp = currentWireDir.Orthogonal();
184  xp.SetMag(1);
185  TVector3 yp = currentWireDir.Cross(xp);
186  yp.SetMag(1);
187 
188  TMatrixD rot(3,3);
189 
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();
193 
194  // smear
195  TVectorD smearVec(3);
196  smearVec(0) = resolution_;
197  smearVec(1) = resolution_;
198  smearVec(2) = resolutionWire_;
199  smearVec *= rot;
200  if (!outlier) {
201  hitCoords(0) += gRandom->Gaus(0, smearVec(0));
202  hitCoords(1) += gRandom->Gaus(0, smearVec(1));
203  hitCoords(2) += gRandom->Gaus(0, smearVec(2));
204  }
205 
206 
207  // rotate cov
208  hitCov.Similarity(rot);
209 
210  measurement = new genfit::ProlateSpacepointMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
211  static_cast<genfit::ProlateSpacepointMeasurement*>(measurement)->setLargestErrorDirection(currentWireDir);
212  retVal.push_back(measurement);
213  }
214  break;
215 
216  case StripU: case StripV: case StripUV : {
217  if (debug_) std::cerr << "create StripHit" << std::endl;
218 
219  TVector3 vU, vV;
220  vU = planeNorm.Cross(z);
221  vV = (planeNorm.Cross(z)).Cross(planeNorm);
222  genfit::SharedPlanePtr plane(new genfit::DetPlane(point, vU, vV));
223 
224  TVectorD hitCoords(1);
225  if (outlier)
226  hitCoords(0) = gRandom->Uniform(-outlierRange_, outlierRange_);
227  else
228  hitCoords(0) = gRandom->Gaus(0,resolution_);
229 
230  TMatrixDSym hitCov(1);
231  hitCov(0,0) = resolution_*resolution_;
232 
233  measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
234  static_cast<genfit::PlanarMeasurement*>(measurement)->setPlane(plane, measurementCounter_);
235  if (type == StripV)
236  static_cast<genfit::PlanarMeasurement*>(measurement)->setStripV();
237  retVal.push_back(measurement);
238 
239 
240  if (type == StripUV) {
241  if (outlier)
242  hitCoords(0) = gRandom->Uniform(-outlierRange_, outlierRange_);
243  else
244  hitCoords(0) = gRandom->Gaus(0,resolution_);
245 
246  hitCov(0,0) = resolution_*resolution_;
247 
248  measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
249  static_cast<genfit::PlanarMeasurement*>(measurement)->setPlane(plane, measurementCounter_);
250  static_cast<genfit::PlanarMeasurement*>(measurement)->setStripV();
251  retVal.push_back(measurement);
252  }
253  }
254  break;
255 
256  case Wire: {
257  if (debug_) std::cerr << "create WireHit" << std::endl;
258 
259  if (outlier) {
260  wirePerp.SetMag(gRandom->Uniform(outlierRange_));
261  }
262 
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();
267 
268  hitCoords(3) = (point-wirePerp+currentWireDir).X();
269  hitCoords(4) = (point-wirePerp+currentWireDir).Y();
270  hitCoords(5) = (point-wirePerp+currentWireDir).Z();
271 
272  if (outlier)
273  hitCoords(6) = gRandom->Uniform(outlierRange_);
274  else
275  hitCoords(6) = gRandom->Gaus(wirePerp.Mag(), resolution_);
276 
277  TMatrixDSym hitCov(7);
278  hitCov(6,6) = resolution_*resolution_;
279 
280 
281  measurement = new genfit::WireMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
282  if (idealLRResolution_){
283  static_cast<genfit::WireMeasurement*>(measurement)->setLeftRightResolution(lr);
284  }
285  ++wireCounter_;
286  retVal.push_back(measurement);
287  }
288  break;
289 
290  case WirePoint: {
291  if (debug_) std::cerr << "create WirePointHit" << std::endl;
292 
293  if (outlier) {
294  wirePerp.SetMag(gRandom->Uniform(outlierRange_));
295  }
296 
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();
301 
302  hitCoords(3) = (point-wirePerp+currentWireDir).X();
303  hitCoords(4) = (point-wirePerp+currentWireDir).Y();
304  hitCoords(5) = (point-wirePerp+currentWireDir).Z();
305 
306  if (outlier) {
307  hitCoords(6) = gRandom->Uniform(outlierRange_);
308  hitCoords(7) = gRandom->Uniform(currentWireDir.Mag()-outlierRange_, currentWireDir.Mag()+outlierRange_);
309  }
310  else {
311  hitCoords(6) = gRandom->Gaus(wirePerp.Mag(), resolution_);
312  hitCoords(7) = gRandom->Gaus(currentWireDir.Mag(), resolutionWire_);
313  }
314 
315 
316  TMatrixDSym hitCov(8);
317  hitCov(6,6) = resolution_*resolution_;
318  hitCov(7,7) = resolutionWire_*resolutionWire_;
319 
320  measurement = new genfit::WirePointMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
321  if (idealLRResolution_){
322  static_cast<genfit::WirePointMeasurement*>(measurement)->setLeftRightResolution(lr);
323  }
324  ++wireCounter_;
325  retVal.push_back(measurement);
326  }
327  break;
328 
329  default:
330  std::cerr << "measurement type not defined!" << std::endl;
331  exit(0);
332  }
333 
334  return retVal;
335 
336 }
337 
338 
339 void MeasurementCreator::reset() {
340  wireCounter_ = 0;
341  measurementCounter_ = 0;
342 }
343 
344 } /* End of namespace genfit */
genfit::SharedPlanePtr
std::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
Definition: SharedPlanePtr.h:40
prepareAsicCrosstalkSimDB.xp
xp
location of points for the plot
Definition: prepareAsicCrosstalkSimDB.py:61
genfit::WirePointMeasurement
Class for measurements in wire detectors (Straw tubes and drift chambers) which can measure the coord...
Definition: WirePointMeasurement.h:51
genfit
Defines for I/O streams used for error and debug printing.
Definition: AlignablePXDRecoHit.h:19
genfit::PlanarMeasurement
Measurement class implementing a planar hit geometry (1 or 2D).
Definition: PlanarMeasurement.h:44
genfit::AbsMeasurement
Contains the measurement and covariance in raw detector coordinates.
Definition: AbsMeasurement.h:42
genfit::ProlateSpacepointMeasurement
Class for measurements implementing a space point hit geometry with a very prolate form of the covari...
Definition: ProlateSpacepointMeasurement.h:46
genfit::SpacepointMeasurement
Class for measurements implementing a space point hit geometry.
Definition: SpacepointMeasurement.h:46
genfit::DetPlane
Detector plane.
Definition: DetPlane.h:59
genfit::WireMeasurement
Class for measurements in wire detectors (Straw tubes and drift chambers) which do not measure the co...
Definition: WireMeasurement.h:52