CMS 3D CMS Logo

LocalTrackFitter.cc
Go to the documentation of this file.
1 /****************************************************************************
2 * Authors:
3 * Jan Kašpar (jan.kaspar@gmail.com)
4 ****************************************************************************/
5 
7 
9 
11 
14 
15 #include <set>
16 
17 #include "TMatrixD.h"
18 #include "TVectorD.h"
19 
20 using namespace std;
21 
22 //----------------------------------------------------------------------------------------------------
23 
25  : verbosity(ps.getUntrackedParameter<unsigned int>("verbosity", 0)),
26  minimumHitsPerProjectionPerRP(ps.getParameter<unsigned int>("minimumHitsPerProjectionPerRP")),
27  maxResidualToSigma(ps.getParameter<double>("maxResidualToSigma")) {}
28 
29 //----------------------------------------------------------------------------------------------------
30 
32  if (verbosity > 5)
33  printf(">> LocalTrackFitter::Fit\n");
34 
35  bool selectionChanged = true;
36  unsigned int loopCounter = 0;
37  while (selectionChanged) {
38  // fit/outlier-removal loop
39  while (selectionChanged) {
40  if (verbosity > 5)
41  printf("* fit loop %u\n", loopCounter++);
42 
43  bool fitFailed = false;
44  fitAndRemoveOutliers(selection, geometry, trackFit, fitFailed, selectionChanged);
45 
46  if (fitFailed) {
47  if (verbosity > 5)
48  printf("\tFIT FAILED\n");
49  return false;
50  }
51  }
52 
53  // remove pots with too few active planes
54  removeInsufficientPots(selection, selectionChanged);
55  }
56 
57  return true;
58 }
59 
60 //----------------------------------------------------------------------------------------------------
61 
64  LocalTrackFit &trackFit,
65  bool &failed,
66  bool &selectionChanged) const {
67  if (verbosity > 5)
68  printf(" - LocalTrackFitter::FitAndRemoveOutliers\n");
69 
70  if (selection.empty()) {
71  failed = true;
72  return;
73  }
74 
75  // check if input size is sufficient
76  if (selection.size() < 4) {
77  failed = true;
78  return;
79  }
80 
81  // build matrices and vectors
82  TMatrixD A(selection.size(), 4);
83  TMatrixD Vi(selection.size(), selection.size());
84  TVectorD measVec(selection.size());
85  unsigned int j = 0;
86  for (auto it = selection.begin(); it != selection.end(); ++it, ++j) {
87  const unsigned int &detId = it->id;
88 
89  const DetGeometry &d = geometry.get(detId);
90  const auto &dirData = d.getDirectionData(it->dirIdx);
91 
92  A(j, 0) = it->z * dirData.dx;
93  A(j, 1) = dirData.dx;
94  A(j, 2) = it->z * dirData.dy;
95  A(j, 3) = dirData.dy;
96 
97  measVec(j) = it->position + dirData.s - (it->z - d.z) * dirData.dz; // in mm
98 
99  Vi(j, j) = 1. / it->sigma / it->sigma;
100  }
101 
102  // evaluate local track parameter estimates (h stands for hat)
103  TMatrixD AT(4, selection.size());
104  AT.Transpose(A);
105  TMatrixD ATViA(4, 4);
106  ATViA = AT * Vi * A;
107  TMatrixD ATViAI(ATViA);
108 
109  try {
110  ATViAI = ATViA.Invert();
111  } catch (cms::Exception &e) {
112  failed = true;
113  return;
114  }
115 
116  TVectorD theta(4);
117  theta = ATViAI * AT * Vi * measVec;
118 
119  // residuals
120  TVectorD R(measVec);
121  R -= A * theta;
122 
123  // save results to trackFit
124  trackFit.ax = theta(0);
125  trackFit.bx = theta(1);
126  trackFit.ay = theta(2);
127  trackFit.by = theta(3);
128  trackFit.z0 = geometry.z0;
129  trackFit.ndf = selection.size() - 4;
130  trackFit.chi_sq = 0;
131  for (int i = 0; i < R.GetNrows(); i++)
132  trackFit.chi_sq += R(i) * R(i) * Vi(i, i);
133 
134  if (verbosity > 5) {
135  printf(" ax = %.3f mrad, bx = %.4f mm, ay = %.3f mrad, by = %.4f mm, z0 = %.3f mm\n",
136  trackFit.ax * 1E3,
137  trackFit.bx,
138  trackFit.ay * 1E3,
139  trackFit.by,
140  trackFit.z0);
141  printf(" ndof = %i, chi^2/ndof/si^2 = %.3f\n", trackFit.ndf, trackFit.chi_sq / trackFit.ndf);
142  }
143 
144  // check residuals
145  selectionChanged = false;
146  TVectorD interpolation(A * theta);
147  j = 0;
148  for (auto it = selection.begin(); it != selection.end(); ++j) {
149  if (verbosity > 5) {
150  printf(" %2u, ", j);
151  printId(it->id);
152  printf(", dirIdx=%u: interpol = %+8.1f um, residual = %+6.1f um, residual / sigma = %+6.2f\n",
153  it->dirIdx,
154  interpolation[j] * 1E3,
155  R[j] * 1E3,
156  R[j] / it->sigma);
157  }
158 
159  double resToSigma = R[j] / it->sigma;
160  if (fabs(resToSigma) > maxResidualToSigma) {
161  selection.erase(it);
162  selectionChanged = true;
163  if (verbosity > 5)
164  printf(" Removed\n");
165  } else
166  ++it;
167  }
168 }
169 
170 //----------------------------------------------------------------------------------------------------
171 
172 void LocalTrackFitter::removeInsufficientPots(HitCollection &selection, bool &selectionChanged) const {
173  if (verbosity > 5)
174  printf(" - RemoveInsufficientPots\n");
175 
176  selectionChanged = false;
177 
178  // map: RP id -> (active planes in projection 1, active planes in projection 2)
179  map<unsigned int, pair<set<unsigned int>, set<unsigned int> > > planeMap;
180  for (auto it = selection.begin(); it != selection.end(); ++it) {
181  CTPPSDetId senId(it->id);
182  const unsigned int rpId = senId.rpId();
183 
184  if (senId.subdetId() == CTPPSDetId::sdTrackingStrip) {
185  const unsigned int plane = TotemRPDetId(it->id).plane();
186  if ((plane % 2) == 0)
187  planeMap[rpId].first.insert(senId);
188  else
189  planeMap[rpId].second.insert(senId);
190  }
191 
192  if (senId.subdetId() == CTPPSDetId::sdTrackingPixel) {
193  planeMap[rpId].first.insert(senId);
194  planeMap[rpId].second.insert(senId);
195  }
196  }
197 
198  // remove RPs with insufficient information
199  selectionChanged = false;
200 
201  for (auto it = planeMap.begin(); it != planeMap.end(); ++it) {
202  if (it->second.first.size() < minimumHitsPerProjectionPerRP ||
203  it->second.second.size() < minimumHitsPerProjectionPerRP) {
204  if (verbosity > 5)
205  printf("\tRP %u: projection1 = %lu, projection 2 = %lu\n",
206  it->first,
207  it->second.first.size(),
208  it->second.second.size());
209 
210  // remove all hits from that RP
211  for (auto dit = selection.begin(); dit != selection.end();) {
212  if (it->first == CTPPSDetId(dit->id).rpId()) {
213  if (verbosity > 5)
214  printf("\t\tremoving %u\n", dit->id);
215  selection.erase(dit);
216  selectionChanged = true;
217  } else {
218  dit++;
219  }
220  }
221  }
222  }
223 }
HIPAlignmentAlgorithm_cfi.verbosity
verbosity
Definition: HIPAlignmentAlgorithm_cfi.py:7
mps_fire.i
i
Definition: mps_fire.py:428
geometry
Definition: geometry.py:1
LocalTrackFit::ndf
signed int ndf
the number of degrees of freedom
Definition: LocalTrackFit.h:27
runEdmFileComparison.failed
failed
Definition: runEdmFileComparison.py:225
TotemRPDetId::plane
uint32_t plane() const
Definition: TotemRPDetId.h:51
year_2016_postTS2_cff.rpId
rpId
Definition: year_2016_postTS2_cff.py:23
LocalTrackFit::by
double by
Definition: LocalTrackFit.h:24
LocalTrackFitter.h
DetGeometry
A structure to hold relevant geometrical information about one detector/sensor.
Definition: AlignmentGeometry.h:19
CTPPSDetId::sdTrackingStrip
Definition: CTPPSDetId.h:44
LocalTrackFitter::fit
bool fit(HitCollection &, const AlignmentGeometry &, LocalTrackFit &) const
Definition: LocalTrackFitter.cc:30
corrVsCorr.selection
selection
main part
Definition: corrVsCorr.py:100
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
LocalTrackFitter::fitAndRemoveOutliers
void fitAndRemoveOutliers(HitCollection &, const AlignmentGeometry &, LocalTrackFit &, bool &failed, bool &selectionChanged) const
Definition: LocalTrackFitter.cc:61
LocalTrackFitter::verbosity
unsigned int verbosity
verbosity level
Definition: LocalTrackFitter.h:35
LocalTrackFitter::removeInsufficientPots
void removeInsufficientPots(HitCollection &, bool &selectionChanged) const
removes the hits of pots with too few planes active
Definition: LocalTrackFitter.cc:171
CTPPSDetId::sdTrackingPixel
Definition: CTPPSDetId.h:44
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum)
Definition: DetId.h:48
LocalTrackFitter::LocalTrackFitter
LocalTrackFitter()
dummy constructor (not to be used)
Definition: LocalTrackFitter.h:22
edm::ParameterSet
Definition: ParameterSet.h:47
LocalTrackFit::chi_sq
double chi_sq
the residual sum of squares
Definition: LocalTrackFit.h:30
LocalTrackFitter::minimumHitsPerProjectionPerRP
unsigned int minimumHitsPerProjectionPerRP
minimum of hits to accept data from a RP
Definition: LocalTrackFitter.h:38
CTPPSDetId
Base class for CTPPS detector IDs.
Definition: CTPPSDetId.h:31
createfilelist.int
int
Definition: createfilelist.py:10
Utilities.h
LocalTrackFit::ax
double ax
slopes in rad
Definition: LocalTrackFit.h:21
MaterialEffects_cfi.A
A
Definition: MaterialEffects_cfi.py:11
std
Definition: JetResolutionObject.h:76
LocalTrackFit::bx
double bx
intercepts in mm
Definition: LocalTrackFit.h:24
printId
void printId(unsigned int id)
Definition: Utilities.cc:25
HitCollection
std::vector< Hit > HitCollection
Definition: HitCollection.h:34
LocalTrackFit
Local (linear) track description (or a fit result). Uses global reference system.
Definition: LocalTrackFit.h:14
TotemRPDetId.h
AlignmentGeometry
Definition: AlignmentGeometry.h:55
CTPPSDetId.h
ztail.d
d
Definition: ztail.py:151
cms::Exception
Definition: Exception.h:70
ParameterSet.h
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
CTPPSDetId::rpId
CTPPSDetId rpId() const
Definition: CTPPSDetId.h:82
LocalTrackFit::z0
double z0
the point where intercepts are measured, in mm
Definition: LocalTrackFit.h:18
LocalTrackFitter::maxResidualToSigma
double maxResidualToSigma
hits with higher ratio residual/sigma will be dropped
Definition: LocalTrackFitter.h:41
A
dttmaxenums::R
Definition: DTTMax.h:29
TotemRPDetId
Detector ID class for TOTEM Si strip detectors.
Definition: TotemRPDetId.h:29
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
LocalTrackFit::ay
double ay
Definition: LocalTrackFit.h:21