utilities/src/MeasurementCreator.cc

0001 /* Copyright 2008-2010, Technische Universitaet Muenchen,
0002    Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
0003
0004    This file is part of GENFIT.
0005
0006    GENFIT is free software: you can redistribute it and/or modify
0007    it under the terms of the GNU Lesser General Public License as published
0008    by the Free Software Foundation, either version 3 of the License, or
0009    (at your option) any later version.
0010
0011    GENFIT is distributed in the hope that it will be useful,
0012    but WITHOUT ANY WARRANTY; without even the implied warranty of
0013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014    GNU Lesser General Public License for more details.
0015
0016    You should have received a copy of the GNU Lesser General Public License
0017    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0018 */
0019
0020 #include "MeasurementCreator.h"
0021
0022 #include <iostream>
0023
0024 #include <PlanarMeasurement.h>
0025 #include <ProlateSpacepointMeasurement.h>
0026 #include <SpacepointMeasurement.h>
0027 #include <WireMeasurement.h>
0028 #include <WirePointMeasurement.h>
0029
0030 #include <TRandom.h>
0031 #include <TMath.h>
0032
0033 #include <assert.h>
0034 #include <math.h>
0035
0036 namespace genfit {
0037
0038 MeasurementCreator::MeasurementCreator() :
0039     trackModel_(nullptr),
0040     resolution_(0.01),
0041     resolutionWire_(0.1),
0042     outlierProb_(0),
0043     outlierRange_(2),
0044     thetaDetPlane_(90),
0045     phiDetPlane_(0),
0046     wireCounter_(0),
0047     wireDir_(0.,0.,1.),
0048     minDrift_(0),
0049     maxDrift_(2),
0050     idealLRResolution_(true),
0051     useSkew_(false),
0052     skewAngle_(5),
0053     nSuperLayer_(5),
0054     measurementCounter_(0),
0055     debug_(false)
0056 {
0057   ;
0058 }
0059
0060
0061 std::vector<genfit::AbsMeasurement*> MeasurementCreator::create(eMeasurementType type, double tracklength, bool& outlier, int& lr) {
0062
0063   outlier = false;
0064   lr = 0;
0065   std::vector<AbsMeasurement*> retVal;
0066   genfit::AbsMeasurement* measurement;
0067
0068   TVector3 point, dir;
0069   trackModel_->getPosDir(tracklength, point, dir);
0070
0071
0072   TVector3 planeNorm(dir);
0073   planeNorm.SetTheta(thetaDetPlane_*TMath::Pi()/180);
0074   planeNorm.SetPhi(planeNorm.Phi()+phiDetPlane_);
0075   static const TVector3 z(0,0,1);
0076   static const TVector3 x(1,0,0);
0077
0078
0079   TVector3 currentWireDir(wireDir_);
0080   TVector3 wirePerp;
0081
0082   if (type == Wire ||
0083       type == WirePoint){
0084
0085     // skew layers
0086     if (useSkew_ && (int)((double)wireCounter_/(double)nSuperLayer_)%2 == 1) {
0087       TVector3 perp(wireDir_.Cross(dir));
0088       if ((int)((double)wireCounter_/(double)nSuperLayer_)%4 == 1){
0089         currentWireDir.Rotate(skewAngle_*TMath::Pi()/180, wireDir_.Cross(perp));
0090       }
0091       else currentWireDir.Rotate(-skewAngle_*TMath::Pi()/180, wireDir_.Cross(perp));
0092     }
0093     currentWireDir.SetMag(1.);
0094
0095     // left/right
0096     lr = 1;
0097     wirePerp = dir.Cross(currentWireDir);
0098     if (gRandom->Uniform(-1,1) >= 0) {
0099       wirePerp *= -1.;
0100       lr = -1;
0101     }
0102     wirePerp.SetMag(gRandom->Uniform(minDrift_, maxDrift_));
0103   }
0104
0105   if (outlierProb_ > gRandom->Uniform(1.)) {
0106     outlier = true;
0107     if(debug_)  std::cerr << "create outlier" << std::endl;
0108   }
0109
0110
0111   switch(type){
0112     case Pixel: {
0113       if (debug_) std::cerr << "create PixHit" << std::endl;
0114
0115       genfit::SharedPlanePtr plane(new genfit::DetPlane(point, planeNorm.Cross(z), (planeNorm.Cross(z)).Cross(planeNorm)));
0116
0117       TVectorD hitCoords(2);
0118       if (outlier) {
0119         hitCoords(0) = gRandom->Uniform(-outlierRange_, outlierRange_);
0120         hitCoords(1) = gRandom->Uniform(-outlierRange_, outlierRange_);
0121       }
0122       else {
0123         hitCoords(0) = gRandom->Gaus(0,resolution_);
0124         hitCoords(1) = gRandom->Gaus(0,resolution_);
0125       }
0126
0127       TMatrixDSym hitCov(2);
0128       hitCov(0,0) = resolution_*resolution_;
0129       hitCov(1,1) = resolution_*resolution_;
0130
0131       measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
0132       static_cast<genfit::PlanarMeasurement*>(measurement)->setPlane(plane, measurementCounter_);
0133       retVal.push_back(measurement);
0134     }
0135     break;
0136
0137     case Spacepoint: {
0138       if (debug_) std::cerr << "create SpacepointHit" << std::endl;
0139
0140       TVectorD hitCoords(3);
0141       if (outlier) {
0142         hitCoords(0) = gRandom->Uniform(point.X()-outlierRange_, point.X()+outlierRange_);
0143         hitCoords(1) = gRandom->Uniform(point.Y()-outlierRange_, point.Y()+outlierRange_);
0144         hitCoords(2) = gRandom->Uniform(point.Z()-outlierRange_, point.Z()+outlierRange_);
0145       }
0146       else {
0147         hitCoords(0) = gRandom->Gaus(point.X(),resolution_);
0148         hitCoords(1) = gRandom->Gaus(point.Y(),resolution_);
0149         hitCoords(2) = gRandom->Gaus(point.Z(),resolution_);
0150       }
0151
0152       TMatrixDSym hitCov(3);
0153       hitCov(0,0) = resolution_*resolution_;
0154       hitCov(1,1) = resolution_*resolution_;
0155       hitCov(2,2) = resolution_*resolution_;
0156
0157       measurement = new genfit::SpacepointMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
0158       retVal.push_back(measurement);
0159     }
0160     break;
0161
0162     case ProlateSpacepoint: {
0163       if (debug_) std::cerr << "create ProlateSpacepointHit" << std::endl;
0164
0165       TVectorD hitCoords(3);
0166       if (outlier) {
0167         hitCoords(0) = gRandom->Uniform(point.X()-outlierRange_, point.X()+outlierRange_);
0168         hitCoords(1) = gRandom->Uniform(point.Y()-outlierRange_, point.Y()+outlierRange_);
0169         hitCoords(2) = gRandom->Uniform(point.Z()-outlierRange_, point.Z()+outlierRange_);
0170       }
0171       else {
0172         hitCoords(0) = point.X();
0173         hitCoords(1) = point.Y();
0174         hitCoords(2) = point.Z();
0175       }
0176
0177       TMatrixDSym hitCov(3);
0178       hitCov(0,0) = resolution_*resolution_;
0179       hitCov(1,1) = resolution_*resolution_;
0180       hitCov(2,2) = resolutionWire_*resolutionWire_;
0181
0182       // rotation matrix
0183       TVector3 xp = currentWireDir.Orthogonal();
0184       xp.SetMag(1);
0185       TVector3 yp = currentWireDir.Cross(xp);
0186       yp.SetMag(1);
0187
0188       TMatrixD rot(3,3);
0189
0190       rot(0,0) = xp.X();  rot(0,1) = yp.X();  rot(0,2) = currentWireDir.X();
0191       rot(1,0) = xp.Y();  rot(1,1) = yp.Y();  rot(1,2) = currentWireDir.Y();
0192       rot(2,0) = xp.Z();  rot(2,1) = yp.Z();  rot(2,2) = currentWireDir.Z();
0193
0194       // smear
0195       TVectorD smearVec(3);
0196       smearVec(0) = resolution_;
0197       smearVec(1) = resolution_;
0198       smearVec(2) = resolutionWire_;
0199       smearVec *= rot;
0200       if (!outlier) {
0201         hitCoords(0) += gRandom->Gaus(0, smearVec(0));
0202         hitCoords(1) += gRandom->Gaus(0, smearVec(1));
0203         hitCoords(2) += gRandom->Gaus(0, smearVec(2));
0204       }
0205
0206
0207       // rotate cov
0208       hitCov.Similarity(rot);
0209
0210       measurement = new genfit::ProlateSpacepointMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
0211       static_cast<genfit::ProlateSpacepointMeasurement*>(measurement)->setLargestErrorDirection(currentWireDir);
0212       retVal.push_back(measurement);
0213     }
0214     break;
0215
0216     case StripU: case StripV: case StripUV : {
0217       if (debug_) std::cerr << "create StripHit" << std::endl;
0218
0219       TVector3 vU, vV;
0220       vU = planeNorm.Cross(z);
0221       vV = (planeNorm.Cross(z)).Cross(planeNorm);
0222       genfit::SharedPlanePtr plane(new genfit::DetPlane(point, vU, vV));
0223
0224       TVectorD hitCoords(1);
0225       if (outlier)
0226         hitCoords(0) = gRandom->Uniform(-outlierRange_, outlierRange_);
0227       else
0228         hitCoords(0) = gRandom->Gaus(0,resolution_);
0229
0230       TMatrixDSym hitCov(1);
0231       hitCov(0,0) = resolution_*resolution_;
0232
0233       measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
0234       static_cast<genfit::PlanarMeasurement*>(measurement)->setPlane(plane, measurementCounter_);
0235       if (type == StripV)
0236         static_cast<genfit::PlanarMeasurement*>(measurement)->setStripV();
0237       retVal.push_back(measurement);
0238
0239
0240       if (type == StripUV) {
0241         if (outlier)
0242           hitCoords(0) = gRandom->Uniform(-outlierRange_, outlierRange_);
0243         else
0244           hitCoords(0) = gRandom->Gaus(0,resolution_);
0245
0246         hitCov(0,0) = resolution_*resolution_;
0247
0248         measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
0249         static_cast<genfit::PlanarMeasurement*>(measurement)->setPlane(plane, measurementCounter_);
0250         static_cast<genfit::PlanarMeasurement*>(measurement)->setStripV();
0251         retVal.push_back(measurement);
0252       }
0253     }
0254     break;
0255
0256     case Wire: {
0257       if (debug_) std::cerr << "create WireHit" << std::endl;
0258
0259       if (outlier) {
0260         wirePerp.SetMag(gRandom->Uniform(outlierRange_));
0261       }
0262
0263       TVectorD hitCoords(7);
0264       hitCoords(0) = (point-wirePerp-currentWireDir).X();
0265       hitCoords(1) = (point-wirePerp-currentWireDir).Y();
0266       hitCoords(2) = (point-wirePerp-currentWireDir).Z();
0267
0268       hitCoords(3) = (point-wirePerp+currentWireDir).X();
0269       hitCoords(4) = (point-wirePerp+currentWireDir).Y();
0270       hitCoords(5) = (point-wirePerp+currentWireDir).Z();
0271
0272       if (outlier)
0273         hitCoords(6) = gRandom->Uniform(outlierRange_);
0274       else
0275         hitCoords(6) = gRandom->Gaus(wirePerp.Mag(), resolution_);
0276
0277       TMatrixDSym hitCov(7);
0278       hitCov(6,6) = resolution_*resolution_;
0279
0280
0281       measurement = new genfit::WireMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
0282       if (idealLRResolution_){
0283         static_cast<genfit::WireMeasurement*>(measurement)->setLeftRightResolution(lr);
0284       }
0285       ++wireCounter_;
0286       retVal.push_back(measurement);
0287     }
0288     break;
0289
0290     case WirePoint: {
0291       if (debug_) std::cerr << "create WirePointHit" << std::endl;
0292
0293       if (outlier) {
0294         wirePerp.SetMag(gRandom->Uniform(outlierRange_));
0295       }
0296
0297       TVectorD hitCoords(8);
0298       hitCoords(0) = (point-wirePerp-currentWireDir).X();
0299       hitCoords(1) = (point-wirePerp-currentWireDir).Y();
0300       hitCoords(2) = (point-wirePerp-currentWireDir).Z();
0301
0302       hitCoords(3) = (point-wirePerp+currentWireDir).X();
0303       hitCoords(4) = (point-wirePerp+currentWireDir).Y();
0304       hitCoords(5) = (point-wirePerp+currentWireDir).Z();
0305
0306       if (outlier) {
0307         hitCoords(6) = gRandom->Uniform(outlierRange_);
0308         hitCoords(7) = gRandom->Uniform(currentWireDir.Mag()-outlierRange_, currentWireDir.Mag()+outlierRange_);
0309       }
0310       else {
0311         hitCoords(6) = gRandom->Gaus(wirePerp.Mag(), resolution_);
0312         hitCoords(7) = gRandom->Gaus(currentWireDir.Mag(), resolutionWire_);
0313       }
0314
0315
0316       TMatrixDSym hitCov(8);
0317       hitCov(6,6) = resolution_*resolution_;
0318       hitCov(7,7) = resolutionWire_*resolutionWire_;
0319
0320       measurement = new genfit::WirePointMeasurement(hitCoords, hitCov, int(type), measurementCounter_, nullptr);
0321       if (idealLRResolution_){
0322         static_cast<genfit::WirePointMeasurement*>(measurement)->setLeftRightResolution(lr);
0323       }
0324       ++wireCounter_;
0325       retVal.push_back(measurement);
0326     }
0327     break;
0328
0329     default:
0330       std::cerr << "measurement type not defined!" << std::endl;
0331       exit(0);
0332   }
0333
0334   return retVal;
0335
0336 }
0337
0338
0339 void MeasurementCreator::reset() {
0340   wireCounter_ = 0;
0341   measurementCounter_ = 0;
0342 }
0343
0344 } /* End of namespace genfit */