File indexing completed on 2024-04-06 12:27:40
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "RecoPPS/Local/interface/TotemRPLocalTrackFitterAlgorithm.h"
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include "TMatrixD.h"
0015
0016
0017
0018 using namespace std;
0019 using namespace edm;
0020
0021
0022
0023 TotemRPLocalTrackFitterAlgorithm::TotemRPLocalTrackFitterAlgorithm(const edm::ParameterSet &) {}
0024
0025
0026
0027 void TotemRPLocalTrackFitterAlgorithm::reset() { det_data_map_.clear(); }
0028
0029
0030
0031 TotemRPLocalTrackFitterAlgorithm::RPDetCoordinateAlgebraObjs
0032 TotemRPLocalTrackFitterAlgorithm::prepareReconstAlgebraData(unsigned int det_id, const CTPPSGeometry &tot_rp_geom) {
0033 RPDetCoordinateAlgebraObjs det_algebra_obj;
0034
0035 det_algebra_obj.centre_of_det_global_position_ = convert3vector(tot_rp_geom.sensorTranslation(det_id));
0036
0037 TVector3 rd_dir = convert3vector(tot_rp_geom.localToGlobalDirection(det_id, rp_topology_.GetStripReadoutAxisDir()));
0038
0039 TVector2 v(rd_dir.X(), rd_dir.Y());
0040 det_algebra_obj.readout_direction_ = v.Unit();
0041 det_algebra_obj.rec_u_0_ = 0.0;
0042 det_algebra_obj.available_ = true;
0043 det_algebra_obj.rec_u_0_ =
0044 -(det_algebra_obj.readout_direction_ * det_algebra_obj.centre_of_det_global_position_.XYvector());
0045
0046 return det_algebra_obj;
0047 }
0048
0049
0050
0051 TotemRPLocalTrackFitterAlgorithm::RPDetCoordinateAlgebraObjs *TotemRPLocalTrackFitterAlgorithm::getDetAlgebraData(
0052 unsigned int det_id, const CTPPSGeometry &tot_rp_geom) {
0053 auto it = det_data_map_.find(det_id);
0054 if (it != det_data_map_.end()) {
0055 return &(it->second);
0056 } else {
0057 det_data_map_[det_id] = prepareReconstAlgebraData(det_id, tot_rp_geom);
0058 return &det_data_map_[det_id];
0059 }
0060 }
0061
0062
0063
0064 bool TotemRPLocalTrackFitterAlgorithm::fitTrack(const edm::DetSetVector<TotemRPRecHit> &hits,
0065 double z_0,
0066 const CTPPSGeometry &tot_geom,
0067 TotemRPLocalTrack &fitted_track) {
0068 fitted_track.setValid(false);
0069
0070
0071 struct HitWithAlg {
0072 unsigned int detId;
0073 const TotemRPRecHit *hit;
0074 RPDetCoordinateAlgebraObjs *alg;
0075 };
0076
0077 vector<HitWithAlg> applicable_hits;
0078
0079 for (auto &ds : hits) {
0080 unsigned int detId = ds.detId();
0081
0082 for (auto &h : ds) {
0083 RPDetCoordinateAlgebraObjs *alg = getDetAlgebraData(detId, tot_geom);
0084 if (alg->available_)
0085 applicable_hits.push_back({detId, &h, alg});
0086 }
0087 }
0088
0089 if (applicable_hits.size() < 5)
0090 return false;
0091
0092 TMatrixD H(applicable_hits.size(), 4);
0093 TVectorD V(applicable_hits.size());
0094 TVectorD V_inv(applicable_hits.size());
0095 TVectorD U(applicable_hits.size());
0096
0097 for (unsigned int i = 0; i < applicable_hits.size(); ++i) {
0098 RPDetCoordinateAlgebraObjs *alg_obj = applicable_hits[i].alg;
0099
0100 H(i, 0) = alg_obj->readout_direction_.X();
0101 H(i, 1) = alg_obj->readout_direction_.Y();
0102 double delta_z = alg_obj->centre_of_det_global_position_.Z() - z_0;
0103 H(i, 2) = alg_obj->readout_direction_.X() * delta_z;
0104 H(i, 3) = alg_obj->readout_direction_.Y() * delta_z;
0105 double var = applicable_hits[i].hit->sigma();
0106 var *= var;
0107 V[i] = var;
0108 V_inv[i] = 1.0 / var;
0109 U[i] = applicable_hits[i].hit->position() - alg_obj->rec_u_0_;
0110 }
0111
0112 TMatrixD H_T_V_inv(TMatrixD::kTransposed, H);
0113 multiplyByDiagonalInPlace(H_T_V_inv, V_inv);
0114 TMatrixD V_a(H_T_V_inv);
0115 TMatrixD V_a_mult(V_a, TMatrixD::kMult, H);
0116 try {
0117 V_a_mult.Invert();
0118 } catch (cms::Exception &e) {
0119 LogError("TotemRPLocalTrackFitterAlgorithm") << "Error in TotemRPLocalTrackFitterAlgorithm::fitTrack > "
0120 << "Fit matrix is singular. Skipping.";
0121 return false;
0122 }
0123
0124 TMatrixD u_to_a(V_a_mult, TMatrixD::kMult, H_T_V_inv);
0125 TVectorD a(U);
0126 a *= u_to_a;
0127
0128 fitted_track.setZ0(z_0);
0129 fitted_track.setParameterVector(a);
0130 fitted_track.setCovarianceMatrix(V_a_mult);
0131
0132 double Chi_2 = 0;
0133 for (unsigned int i = 0; i < applicable_hits.size(); ++i) {
0134 RPDetCoordinateAlgebraObjs *alg_obj = applicable_hits[i].alg;
0135 TVector2 readout_dir = alg_obj->readout_direction_;
0136 double det_z = alg_obj->centre_of_det_global_position_.Z();
0137 double sigma_str = applicable_hits[i].hit->sigma();
0138 double sigma_str_2 = sigma_str * sigma_str;
0139 TVector2 fited_det_xy_point = fitted_track.trackPoint(det_z);
0140 double U_readout = applicable_hits[i].hit->position() - alg_obj->rec_u_0_;
0141 double U_fited = (readout_dir *= fited_det_xy_point);
0142 double residual = U_fited - U_readout;
0143 TMatrixD V_T_Cov_X_Y(1, 2);
0144 V_T_Cov_X_Y(0, 0) = readout_dir.X();
0145 V_T_Cov_X_Y(0, 1) = readout_dir.Y();
0146 TMatrixD V_T_Cov_X_Y_mult(V_T_Cov_X_Y, TMatrixD::kMult, fitted_track.trackPointInterpolationCovariance(det_z));
0147 double fit_strip_var = V_T_Cov_X_Y_mult(0, 0) * readout_dir.X() + V_T_Cov_X_Y_mult(0, 1) * readout_dir.Y();
0148 double pull_normalization = sqrt(sigma_str_2 - fit_strip_var);
0149 double pull = residual / pull_normalization;
0150
0151 Chi_2 += residual * residual / sigma_str_2;
0152
0153 TotemRPLocalTrack::FittedRecHit hit_point(
0154 *(applicable_hits[i].hit), TVector3(fited_det_xy_point.X(), fited_det_xy_point.Y(), det_z), residual, pull);
0155 fitted_track.addHit(applicable_hits[i].detId, hit_point);
0156 }
0157
0158 fitted_track.setChiSquared(Chi_2);
0159 fitted_track.setValid(true);
0160 return true;
0161 }
0162
0163
0164
0165 void TotemRPLocalTrackFitterAlgorithm::multiplyByDiagonalInPlace(TMatrixD &mt, const TVectorD &diag) {
0166 for (int i = 0; i < mt.GetNrows(); ++i) {
0167 for (int j = 0; j < mt.GetNcols(); ++j) {
0168 mt[i][j] *= diag[j];
0169 }
0170 }
0171 }
