Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MeasurementCreator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 */