ATLAS Offline Software
JetFitterHelper.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /***************************************************************************
6  JetFitterHelper.h - Description
7  -------------------
8 
9  begin : Februar 2007
10  authors: Giacinto Piacquadio (University of Freiburg),
11  Christian Weiser (University of Freiburg)
12  email : nicola.giacinto.piacquadio@cern.ch,
13  christian.weiser@cern.ch
14  changes: new!
15 
16  2007 (c) Atlas Detector Software
17 
18  Look at the header file for more information.
19 
20  ***************************************************************************/
21 
22 #include "TrkJetVxFitter/JetFitterHelper.h"
23 #include "VxVertex/RecVertex.h"
24 #include "VxJetVertex/VxJetCandidate.h"
25 #include "VxVertex/VxTrackAtVertex.h"
26 #include "VxJetVertex/VxVertexOnJetAxis.h"
27 #include "VxJetVertex/RecVertexPositions.h"
28 #include "VxJetVertex/JetVtxParamDefs.h"
29 
30 #include "TrkJetVxFitter/Utilities.h"
31 
32 
33 
34 namespace Trk
35 {
36 
37  JetFitterHelper::JetFitterHelper(const std::string& t, const std::string& n, const IInterface* p) :
38  AthAlgTool(t,n,p)
39  {
40  declareInterface< JetFitterHelper >(this) ;
41  }
42 
43 
44  JetFitterHelper::~JetFitterHelper() = default;
45 
46 
47  void JetFitterHelper::addTracksOfFirstVertexToSecondVertex(const VxVertexOnJetAxis & first,VxVertexOnJetAxis & second) const {
48 
49 
50  //disclaimer: no real copy of tracks, only copy of pointers...
51 
52  //to add the tracks to the second vertex you *need* to retrieve a copy
53  //this is safer and I would like to keep it...
54 
55  std::vector<VxTrackAtVertex*> tracksOfSecondVertex(second.getTracksAtVertex());
56 
57  const std::vector<VxTrackAtVertex*>& tracksOfFirstVertex(first.getTracksAtVertex());
58 
59  const std::vector<VxTrackAtVertex*>::const_iterator tracksBegin=tracksOfFirstVertex.begin();
60  const std::vector<VxTrackAtVertex*>::const_iterator tracksEnd=tracksOfFirstVertex.end();
61 
62  for (std::vector<VxTrackAtVertex*>::const_iterator tracksIter=tracksBegin;
63  tracksIter!=tracksEnd;
64  ++tracksIter) {
65 
66  if (*tracksIter!=nullptr) {
67 
68  tracksOfSecondVertex.push_back(*tracksIter);
69 
70  } else {
71  ATH_MSG_WARNING( " addTracksOfFirstVertexToSecondVertex method found an empty pointer for a track of the first vertex. Skipping it..." );
72  // continue;
73  }
74 
75  }
76 
77  second.setTracksAtVertex(std::move(tracksOfSecondVertex));
78 
79  }
80 
81  void JetFitterHelper::performKalmanConstraintToBePrimaryVertex(RecVertexPositions & positionToUpdate,
82  const VxVertexOnJetAxis & vertexToConstraint) const {
83 
84  //you need to copy the object... safeness here...
85  const Amg::VectorX position(positionToUpdate.position());
86  const Amg::MatrixX Ckm1(positionToUpdate.covariancePosition());
87 
88  if (position.rows()!=Ckm1.rows()) {
89  ATH_MSG_WARNING ("Error in performKalmanConstraintToBePrimaryVertex method:" <<
90  " covariance matrix of size " << Ckm1.rows() <<
91  " and vector of VertexPositions do not match. Not performing " <<
92  "any constraint... ");
93  return;
94  }
95 
96  int numbRow(numRow(vertexToConstraint.getNumVertex()));
97 
98  //create the linearization matrix for primary vertex constraint
99  Amg::MatrixX Gk(1,position.rows());Gk.setZero();
100  Gk(0,numbRow)=1;
101 
102  //applying *EXACT* constraint without error
103  Amg::MatrixX Rk_km1=Ckm1.similarity(Gk);//+ApparentError;
104 
105  if(Rk_km1.determinant()==0) {
106  ATH_MSG_WARNING ("The Rk_k-1 matrix inversion will fail - det(Rk)=0 detected."
107  " Constraint to primary vertex failed...");
108  return;
109  }
110  Amg::MatrixX Rk_km1_inverse=Rk_km1.inverse().eval();
111 
112  //kalman gain matrix
113  Amg::MatrixX Kk=Ckm1*Gk.transpose()*Rk_km1_inverse;
114 
115  //in this case the linearized constraint is just a scalar!
116  AmgVector(1) ConstraintInOldPosition;
117  ConstraintInOldPosition(0)=position[numbRow];
118 
119  Amg::VectorX newPosition=position-Kk*ConstraintInOldPosition;
120 
121  Amg::MatrixX newCovPosition=Ckm1-Kk*(2*Gk*Ckm1-Rk_km1*Kk.transpose());
122 
123  double newChi2=positionToUpdate.fitQuality().chiSquared()+Rk_km1_inverse.similarityT(ConstraintInOldPosition).determinant();
124 
125  double newNDF=positionToUpdate.fitQuality().numberDoF()+1;
126  //double newNDF=positionToUpdate.doubleNdf()+1.;//fitQuality().numberDoF()+1;
127 
128  const Amg::MatrixX& symmetrizedMatrix(newCovPosition);
129  positionToUpdate=Trk::RecVertexPositions(newPosition, symmetrizedMatrix, newNDF, newChi2);
130  }
131 
132  void JetFitterHelper::performKalmanConstraintToMergeVertices(RecVertexPositions & positionToUpdate,
133  const VxVertexOnJetAxis & vertexToConstraint1,
134  const VxVertexOnJetAxis & vertexToConstraint2) const {
135 
136  //you need to copy the object... safeness here...
137  const Amg::VectorX & position(positionToUpdate.position());
138  const Amg::MatrixX & Ckm1(positionToUpdate.covariancePosition());
139 
140  if (position.rows()!=Ckm1.rows()) {
141  ATH_MSG_WARNING("Error in performKalmanConstraintToBePrimaryVertex method: " <<
142  "covariance matrix dim " << Ckm1.rows() << " and vector of " <<
143  "VertexPositions dim " << position.rows() << " do not match."<<
144  " Not performing any constraint... ");
145  return;
146  }
147 
148  int numbRow1(numRow(vertexToConstraint1.getNumVertex()));
149  int numbRow2(numRow(vertexToConstraint2.getNumVertex()));
150 
151  //create the linearization matrix for primary vertex constraint
152  Amg::MatrixX Gk(1,position.rows()); Gk.setZero();
153  Gk(0,numbRow1)=1;
154  Gk(0,numbRow2)=-1;
155 
156  //applying *EXACT* constraint without error
157  AmgMatrix(1,1) Rk_km1=Ckm1.similarity(Gk);//+ApparentError;
158 
159  if(Rk_km1.determinant() ==0) {
160  ATH_MSG_WARNING ("The Rk_k-1 matrix inversion will fail - Rk_k=0 detected."
161  " Constraint to primary vertex failed...");
162  return;
163  }
164  AmgMatrix(1,1) Rk_km1_inverse=Rk_km1.inverse();
165 
166  //kalman gain matrix
167  Amg::MatrixX Kk=Ckm1*Gk.transpose()*Rk_km1_inverse;
168 
169  //in this case the linearized constraint is just a scalar!
170  AmgVector(1) ConstraintInOldPosition(1);
171  ConstraintInOldPosition(0)=position[numbRow1]-position[numbRow2];
172 
173  Amg::VectorX newPosition=position-Kk*ConstraintInOldPosition;
174 
175  Amg::MatrixX newCovPosition=Ckm1-Kk*(2*Gk*Ckm1-Rk_km1*Kk.transpose());
176 
177  double newChi2=positionToUpdate.fitQuality().chiSquared()+Rk_km1_inverse.similarityT(ConstraintInOldPosition).determinant();
178  double newNDF=positionToUpdate.fitQuality().numberDoF()+1;
179 
180  positionToUpdate=Trk::RecVertexPositions(newPosition,newCovPosition,newNDF,newChi2);
181  }
182 
183  VxVertexOnJetAxis & JetFitterHelper::mergeVerticesInJetCandidate(VxVertexOnJetAxis & vertex1,
184  VxVertexOnJetAxis & vertex2,
185  VxJetCandidate & myJetCandidate) const {
186 
187  if (&vertex1==myJetCandidate.getPrimaryVertex()) {
188  mergeVertexToPrimaryInJetCandidate(vertex2,myJetCandidate);
189  return *myJetCandidate.getPrimaryVertex();
190  }
191  if (&vertex2 == myJetCandidate.getPrimaryVertex()) {
192  mergeVertexToPrimaryInJetCandidate(vertex1, myJetCandidate);
193  return *myJetCandidate.getPrimaryVertex();
194  }
195 
196  addTracksOfFirstVertexToSecondVertex(vertex2,vertex1);
197 
198  //now you need to *delete* the second vertex in copyOfRecVertexPositions and
199  //in copyOfLinearizationPositions
200  const Amg::VectorX & recPosition=myJetCandidate.getRecVertexPositions().position();
201  Amg::VectorX linPosition=myJetCandidate.getLinearizationVertexPositions().position();
202  const Amg::VectorX & constraintPosition=myJetCandidate.getConstraintVertexPositions().position();
203  const Amg::MatrixX & covPosition=myJetCandidate.getRecVertexPositions().covariancePosition();
204  const Amg::MatrixX & covConstraintPosition=myJetCandidate.getConstraintVertexPositions().covariancePosition();
205 
206 
207  int numbVertex2=numRow(vertex2.getNumVertex());
208 
209  //set the linearized position of the merged vertex to the weighted mean of the previous positions
210  int numbVertex1=numRow(vertex1.getNumVertex());
211  if (covPosition(numbVertex1,numbVertex1)!=0&&covPosition(numbVertex2,numbVertex2)!=0) {
212  linPosition(numbVertex1)=(recPosition(numbVertex1)/covPosition(numbVertex1,numbVertex1)+
213  recPosition(numbVertex2)/covPosition(numbVertex2,numbVertex2))/
214  (1./covPosition(numbVertex1,numbVertex1)+1./covPosition(numbVertex2,numbVertex2));
215  } else {
216  ATH_MSG_WARNING ("one of the vertices to merge has error on position 0. should not happen, however: --> proceeding...");
217  }
218 
219  //call the function which deletes a single row (they're in the anonymous namespace, in Utilities.h)
220  Amg::VectorX reducedRecPositions=deleteRowFromVector(recPosition,numbVertex2);
221  Amg::VectorX reducedLinPositions=deleteRowFromVector(linPosition,numbVertex2);
222  Amg::VectorX reducedConstraintPositions=deleteRowFromVector(constraintPosition,numbVertex2);
223  Amg::MatrixX reducedCovPositions=deleteRowFromSymMatrix(covPosition,numbVertex2);
224  Amg::MatrixX reducedConstraintCovPositions=deleteRowFromSymMatrix(covConstraintPosition,numbVertex2);
225 
226  myJetCandidate.setRecVertexPositions(RecVertexPositions(reducedRecPositions,
227  reducedCovPositions,
228  0.,0.));
229  myJetCandidate.setConstraintVertexPositions(RecVertexPositions(reducedConstraintPositions,
230  reducedConstraintCovPositions,
231  0.,0.));
232 
233 
234  myJetCandidate.setLinearizationVertexPositions(VertexPositions(reducedLinPositions));
235 
236  //now as a last step you need to delete the vertex you're not using anymore...
237 
238  std::vector<VxVertexOnJetAxis*> copyOfVerticesAtJetCandidate=myJetCandidate.getVerticesOnJetAxis();
239 
240  const std::vector<VxVertexOnJetAxis*>::iterator VerticesBegin=copyOfVerticesAtJetCandidate.begin();
241  std::vector<VxVertexOnJetAxis*>::iterator VerticesEnd=copyOfVerticesAtJetCandidate.end();
242 
243  bool found=false;
244 
245  for (std::vector<VxVertexOnJetAxis*>::iterator VerticesIter=VerticesBegin;VerticesIter!=VerticesEnd;) {
246  if ((*VerticesIter)==&vertex2) {
247  delete *VerticesIter;
248  VerticesIter=copyOfVerticesAtJetCandidate.erase(VerticesIter);
249  VerticesEnd=copyOfVerticesAtJetCandidate.end();
250  found=true;
251  } else {
252  ++VerticesIter;
253  }
254  }
255 
256  if (!found) {
257  ATH_MSG_WARNING ("Could not find second vertex from which a cluster was created to delete it... Very strange... Check!!! ");
258  }
259 
260  //update myJetCandidate with the new vector of tracks
261  myJetCandidate.setVerticesOnJetAxis(copyOfVerticesAtJetCandidate);
262 
263  //now you should update the numbering scheme... CHECK CHECK... WHEN TO DO IT???????
264  return vertex1;//return the merged vertex
265 
266  }
267 
268  void JetFitterHelper::mergeVertexToPrimaryInJetCandidate(VxVertexOnJetAxis & vertex,VxJetCandidate & myJetCandidate) const {
269 
270  VxVertexOnJetAxis* primaryVertexPtr(myJetCandidate.getPrimaryVertex());
271 
272  if (primaryVertexPtr==nullptr) {
273  ATH_MSG_WARNING ("Pointer to the primary vertex is 0. No merging with primary vertex possible.");
274  return;
275  }
276 
277  VxVertexOnJetAxis primaryVertex(*primaryVertexPtr);
278 
279  addTracksOfFirstVertexToSecondVertex(vertex,primaryVertex);
280 
281  //now you need to *delete* the second vertex in copyOfRecVertexPositions and
282  //in copyOfLinearizationPositions
283  const Amg::VectorX & recPosition=myJetCandidate.getRecVertexPositions().position();
284  Amg::VectorX linPosition=myJetCandidate.getLinearizationVertexPositions().position();
285  const Amg::VectorX & constraintPosition=myJetCandidate.getConstraintVertexPositions().position();
286  const Amg::MatrixX & covPosition=myJetCandidate.getRecVertexPositions().covariancePosition();
287  const Amg::MatrixX & covConstraintPosition=myJetCandidate.getConstraintVertexPositions().covariancePosition();
288 
289  int numbVertex2=numRow(vertex.getNumVertex());
290 
291  //call the function which deletes a single row (they're in the anonymous namespace, in Utilities.h)
292  Amg::VectorX reducedRecPositions=deleteRowFromVector(recPosition,numbVertex2);
293  Amg::VectorX reducedLinPositions=deleteRowFromVector(linPosition,numbVertex2);
294  Amg::VectorX reducedConstraintPositions=deleteRowFromVector(constraintPosition,numbVertex2);
295  Amg::MatrixX reducedCovPositions=deleteRowFromSymMatrix(covPosition,numbVertex2);
296  Amg::MatrixX reducedConstraintCovPositions=deleteRowFromSymMatrix(covConstraintPosition,numbVertex2);
297 
298  myJetCandidate.setRecVertexPositions(RecVertexPositions(reducedRecPositions,
299  reducedCovPositions,
300  0.,0.));
301  myJetCandidate.setConstraintVertexPositions(RecVertexPositions(reducedConstraintPositions,
302  reducedConstraintCovPositions,
303  0.,0.));
304 
305  myJetCandidate.setLinearizationVertexPositions(VertexPositions(reducedLinPositions));
306 
307  //now as a last step you need to delete the vertex you're not using anymore...
308 
309  std::vector<VxVertexOnJetAxis*> copyOfVerticesAtJetCandidate=myJetCandidate.getVerticesOnJetAxis();
310 
311  const std::vector<VxVertexOnJetAxis*>::iterator VerticesBegin=copyOfVerticesAtJetCandidate.begin();
312  std::vector<VxVertexOnJetAxis*>::iterator VerticesEnd=copyOfVerticesAtJetCandidate.end();
313 
314  bool found=false;
315 
316  for (std::vector<VxVertexOnJetAxis*>::iterator VerticesIter=VerticesBegin;VerticesIter!=VerticesEnd;) {
317  if ((*VerticesIter)==&vertex) {
318  delete *VerticesIter;
319  VerticesIter=copyOfVerticesAtJetCandidate.erase(VerticesIter);
320  VerticesEnd=copyOfVerticesAtJetCandidate.end();
321  found=true;
322  } else {
323  ++VerticesIter;
324  }
325  }
326 
327  if (!found) {
328  ATH_MSG_WARNING ("Could not find second vertex from which a cluster with the primary vertex was created to delete it... Very strange... Check!!! ");
329  }
330 
331  //update myJetCandidate with the new vector of tracks
332  myJetCandidate.setVerticesOnJetAxis(copyOfVerticesAtJetCandidate);
333  myJetCandidate.setPrimaryVertex(&primaryVertex);
334 
335  //now you should update the numbering scheme... CHECK CHECK... WHEN TO DO IT???????
336  }
337 
338 }
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
RecVertex.h
python.SystemOfUnits.second
int second
Definition: SystemOfUnits.py:120
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition: EventPrimitives.h:32
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:29
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
JetVtxParamDefs.h
VxVertexOnJetAxis.h
VxJetCandidate.h
Trk::VxVertexOnJetAxis
VxVertexOnJetAxis inherits from Vertex.
Definition: VxVertexOnJetAxis.h:79
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
const
bool const RAWDATA *ch2 const
Definition: LArRodBlockPhysicsV0.cxx:562
AmgMatrix
#define AmgMatrix(rows, cols)
Definition: EventPrimitives.h:51
JetFitterHelper.h
similarity
Matrix< Scalar, OtherDerived::RowsAtCompileTime, OtherDerived::RowsAtCompileTime > similarity(const MatrixBase< OtherDerived > &m) const
Definition: AmgMatrixBasePlugin.h:130
Trk::JetFitterHelper
Februar 2007 (c) Atlas Detector Reconstruction Software.
Definition: JetFitterHelper.h:44
Trk::JetFitterHelper::JetFitterHelper
JetFitterHelper(const std::string &t, const std::string &n, const IInterface *p)
Constructor.
Definition: JetFitterHelper.cxx:52
beamspotman.n
n
Definition: beamspotman.py:731
AmgVector
AmgVector(4) T2BSTrackFilterTool
Definition: T2BSTrackFilterTool.cxx:114
VxTrackAtVertex.h
Utilities.h
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Amg
Definition of ATLAS Math & Geometry primitives (Amg)
Definition: AmgStringHelpers.h:19
Trk::VxJetCandidate
Definition: VxJetCandidate.h:72
Trk::RecVertexPositions
Definition: RecVertexPositions.h:34
IDTPM::chiSquared
float chiSquared(const U &p)
Definition: TrackParametersHelper.h:136
Trk::GsfMeasurementUpdator::fitQuality
FitQualityOnSurface fitQuality(const MultiComponentState &, const MeasurementBase &)
Method for determining the chi2 of the multi-component state and the number of degrees of freedom.
Definition: GsfMeasurementUpdator.cxx:845
similarityT
Matrix< Scalar, OtherDerived::RowsAtCompileTime, OtherDerived::RowsAtCompileTime > similarityT(const MatrixBase< OtherDerived > &m) const
similarityT method : yields ms = m^T*s*m
Definition: AmgMatrixBasePlugin.h:138
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
CondAlgsOpts.found
int found
Definition: CondAlgsOpts.py:101
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
DeMoScan.first
bool first
Definition: DeMoScan.py:534
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:569
Trk::VertexPositions
VertexPositions class to represent and store a vertex.
Definition: VertexPositions.h:36
AthAlgTool
Definition: AthAlgTool.h:26
RecVertexPositions.h
Generated on Sun Jun 30 2024 21:18:48 for ATLAS Offline Software by   doxygen 1.8.18