ATLAS Offline Software
JetFitterRoutines.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /***************************************************************************
6  JetFitterRoutines.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 <iostream>
23 
24 // to get AmgMatrix plugin:
25 #include "GeoPrimitives/GeoPrimitives.h"
26 
27 #include "TrkJetVxFitter/JetFitterRoutines.h"
28 #include "VxVertex/RecVertex.h"
29 #include "VxJetVertex/VxJetCandidate.h"
30 #include "VxVertex/VxTrackAtVertex.h"
31 #include "VxJetVertex/VxVertexOnJetAxis.h"
32 #include "VxJetVertex/RecVertexPositions.h"
33 #include "VxJetVertex/VxClusteringTable.h"
34 #include "VxJetVertex/JetVtxParamDefs.h"
35 #include "VxJetVertex/PairOfVxVertexOnJetAxis.h"
36 #include "TrkJetVxFitter/KalmanVertexOnJetAxisUpdator.h"
37 #include "TrkJetVxFitter/KalmanVertexOnJetAxisSmoother.h"
38 #include "TrkJetVxFitter/JetFitterHelper.h"
39 #include "TrkJetVxFitter/JetFitterInitializationHelper.h"
40 #include "VxJetVertex/VxJetFitterDebugInfo.h"
41 
42 #include "TrkNeutralParameters/NeutralParameters.h"
43 #include "TrkJetVxFitter/Utilities.h"
44 
45 #include <TMath.h>
46 #include <cmath>
47 #include <sstream>
48 
49 #include "TrkJetVxFitter/TrkDistanceFinderNeutralCharged.h"
50 #include "TrkJetVxFitter/TrkDistanceFinderNeutralNeutral.h"
51 
52 
53 //#define JetFitterRoutines_DEBUG2
54 //#define JetFitterRoutines_DEBUG
55 
56 namespace Trk
57 {
58 
59  JetFitterRoutines::JetFitterRoutines(const std::string& t, const std::string& n, const IInterface* p) :
60  AthAlgTool(t,n,p),
61  m_initializationHelper("Trk::JetFitterInitializationHelper", this),
62  m_helper("Trk::JetFitterHelper", this),
63  m_updator("Trk::KalmanVertexOnJetAxisUpdator", this),
64  m_smoother("Trk::KalmanVertexOnJetAxisSmoother", this),
65  m_minDistanceFinder("Trk::TrkDistanceFinderNeutralCharged/TrkDistanceFinderNeutralCharged", this),
66  m_minDistanceFinderNeutral("Trk::TrkDistanceFinderNeutralNeutral/TrkDistanceFinderNeutralNeutral", this),
67  m_fast(false),
68  m_maxDRshift(.5),
69  m_noPrimaryVertexRefit(false),
70  m_maxR(1150.),//max R of ID
71  m_maxZ(2727.)//max Z of ID
72  {
73  declareProperty("KalmanVertexOnJetAxisUpdator",m_updator);
74  declareProperty("KalmanVertexOnJetAxisSmoother",m_smoother);
75  declareProperty("JetFitterHelper",m_helper);
76  declareProperty("JetFitterInitializationHelper",m_initializationHelper);
77  declareProperty("JetFitterMinimumDistanceFinder",m_minDistanceFinder);
78  declareProperty("JetFitterMinimumDistanceFinderNeutral",m_minDistanceFinderNeutral);
79  declareProperty("ID_maxR",m_maxR);
80  declareProperty("ID_maxZ",m_maxZ);
81  declareInterface< JetFitterRoutines >(this) ;
82  declareProperty("BeFast",m_fast);
83  declareProperty("maxDRshift",m_maxDRshift);
84  declareProperty("noPrimaryVertexRefit",m_noPrimaryVertexRefit);
85 
86  }
87 
88 
89 
90  JetFitterRoutines::~JetFitterRoutines() = default;
91 
92 
93  StatusCode JetFitterRoutines::initialize() {
94 
95  AthAlgTool::initialize().ignore();
96  ATH_CHECK( m_fieldCacheCondObjInputKey.initialize() );
97 
98  //retrieving the udator itself
99  ATH_CHECK( m_helper.retrieve() );
100 
101  ATH_CHECK( m_initializationHelper.retrieve() );
102 
103  ATH_CHECK( m_minDistanceFinder.retrieve() );
104 
105  ATH_CHECK( m_minDistanceFinderNeutral.retrieve() );
106 
107  ATH_CHECK( m_updator.retrieve() );
108 
109  ATH_CHECK( m_smoother.retrieve() );
110 
111  return StatusCode::SUCCESS;
112  }
113 
114 
115  void JetFitterRoutines::initializeToMinDistancesToJetAxis(VxJetCandidate* myJetCandidate) const {
116 
117  ATH_MSG_DEBUG ("initializingToMinDistancesToJetAxis is now implemented! Will converge faster!!! Neutrals are fully supported...");
118 
119  VertexPositions & linVertexPositions=myJetCandidate->getLinearizationVertexPositions();
120  Amg::VectorX linPositions=linVertexPositions.position();
121 
122 
123  const Trk::RecVertexPositions & recVtxPosition=myJetCandidate->getRecVertexPositions();
124  const Amg::VectorX & recPosition=recVtxPosition.position();
125 
126  Amg::Vector3D a(recPosition[Trk::jet_xv],
127  recPosition[Trk::jet_yv],
128  recPosition[Trk::jet_zv]);
129 
130  auto const sinRecJetTheta = std::sin(recPosition[Trk::jet_theta]);
131  auto const sinRecJetPhi = std::sin(recPosition[Trk::jet_phi]);
132  auto const cosRecJetTheta = std::cos(recPosition[Trk::jet_theta]);
133  auto const cosRecJetPhi = std::cos(recPosition[Trk::jet_phi]);
134 
135  auto const absRecJetTheta = std::abs(recPosition[Trk::jet_theta]);
136  auto const abssinRecJetTheta = std::abs(sinRecJetTheta);
137  auto const abscosRecJetTheta = std::abs(cosRecJetTheta);
138 
139  Amg::Vector3D b(cosRecJetPhi*sinRecJetTheta,
140  sinRecJetPhi*sinRecJetTheta,
141  cosRecJetTheta);
142 
143  NeutralTrack myJetAxis(a,b);
144 
145  //Yes, but the seeding is more than just speed!
146  const std::vector<VxVertexOnJetAxis*> & associatedVertices=myJetCandidate->getVerticesOnJetAxis();
147 
148  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxBegin=associatedVertices.begin();
149  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxEnd=associatedVertices.end();
150 
151  if (associatedVertices.empty()) {//Was that your intention? to be checked... 15.03.2007
152  SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCacheCondObjInputKey};
153  const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
154  if (!readHandle.isValid()) {
155  std::stringstream msg;
156  msg << "Failed to retrieve magmnetic field conditions data " << m_fieldCacheCondObjInputKey.key() << ".";
157  throw std::runtime_error(msg.str());
158  }
159  MagField::AtlasFieldCache fieldCache;
160  fieldCondObj->getInitializedCache (fieldCache);
161  for (std::vector<VxVertexOnJetAxis*>::const_iterator VtxIter=VtxBegin;VtxIter!=VtxEnd;++VtxIter) {
162  VxVertexOnJetAxis* myVertex=(*VtxIter);
163  if (myVertex!=nullptr) {
164 
165  const std::vector<VxTrackAtVertex*> & tracksAtVertex=myVertex->getTracksAtVertex();
166  if (tracksAtVertex.size()>1) {
167  ATH_MSG_DEBUG( "Warning in JetFitterInitializationHelper.Number of tracks at vertex is bigger than one, "
168  << "even during initialization phase. Skipping this vertex (already initialized)..." );
169  }
170  else if (tracksAtVertex.empty())
171  {
172  ATH_MSG_WARNING( "No track at vertex. Internal fitter error. Contact author (GP) ... " );
173  }
174  else
175  {
176  ATH_MSG_VERBOSE( " TrackAtVertexSize is: " << tracksAtVertex.size() );
177  ATH_MSG_VERBOSE( " pointer to first element: " << tracksAtVertex[0] );
178  ATH_MSG_VERBOSE( " pointer to initialPerigee: " << tracksAtVertex[0]->initialPerigee() );
179 
180  // RS 19.04.2011 try to fix coverity defect 22750
181  const Trk::Perigee* ptr = dynamic_cast<const Trk::Perigee*>((tracksAtVertex[0]->initialPerigee()));
182  if (ptr)
183  {
184 
185  double distOnAxis=-999.;
186  std::pair<Amg::Vector3D,double> result;
187  try {
188  result=m_minDistanceFinder->getPointAndDistance(myJetAxis,*ptr,distOnAxis, fieldCache);
189 
190  double R=distOnAxis*sinRecJetTheta;
191  double Z=distOnAxis*cosRecJetTheta;
192 
193  if (std::abs(R)>m_maxR)
194  {
195 
196  if (absRecJetTheta>1e-8)
197  {
198  ATH_MSG_DEBUG( " Closest distance of track to jet axis is outside ID envelope, R=" << R << ", setting to R= " << m_maxR );
199  distOnAxis=m_maxR /abssinRecJetTheta;
200  }
201  }
202 
203  Z=distOnAxis*cosRecJetTheta;
204  if (std::abs(Z)>m_maxZ)
205  {
206  if (abscosRecJetTheta>1e-8)
207  {
208  ATH_MSG_DEBUG( " Closest distance of track to jet axis is outside ID envelope, Z=" << Z << ", setting to Z= " << m_maxZ );
209  distOnAxis=m_maxZ / cosRecJetTheta;
210  }
211  }
212 
213  linPositions[numRow(myVertex->getNumVertex())]=distOnAxis;
214 
215 
216  } catch (Error::NewtonProblem e) {
217  ATH_MSG_WARNING( "Problem with Newton finder " << e.p );
218  } catch (...) {
219  ATH_MSG_ERROR( "Could not catch error " );
220  }
221  ATH_MSG_DEBUG ("initializingToMinDistancesToJetAxis for vertex number... " <<
222  myVertex->getNumVertex() << " to distance " << distOnAxis <<
223  " distance to axis " << result.second);
224  }
225  // FIXME THIS PART IS DEAD IN NEW TRACKING EDM:
226  else if (dynamic_cast<const Trk::NeutralPerigee*>((tracksAtVertex[0]->initialPerigee()))!=nullptr)
227  {
228  double distOnAxis=-999.;
229  std::pair<Amg::Vector3D,double> result;
230 
231  const Trk::NeutralPerigee * neutralTrack =
232  dynamic_cast<const Trk::NeutralPerigee*>((tracksAtVertex[0]->initialPerigee()));
233 
234  NeutralTrack myNeutralTrack(neutralTrack->position(),neutralTrack->momentum());
235 
236  result=m_minDistanceFinderNeutral->getPointAndDistance(myJetAxis,myNeutralTrack,distOnAxis);
237 
238  double R=distOnAxis*sinRecJetTheta;
239  double Z=distOnAxis*cosRecJetTheta;
240 
241  if (std::abs(R)>m_maxR)
242  {
243 
244  if (absRecJetTheta>1e-8)
245  {
246  ATH_MSG_DEBUG( " Closest distance of track to jet axis is outside ID envelope, R=" << R << ", setting to R= " << m_maxR );
247  distOnAxis=m_maxR / abssinRecJetTheta;
248  }
249  }
250 
251  Z=distOnAxis*cosRecJetTheta;
252  if (std::abs(Z)>m_maxZ)
253  {
254  if (abscosRecJetTheta>1e-8)
255  {
256  ATH_MSG_DEBUG( " Closest distance of track to jet axis is outside ID envelope, Z=" << Z << ", setting to Z= " << m_maxZ );
257  distOnAxis=m_maxZ / cosRecJetTheta;
258  }
259  }
260 
261  linPositions[numRow(myVertex->getNumVertex())]=distOnAxis;
262  ATH_MSG_DEBUG( "initializingToMinDistancesToJetAxis for vertex from NEUTRAL number... " <<
263  myVertex->getNumVertex() << " to distance " <<
264  distOnAxis << " distance to axis " << result.second );
265 
266  }
267  else
268  {
269  ATH_MSG_WARNING( "Could not cast to neither CHARGED or NEUTRAL! This error is not FATAL" );
270  }
271  }
272  } else {
273  ATH_MSG_WARNING( "Warning in JetFitterInitializationHelper.Inconsistency found. Pointer to VxVertexOnJetAxis should be different from zero. Skipping track..." );
274  throw ("Warning in JetFitterInitializationHelper.Inconsistency found. Pointer to VxVertexOnJetAxis should be different from zero. Skipping track...");
275  }
276  }
277 
278  linVertexPositions.setPosition(linPositions);
279 
280  } else {
281  ATH_MSG_DEBUG ("No Associated Vertices found! no initialization to minimum distance is possible.");
282  }
283  }
284 
285  void JetFitterRoutines::performTheFit(VxJetCandidate* myJetCandidate,
286  int num_maxiterations,//default 20
287  bool treat_sign_flip,//default true
288  int num_signflip_maxiterations,//default 10
289  double deltachi2_convergence) const {//default 0.001
290 
291  bool isInitialized=this->checkJetCandidate(myJetCandidate);
292 
293  if (!isInitialized) {
294  ATH_MSG_DEBUG ("JetFitter found the candidate was not correctly initialized. Not proceeding with the fitt...");
295  return;
296  }
297 
298 
299 
300  //linearization is DONE. now do fit with all the vertices (updating them through the knowledge coming from all the tracks parameters)
301 
302  int num_iteration_signflip=0;
303  double lastchi2=-99.;
304  bool converged=false;
305 
306  if (treat_sign_flip) {
307  do {
308 
309  myJetCandidate->setRecVertexPositions(RecVertexPositions(myJetCandidate->getConstraintVertexPositions()));
310 
311  ATH_MSG_DEBUG( myJetCandidate->getRecVertexPositions() );
312 
313  //linearization position is now stored in a different VertexPositions attached to the VxJetCandidate...
314  m_initializationHelper->linearizeAllTracks(myJetCandidate,true);
315 
316  //now update all vertices with the new information
317  updateAllVertices(myJetCandidate);
318 
319  //now do the chi2 update of all tracks after end of fit iterations
320  //(this permits to invert the big weight matrix only once at the end!!)
321  // updateChi2AllVertices(myJetCandidate);
322 
323  num_iteration_signflip+=1;
324 
325  const RecVertexPositions & myRecPosition=myJetCandidate->getRecVertexPositions();
326 
327  const FitQuality & myFitQuality=myRecPosition.fitQuality();
328 
329  double actualchi2=myFitQuality.chiSquared();
330 
331  auto const absActLastChi2 = std::abs(actualchi2-lastchi2);
332 
333 #ifdef JetFitterRoutines_DEBUG
334  ATH_MSG_DEBUG( " last chi2 " << lastchi2 << " actual chi2 " << actualchi2 << " difference " <<
335  absActLastChi2<< " < " << deltachi2_convergence << " ? " << " ndf " << myFitQuality.numberDoF() );
336 #endif
337 
338  if (absActLastChi2<deltachi2_convergence) {
339  converged=true;
340  } else {
341 
342  //GP 12.3.2012 Check that no vertex is outside the ID acceptance
343  copyRecoPositionsToLinearizationPositions(*myJetCandidate);
344 
345  }
346 
347  lastchi2=actualchi2;
348 
349  } while ((!converged)&&num_iteration_signflip<num_signflip_maxiterations);
350 
351 
352  if (converged) {
353  ATH_MSG_VERBOSE( " For sign flip treatment there was convergence after " << num_iteration_signflip );
354  }
355  }
356 
357  // if (msgSvc()->outputLevel()== MSG::DEBUG||msgSvc()->outputLevel()== MSG::VERBOSE) {
358 #ifdef JetFitterRoutines_DEBUG
359  ATH_MSG_DEBUG( "JetFitterRoutines: after convergence with sign flip treatment: " << myJetCandidate->getRecVertexPositions() );
360 #endif
361  // }
362 
363  int num_iteration=0;
364  lastchi2=-99.;
365  converged=false;
366 
367  do {
368 
369 
370  myJetCandidate->setRecVertexPositions(RecVertexPositions(myJetCandidate->getConstraintVertexPositions()));
371 
372 
373  //linearization position is now stored in a different VertexPositions attached to the VxJetCandidate...
374  m_initializationHelper->linearizeAllTracks(myJetCandidate,false);
375 
376  //now update all vertices with the new information
377  updateAllVertices(myJetCandidate);
378 
379  //now do the chi2 update of all tracks after end of fit iterations
380  //(this permits to invert the big weight matrix only once at the end!!)
381  // updateChi2AllVertices(myJetCandidate);
382 
383  num_iteration+=1;
384 
385  const RecVertexPositions & myRecPosition=myJetCandidate->getRecVertexPositions();
386 
387  ATH_MSG_VERBOSE( " num_iteration (full fit): " << num_iteration << " det " << myRecPosition.covariancePosition().determinant() << " recVertex " << myJetCandidate->getRecVertexPositions() );
388 
389  const FitQuality & myFitQuality=myRecPosition.fitQuality();
390  double actualchi2=myFitQuality.chiSquared();
391 
392  auto const absActLastChi2 = std::abs(actualchi2-lastchi2);
393 
394  // if (msgSvc()->outputLevel()== MSG::VERBOSE) {
395 #ifdef JetFitterRoutines_DEBUG
396  ATH_MSG_DEBUG( " without sign flip: last chi2 " << lastchi2 << " actual chi2 " << actualchi2 << " difference " <<
397  absActLastChi2 << " < " << deltachi2_convergence << " ? " << " ndf " << myFitQuality.numberDoF() );
398 #endif
399  // }
400 
401  if (absActLastChi2<deltachi2_convergence) {
402  converged=true;
403  } else {
404  //now set the linearization position for the next step to the actual fitted vertex
405  //GP 12.3.2012 Check that no vertex is outside the ID acceptance
406  copyRecoPositionsToLinearizationPositions(*myJetCandidate);
407  //OLD myJetCandidate->setLinearizationVertexPositions(myJetCandidate->getRecVertexPositions());
408  }
409  lastchi2=actualchi2;
410 
411  } while ((!converged)&&num_iteration<num_maxiterations);
412 
413  if (converged) {
414  ATH_MSG_VERBOSE( " Fit without sign flip treatment there was convergence after " << num_iteration );
415  }
416 
417  // if (msgSvc()->outputLevel()== MSG::VERBOSE) {
418  #ifdef JetFitterRoutines_DEBUG
419  ATH_MSG_DEBUG( "JetFitterRoutines: after convergence without sign flip treatment: " << myJetCandidate->getRecVertexPositions() );
420  #endif
421  // }
422  if (num_iteration>=num_maxiterations)
423  {
424  ATH_MSG_DEBUG( "There wasn't convergence in JetFitter after: " << num_maxiterations );
425  }
426 
427  //now only the smoothing is missing as a last step... (updated momenta, chi2 + ndf of clusters,...)
428 
429  smoothAllVertices(myJetCandidate);
430 
431  Trk::VxJetFitterDebugInfo * & myDebugInfo=myJetCandidate->getDebugInfo();
432  if (myDebugInfo!=nullptr) {
433  delete myDebugInfo;
434  }
435  myDebugInfo=new VxJetFitterDebugInfo();
436  myDebugInfo->setNumFitIterations(num_iteration);
437  myDebugInfo->setSignFlipNumFitIterations(num_iteration_signflip);
438 
439  }
440 
441 
442  void JetFitterRoutines::updateAllVertices(VxJetCandidate* myJetCandidate) const {
443 
444 // ATH_MSG_DEBUG( " Updating PV " );
445 
446  int n_iteration=0;
447 
448 
449  if (!m_noPrimaryVertexRefit) {
450  //new iteration
451  VxVertexOnJetAxis* myPrimary=myJetCandidate->getPrimaryVertex();
452  const std::vector<VxTrackAtVertex*> & primaryVectorTracks=myPrimary->getTracksAtVertex();
453 
454  const std::vector<VxTrackAtVertex*>::const_iterator primaryVectorTracksBegin=primaryVectorTracks.begin();
455  const std::vector<VxTrackAtVertex*>::const_iterator primaryVectorTracksEnd=primaryVectorTracks.end();
456 
457  for (std::vector<VxTrackAtVertex*>::const_iterator primaryVectorIter=primaryVectorTracksBegin;
458  primaryVectorIter!=primaryVectorTracksEnd;++primaryVectorIter) {
459 
460  ATH_MSG_VERBOSE( " RecVertexPositions before update " << myJetCandidate->getRecVertexPositions() );
461 
462  if ((!m_fast)) {
463  m_updator->add(*primaryVectorIter,myPrimary,myJetCandidate);
464  } else {
465  m_updator->addWithFastUpdate(*primaryVectorIter,myPrimary,myJetCandidate);
466  }
467 
468  const RecVertexPositions & myRecPosition=myJetCandidate->getRecVertexPositions();
469 
470  ATH_MSG_VERBOSE( " Determinant after PRIMARY VTX update: " << n_iteration << " det: " << myRecPosition.covariancePosition().determinant() << " recVertex " << myJetCandidate->getRecVertexPositions() );
471  }
472  }
473 
474  n_iteration=0;
475 
476  // const RecVertexPositions & myRecPositionBeg=myJetCandidate->getRecVertexPositions();
477 
478  const std::vector<VxVertexOnJetAxis*> & associatedVertices=myJetCandidate->getVerticesOnJetAxis();
479 
480  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxBegin=associatedVertices.begin();
481  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxEnd=associatedVertices.end();
482 
483  for (std::vector<VxVertexOnJetAxis*>::const_iterator VtxIter=VtxBegin;VtxIter!=VtxEnd;++VtxIter) {
484 
485 // ATH_MSG_DEBUG( " Updating an SV along jet axis " );
486 
487  const std::vector<VxTrackAtVertex*> & tracksAtVertex=(*VtxIter)->getTracksAtVertex();
488 
489  const std::vector<VxTrackAtVertex*>::const_iterator TracksBegin=tracksAtVertex.begin();
490  const std::vector<VxTrackAtVertex*>::const_iterator TracksEnd=tracksAtVertex.end();
491 
492  for (std::vector<VxTrackAtVertex*>::const_iterator TrackVectorIter=TracksBegin;
493  TrackVectorIter!=TracksEnd;++TrackVectorIter) {
494 
495  ATH_MSG_VERBOSE( " RecVertexPositions before update " << myJetCandidate->getRecVertexPositions() );
496 
497  if (!m_fast) {
498  m_updator->add(*TrackVectorIter,*VtxIter,myJetCandidate);
499  } else {
500  m_updator->addWithFastUpdate(*TrackVectorIter,*VtxIter,myJetCandidate);
501  }
502 
503  n_iteration+=1;
504 
505  const RecVertexPositions & myRecPosition=myJetCandidate->getRecVertexPositions();
506 
507  ATH_MSG_VERBOSE( " Determinant after sec update: " << n_iteration << " det: " << myRecPosition.covariancePosition().determinant() << " recVertex " << myJetCandidate->getRecVertexPositions() );
508  }
509  }
510 
511  myJetCandidate->getRecVertexPositions().finalizePosition();
512  }
513 
514  void JetFitterRoutines::smoothAllVertices(VxJetCandidate* myJetCandidate) const {
515  //new iteration
516  VxVertexOnJetAxis* myPrimary=myJetCandidate->getPrimaryVertex();
517 
518  if (!m_fast) {
519  m_smoother->update(myPrimary,myJetCandidate);
520  } else {
521  m_smoother->fastUpdate(myPrimary,myJetCandidate);
522  }
523 
524  const std::vector<VxVertexOnJetAxis*> & associatedVertices=myJetCandidate->getVerticesOnJetAxis();
525 
526  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxBegin=associatedVertices.begin();
527  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxEnd=associatedVertices.end();
528 
529  for (std::vector<VxVertexOnJetAxis*>::const_iterator VtxIter=VtxBegin;VtxIter!=VtxEnd;++VtxIter) {
530  if (!m_fast) {
531  m_smoother->update(*VtxIter,myJetCandidate);
532  } else {
533  m_smoother->fastUpdate(*VtxIter,myJetCandidate);
534  }
535  }
536  }
537 
538 
539 
540 
541  bool JetFitterRoutines::checkJetCandidate(VxJetCandidate* myJetCandidate) const {
542 
543  int sizeprimary(0);
544 
545  //check if primaryvertex exists
546  const VxVertexOnJetAxis* myPrimary=myJetCandidate->getPrimaryVertex();
547  if (myPrimary==nullptr) {
548  ATH_MSG_WARNING( "No primary vertex found in VxJetCandidate class. Initialization was not done correctly..." );
549  return false;
550  }
551 
552  bool ok = true;
553  if (myPrimary->getNumVertex() != -10) {
554  ok = false;
555  ATH_MSG_WARNING("Numvertex of primary vertex not correctly initialized. "
556  "Not proceeding with the fit!");
557  }
558 
559  const std::vector<VxTrackAtVertex*>& primaryVectorTracks =
560  myPrimary->getTracksAtVertex();
561 
562  sizeprimary = primaryVectorTracks.size();
563 
564  ok = (std::find(primaryVectorTracks.begin(),
565  primaryVectorTracks.end(),
566  nullptr) == primaryVectorTracks.end());
567  if (not ok)
568  ATH_MSG_WARNING("One of the VxTrackAtVertex is a null pointer. Not "
569  "proceeding with the fit!");
570 
571  if (!ok) {
572  return false;
573  }
574  // end if else rimary==0
575 
576  // check std::vector<VxVertexOnJetAxis*> (if pointers are not empty and if
577  // all associated tracks are not empty)
578  const std::vector<VxVertexOnJetAxis*>& tracksOfVertex =
579  myJetCandidate->getVerticesOnJetAxis();
580 
581  auto badVertex = [](VxVertexOnJetAxis* pVertex) {
582  return (pVertex == nullptr) or (pVertex->getNumVertex() < 0);
583  };
584  ok =
585  (std::find_if(tracksOfVertex.begin(), tracksOfVertex.end(), badVertex) ==
586  tracksOfVertex.end());
587  if (not ok) {
588  ATH_MSG_WARNING(
589  "One of the VxTrackAtVertex is a null pointer or uninitialized. Not "
590  "proceeding with the fit!"); // Two error messages combined into one
591  }
592  if (not ok) {
593  return false;
594  }
595 
596  //now check if there is some track at least to do the fit...
597 
598  if (tracksOfVertex.empty()&&sizeprimary==0) {
599  ATH_MSG_DEBUG( "No tracks at primary, no tracks on jet axis. Not proceeding with the fit!" );
600  return false;
601  }
602 
603  //now check if the number of tracks corrisponds to the number of components of the recVertexPositions
604  const Trk::RecVertexPositions & myRecVertexPositions=myJetCandidate->getRecVertexPositions();
605 
606  const Amg::VectorX& myPosition=myRecVertexPositions.position();
607 
608  if (static_cast<unsigned int>(tracksOfVertex.size()+5)!=static_cast<unsigned int>(myPosition.rows())) {
609  ATH_MSG_WARNING ( "The position matrix has " << myPosition.rows()
610  << " components while " << tracksOfVertex.size()+5
611  << " are expected. Not proceeding with the fit " );
612  }
613 
614  const Amg::MatrixX & myErrorMatrix=myRecVertexPositions.covariancePosition();
615 
616  //check if symMatrix and myPosition are compatible (size matches)
617  if (myPosition.rows()!=myErrorMatrix.rows()) {
618  ATH_MSG_WARNING ("The dimension of the position vector and the covariance matrix does not match. Not performing fit...");
619  return false;
620  }
621 
622  //check if all the diagonal values of the covariance matrix are not zero
623  for (int i=0;i<myPosition.rows();i++) {
624  if (std::abs(myErrorMatrix(i,i))<1e-20) {
625  ATH_MSG_WARNING ("Value of cov matrix component n. " << i << " has a value smaller than 1e-8. Not considered as possible. Not performing fit...");
626  return false;
627  }
628  }
629 
630  //well, more checks in the future!!!
631 
632  return true;
633 
634  }//end method
635 
636 
637  std::pair<double,bool> JetFitterRoutines::fastProbabilityOfMerging(const VxVertexOnJetAxis* firstVertex,
638  const VxVertexOnJetAxis* secondVertex,
639  const VxJetCandidate* myJetCandidate) const {
640  //this method is for evaluating the probability of merging in a very fast (and rough) way
641  //if above a threshold you should think about evaluating the "fullProbabilityOfMerging" :-)
642 
643  const VxVertexOnJetAxis * PrimaryVertex=myJetCandidate->getPrimaryVertex();
644 
645  if (firstVertex==PrimaryVertex) {
646  return fastProbabilityOfMergingWithPrimary(secondVertex,myJetCandidate);
647  }
648  if (secondVertex==PrimaryVertex) {
649  return fastProbabilityOfMergingWithPrimary(firstVertex,myJetCandidate);
650  }
651 
652  return fastProbabilityOfMergingNoPrimary(firstVertex,secondVertex,myJetCandidate);
653 
654  }
655 
656 
657  std::pair<double,bool> JetFitterRoutines::fastProbabilityOfMergingWithPrimary(const VxVertexOnJetAxis* otherVertex,
658  const VxJetCandidate* myJetCandidate) const {
659 
660 
661  //first get a copy of all vertex positions (this you can't avoid I fear...)
662  RecVertexPositions copyOfRecVertexPositions(myJetCandidate->getRecVertexPositions());
663 
664 
665  double oldchi2=copyOfRecVertexPositions.fitQuality().chiSquared();
666  double oldndf=copyOfRecVertexPositions.fitQuality().numberDoF();
667 
668  //#ifdef JetFitterRoutines_DEBUG2
669  const Amg::VectorX & positionVector=copyOfRecVertexPositions.position();
670  const Amg::MatrixX & positionCov=copyOfRecVertexPositions.covariancePosition();
671  double phiold=positionVector(Trk::jet_phi);
672  double thetaold=positionVector(Trk::jet_theta);
673  double phierr=std::sqrt(positionCov(Trk::jet_phi,Trk::jet_phi));
674  double thetaerr=std::sqrt(positionCov(Trk::jet_theta,Trk::jet_theta));
675 
676  //now do the merging of the second cluster to the primary vertex...
677  m_helper->performKalmanConstraintToBePrimaryVertex(copyOfRecVertexPositions,
678  *otherVertex);
679 
680  double phinew=positionVector[Trk::jet_phi];
681  double thetanew=positionVector[Trk::jet_theta];
682  bool isshifted=
683  std::pow((phinew-phiold)/phierr,2)+std::pow((thetanew-thetaold)/thetaerr,2)>m_maxDRshift*m_maxDRshift;
684 
685  double chi2 = myJetCandidate->getPrimaryVertex()->fitQuality().chiSquared() +
686  otherVertex->fitQuality().chiSquared() +
687  copyOfRecVertexPositions.fitQuality().chiSquared()-oldchi2;
688 
689  double ndf = copyOfRecVertexPositions.fitQuality().numberDoF()-oldndf +
690  myJetCandidate->getPrimaryVertex()->fitQuality().numberDoF() +
691  otherVertex->fitQuality().numberDoF();
692 
693  if (chi2<0||ndf<0) {
694  ATH_MSG_WARNING (" In the compatibility estimation chi2: " << chi2 << " ndf " << ndf << " giving back 0 prob ");
695  return std::pair<double,bool>(0,isshifted);
696  }
697 
698  return std::pair<double,bool>(TMath::Prob(chi2,(int)std::floor(ndf+0.5)),isshifted);
699 
700  }
701 
702  std::pair<double,bool> JetFitterRoutines::fastProbabilityOfMergingNoPrimary(const VxVertexOnJetAxis* firstVertex,
703  const VxVertexOnJetAxis* secondVertex,
704  const VxJetCandidate* myJetCandidate) const {
705 
706  //first get a copy of all vertex positions (this you can't avoid I fear...)
707  RecVertexPositions copyOfRecVertexPositions(myJetCandidate->getRecVertexPositions());
708  const FitQuality & copyOfRecVertexQuality=copyOfRecVertexPositions.fitQuality();
709 
710  double oldchi2=copyOfRecVertexQuality.chiSquared();
711  double oldndf=copyOfRecVertexQuality.numberDoF();
712 
713  //#ifdef JetFitterRoutines_DEBUG2
714  const Amg::VectorX & positionVector=copyOfRecVertexPositions.position();
715  const Amg::MatrixX & positionCov=copyOfRecVertexPositions.covariancePosition();
716  double phiold=positionVector(Trk::jet_phi);
717  double thetaold=positionVector(Trk::jet_theta);
718  double phierr=std::sqrt(positionCov(Trk::jet_phi,Trk::jet_phi));
719  double thetaerr=std::sqrt(positionCov(Trk::jet_theta,Trk::jet_theta));
720 
721 
722  //now do the merging of the second cluster to the primary vertex...
723  m_helper->performKalmanConstraintToMergeVertices(copyOfRecVertexPositions,
724  *firstVertex,
725  *secondVertex);
726 
727  double phinew=positionVector[Trk::jet_phi];
728  double thetanew=positionVector[Trk::jet_theta];
729  bool isshifted=std::pow((phinew-phiold)/phierr,2)+std::pow((thetanew-thetaold)/thetaerr,2)>m_maxDRshift* m_maxDRshift;
730 
731  double chi2 = firstVertex->fitQuality().chiSquared() +
732  secondVertex->fitQuality().chiSquared() +
733  copyOfRecVertexPositions.fitQuality().chiSquared()-oldchi2;
734 
735 
736  double ndf = copyOfRecVertexPositions.fitQuality().numberDoF()-oldndf +
737  firstVertex->fitQuality().numberDoF()+
738  secondVertex->fitQuality().numberDoF();
739 
740  if (chi2<0||ndf<0) {
741  ATH_MSG_WARNING ("In the compatibility estimation chi2: " << chi2 << " ndf " << ndf << " giving back 0 prob");
742  return std::pair<double,bool>(0,isshifted);
743  }
744 
745  return std::pair<double,bool>(TMath::Prob(chi2,(int)std::floor(ndf+0.5)),isshifted);
746 
747 
748  }
749 
750  double JetFitterRoutines::fullProbabilityOfMerging(const VxVertexOnJetAxis* firstVertex,
751  const VxVertexOnJetAxis* secondVertex,
752  const VxJetCandidate* myJetCandidate,
753  int num_maxiterations,
754  bool treat_sign_flip,
755  int num_signflip_maxiterations,
756  double deltachi2_convergence) const {
757 
758 
759  if (firstVertex==nullptr||secondVertex==nullptr||myJetCandidate==nullptr) {
760  ATH_MSG_WARNING ("zero pointer given to the full probability estimation. No estimation performed, zero prob returned ");
761  return 0;
762  }
763 
764  //copy the VxJetCandidate into a new object (this is expensive...)
765  VxJetCandidate newJetCandidate(*myJetCandidate);
766 
767 
768  //first find correspondence between old first and secondVertex and new ones (horrible,but...)
769 
770  std::map<const VxVertexOnJetAxis*,VxVertexOnJetAxis*> oldToNewVtxPointers;
771 
772 
773  const std::vector<VxVertexOnJetAxis*> vectorOfOldJetCand=myJetCandidate->getVerticesOnJetAxis();
774  const std::vector<VxVertexOnJetAxis*> vectorOfNewJetCand=newJetCandidate.getVerticesOnJetAxis();
775 
776  const VxVertexOnJetAxis* primaryOfFirst=myJetCandidate->getPrimaryVertex();
777  VxVertexOnJetAxis* primaryOfSecond=newJetCandidate.getPrimaryVertex();
778 
779  if (primaryOfFirst==nullptr||primaryOfSecond==nullptr) {
780  ATH_MSG_WARNING ("Empty primary vertex found when estimating fullProbOfMerging. 0 prob returned.");
781  return 0;
782  }
783 
784  oldToNewVtxPointers[primaryOfFirst]=primaryOfSecond;
785 
786  unsigned int sizeOfVertices=vectorOfOldJetCand.size();
787  if (vectorOfNewJetCand.size()!=sizeOfVertices) {
788  ATH_MSG_WARNING ("Old and new track of vertices do not match during fullProbOfMerging. 0 prob returned.");
789  return 0;
790  }
791 
792  for (unsigned int s=0;s<sizeOfVertices;s++) {
793  const VxVertexOnJetAxis* pointer1=vectorOfOldJetCand[s];
794  VxVertexOnJetAxis* pointer2=vectorOfNewJetCand[s];
795  if (pointer1==nullptr||pointer2==nullptr) {
796  ATH_MSG_WARNING ("One of the pointers of the original or copied vector of vertices is empty during fullProbOfMerging. Skipping it...");
797  } else {
798  oldToNewVtxPointers[pointer1]=pointer2;
799  }
800  }
801 
802  //now merge firstVertex and secondVertex
803  VxVertexOnJetAxis* newFirstVertex=oldToNewVtxPointers[firstVertex];
804  VxVertexOnJetAxis* newSecondVertex=oldToNewVtxPointers[secondVertex];
805 
806  if (newFirstVertex==nullptr||newSecondVertex==nullptr) {
807  ATH_MSG_WARNING ("No correspondence to the given firstVertex or secondVertex in fullProbOfMerging. Returning 0 prob.");
808  return 0.;
809  }
810 
811  VxVertexOnJetAxis & commonVertex=m_helper->mergeVerticesInJetCandidate(*newFirstVertex,
812  *newSecondVertex,
813  newJetCandidate);
814 
815  //now you need to update the numbering scheme
816  m_initializationHelper->updateTrackNumbering(&newJetCandidate);
817 
818  performTheFit(&newJetCandidate,num_maxiterations,
819  treat_sign_flip,
820  num_signflip_maxiterations,
821  deltachi2_convergence);
822 
823  // const FitQuality & qualityOfMergedVertex=commonVertex.fitQuality();
824 
825 #ifdef JetFitterRoutines_DEBUG
826  ATH_MSG_DEBUG( " End estimating full prob of merging... chi2 " << commonVertex.fitQuality().chiSquared() << " ndf " << commonVertex.fitQuality().numberDoF() );
827 #endif
828 
829  return (double)TMath::Prob(commonVertex.fitQuality().chiSquared(),
830  (int)std::floor(commonVertex.fitQuality().numberDoF()+0.5));
831 
832  }
833 
834  void JetFitterRoutines::fillTableWithFullProbOfMerging(VxJetCandidate* myJetCandidate,
835  int num_maxiterations,
836  bool treat_sign_flip,
837  int num_signflip_maxiterations,
838  double deltachi2_convergence,
839  double threshold_probability) const {
840  fillTableWithProbOfMerging(myJetCandidate,
841  true,
842  num_maxiterations,
843  treat_sign_flip,
844  num_signflip_maxiterations,
845  deltachi2_convergence,
846  threshold_probability);
847  }
848 
849  void JetFitterRoutines::fillTableWithFastProbOfMerging(VxJetCandidate* myJetCandidate) const {
850  fillTableWithProbOfMerging(myJetCandidate,false);
851  }
852 
853 
854  void JetFitterRoutines::fillTableWithProbOfMerging(VxJetCandidate* myJetCandidate,
855  bool fullcomputation,
856  int num_maxiterations,
857  bool treat_sign_flip,
858  int num_signflip_maxiterations,
859  double deltachi2_convergence,
860  double threshold_probability) const {
861 
862  if (myJetCandidate==nullptr) {
863  ATH_MSG_WARNING( "VxJetCandidate provided is a zero pointer. No compatibility table calculated." );
864  return;
865  }
866 
867  //first create the compatibility table object...
868  Trk::VxClusteringTable* & clusteringTablePtr(myJetCandidate->getClusteringTable());
869  if (clusteringTablePtr!=nullptr) {
870  delete clusteringTablePtr;
871  }
872  clusteringTablePtr=new Trk::VxClusteringTable();
873 
874  double highestprobability(0.);
875 
876  VxVertexOnJetAxis* primaryVertex=myJetCandidate->getPrimaryVertex();
877 
878  if (primaryVertex==nullptr) {
879  ATH_MSG_WARNING( "VxJetCandidate provided has no primary vertex. No compatibility table calculated." );
880  return;
881  }
882 
883  primaryVertex->setCompatibilityToPrimaryVtx(1);//stupid but assign prob 1 to primary vtx for consistency
884 
885  const std::vector<VxVertexOnJetAxis*> & tracksOnJetAxis=myJetCandidate->getVerticesOnJetAxis();
886 
887  //now evaluate probability of cluster forming with primary vertex
888  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxBegin=tracksOnJetAxis.begin();
889  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxEnd=tracksOnJetAxis.end();
890 
891  for (std::vector<VxVertexOnJetAxis*>::const_iterator VtxIter=VtxBegin;
892  VtxIter!=VtxEnd;++VtxIter) {
893 
894  std::pair<double,bool> fastProbabilityAndNonLinearity=fastProbabilityOfMerging(primaryVertex,
895  *VtxIter,
896  myJetCandidate);
897 
898 
899  ATH_MSG_VERBOSE ("Fast probability of merging between primary and " <<
900  (*VtxIter)->getNumVertex() << " is " << fastProbabilityAndNonLinearity.first);
901 
902 #ifdef JetFitterRoutines_DEBUG2
903  ATH_MSG_DEBUG( "Fast probability of merging between primary and " <<
904  (*VtxIter)->getNumVertex() << " is " << fastProbabilityAndNonLinearity.first << " and is max DR " <<
905  fastProbabilityAndNonLinearity.second );
906 #endif
907 
908  if (fullcomputation) {
909  if (fastProbabilityAndNonLinearity.first>threshold_probability) {
910  if (fastProbabilityAndNonLinearity.first>highestprobability/100.&&fastProbabilityAndNonLinearity.second) {
911 
912  double fullProbability=fullProbabilityOfMerging(primaryVertex,*VtxIter,
913  myJetCandidate,num_maxiterations,
914  treat_sign_flip,
915  num_signflip_maxiterations,
916  deltachi2_convergence);
917 
918 
919 #ifdef JetFitterRoutines_DEBUG2
920  ATH_MSG_DEBUG( "Full probability of merging with primary is " << fullProbability );
921 #endif
922 
923  ATH_MSG_DEBUG ("Full probability of merging with primary is " << fullProbability);
924 
925  //store this fullProbability into the VxClusteringTable of myJetCandidate
926 
927  clusteringTablePtr->setCompatibilityOfTo(PairOfVxVertexOnJetAxis(primaryVertex,*VtxIter),fullProbability);
928  if (fullProbability>highestprobability) {
929  highestprobability=fullProbability;
930  }
931  (*VtxIter)->setCompatibilityToPrimaryVtx(fullProbability);
932  } else {
933  clusteringTablePtr->setCompatibilityOfTo(PairOfVxVertexOnJetAxis(primaryVertex,*VtxIter),fastProbabilityAndNonLinearity.first);
934  (*VtxIter)->setCompatibilityToPrimaryVtx(fastProbabilityAndNonLinearity.first);
935  }
936  } else {
937  (*VtxIter)->setCompatibilityToPrimaryVtx(fastProbabilityAndNonLinearity.first);
938  }
939  } else {
940  if (fastProbabilityAndNonLinearity.first>threshold_probability) {
941  clusteringTablePtr->setCompatibilityOfTo(PairOfVxVertexOnJetAxis(primaryVertex,*VtxIter),fastProbabilityAndNonLinearity.first);
942  if (fastProbabilityAndNonLinearity.first>highestprobability) {
943  highestprobability=fastProbabilityAndNonLinearity.first;
944  }
945  }
946  (*VtxIter)->setCompatibilityToPrimaryVtx(fastProbabilityAndNonLinearity.first);
947  }
948  }//end first iteration over tracks for compatibility with primary vertex
949 
950  //now all the remaining combination have to be tested...
951 
952  for (std::vector<Trk::VxVertexOnJetAxis*>::const_iterator VtxIter2=VtxBegin;
953  VtxIter2!=VtxEnd;++VtxIter2) {
954  for (std::vector<Trk::VxVertexOnJetAxis*>::const_iterator VtxIter1=VtxBegin;
955  VtxIter1!=VtxIter2;++VtxIter1) {
956 
957  std::pair<double,bool> fastProbabilityAndNonLinearity=fastProbabilityOfMerging(*VtxIter1,
958  *VtxIter2,
959  myJetCandidate);
960 
961 
962 
963  ATH_MSG_VERBOSE ("Fast probability of merging between vtx n " << (*VtxIter1)->getNumVertex()<<
964  " and " << (*VtxIter2)->getNumVertex() << " is " <<
965  fastProbabilityAndNonLinearity.first << " and is max DR " <<
966  fastProbabilityAndNonLinearity.second);
967 
968 #ifdef JetFitterRoutines_DEBUG2
969  ATH_MSG_DEBUG( "Fast probability of merging between vtx n " << (*VtxIter1)->getNumVertex() << " and " <<
970  (*VtxIter2)->getNumVertex() << " is " << fastProbabilityAndNonLinearity.first << " and is max DR " <<
971  fastProbabilityAndNonLinearity.second );
972 #endif
973  if (fullcomputation) {
974  if (fastProbabilityAndNonLinearity.first>threshold_probability) {
975  if (fastProbabilityAndNonLinearity.first>highestprobability/100.&&fastProbabilityAndNonLinearity.second) {
976  double fullProbability=fullProbabilityOfMerging(*VtxIter1,*VtxIter2,
977  myJetCandidate,num_maxiterations,
978  treat_sign_flip,
979  num_signflip_maxiterations,
980  deltachi2_convergence);
981 
982 
983 #ifdef JetFitterRoutines_DEBUG2
984  ATH_MSG_DEBUG( "Full probability of merging is " << fullProbability );
985 #endif
986 
987  ATH_MSG_VERBOSE ("Full probability of merging is " << fullProbability);
988 
989  //store this fullProbability into the VxClusteringTable of myJetCandidate
990 
991  clusteringTablePtr->setCompatibilityOfTo(PairOfVxVertexOnJetAxis(*VtxIter1,*VtxIter2),fullProbability);
992  if (fullProbability>highestprobability) {
993  highestprobability=fullProbability;
994  }
995  } else {
996  clusteringTablePtr->setCompatibilityOfTo(PairOfVxVertexOnJetAxis(*VtxIter1,*VtxIter2),fastProbabilityAndNonLinearity.first);
997  }
998  }
999  } else {
1000  if (fastProbabilityAndNonLinearity.first>threshold_probability) {
1001  clusteringTablePtr->setCompatibilityOfTo(PairOfVxVertexOnJetAxis(*VtxIter1,*VtxIter2),fastProbabilityAndNonLinearity.first);
1002  if (fastProbabilityAndNonLinearity.first>highestprobability) {
1003  highestprobability=fastProbabilityAndNonLinearity.first;
1004  }
1005  }
1006  }
1007  }//end first vtx loop
1008  }//end second vtx loop
1009 
1010  }//end fillTableWithFullProbOfMerging() method
1011 
1012  void JetFitterRoutines::deleteVertexFromJetCandidate(VxVertexOnJetAxis* vertexToDelete,
1013  VxJetCandidate* myJetCandidate) const {
1014  if (vertexToDelete==myJetCandidate->getPrimaryVertex()) {
1015  ATH_MSG_WARNING ("YOU ARE deleting the primary vertex. This is not possible... ");
1016  return;
1017  }
1018 
1019  // now you need to *delete* the second vertex in copyOfRecVertexPositions and
1020  // in copyOfLinearizationPositions
1021  const Amg::VectorX & recPosition=myJetCandidate->getRecVertexPositions().position();
1022  const Amg::VectorX & linPosition=myJetCandidate->getLinearizationVertexPositions().position();
1023  const Amg::VectorX & constraintPosition=myJetCandidate->getConstraintVertexPositions().position();
1024  const Amg::MatrixX & covPosition=myJetCandidate->getRecVertexPositions().covariancePosition();
1025  const Amg::MatrixX & covConstraintPosition=myJetCandidate->getConstraintVertexPositions().covariancePosition();
1026 
1027  int numbVertex=numRow(vertexToDelete->getNumVertex());
1028  //call the function which deletes a single row (they're in the anonymous namespace, in Utilities.h)
1029  Amg::VectorX reducedRecPositions=deleteRowFromVector(recPosition,numbVertex);
1030  Amg::VectorX reducedLinPositions=deleteRowFromVector(linPosition,numbVertex);
1031  Amg::VectorX reducedConstraintPositions=deleteRowFromVector(constraintPosition,numbVertex);
1032  Amg::MatrixX reducedCovPositions=deleteRowFromSymMatrix(covPosition,numbVertex);
1033  Amg::MatrixX reducedConstraintCovPositions=deleteRowFromSymMatrix(covConstraintPosition,numbVertex);
1034 
1035  myJetCandidate->setRecVertexPositions(RecVertexPositions(reducedRecPositions,
1036  reducedCovPositions,
1037  0.,0.));
1038  myJetCandidate->setConstraintVertexPositions(RecVertexPositions(reducedConstraintPositions,
1039  reducedConstraintCovPositions,
1040  0.,0.));
1041  myJetCandidate->setLinearizationVertexPositions(VertexPositions(reducedLinPositions));
1042 
1043 
1044 
1045  //make a copy of the VxTrackAtVertex vector in myJetCandidate to delete some of the elements
1046 
1047  //PAY ATTENTION: here you are modifying the vector of tracks DIRECTLY...
1048  std::vector<VxTrackAtVertex*>* tracksAtJetCandidate(myJetCandidate->vxTrackAtVertex());
1049  // RS 19.04.2011 attempt to fix coverity defect 22750
1050  // const std::vector<VxTrackAtVertex*>::iterator TracksBegin=tracksAtJetCandidate->begin();
1051  // std::vector<VxTrackAtVertex*>::iterator TracksEnd=tracksAtJetCandidate->end();
1052 
1053 
1054  const std::vector<VxTrackAtVertex*> & tracksAtVtx(vertexToDelete->getTracksAtVertex());
1055 
1056  const std::vector<VxTrackAtVertex*>::const_iterator TracksAtVtxBegin=tracksAtVtx.begin();
1057  const std::vector<VxTrackAtVertex*>::const_iterator TracksAtVtxEnd=tracksAtVtx.end();
1058 
1059  int numberOfTracksBefore=tracksAtJetCandidate->size();
1060  int numberOfTracksToDelete=tracksAtVtx.size();
1061 
1062  for (std::vector<VxTrackAtVertex*>::const_iterator TracksAtVtxIter=TracksAtVtxBegin;
1063  TracksAtVtxIter!=TracksAtVtxEnd;
1064  ++TracksAtVtxIter) {
1065 
1066  // RS 19.04.2011 attempt to fix coverity defect 22750
1067  std::vector<VxTrackAtVertex*>::iterator TracksBegin=tracksAtJetCandidate->begin();
1068  std::vector<VxTrackAtVertex*>::iterator TracksEnd=tracksAtJetCandidate->end();
1069  for (std::vector<VxTrackAtVertex*>::iterator TracksIter=TracksBegin;TracksIter!=TracksEnd;) {
1070 
1071  if (*TracksIter==*TracksAtVtxIter) {
1072  Trk::VxTrackAtVertex* oldPointer=*TracksIter;
1073  TracksIter=tracksAtJetCandidate->erase(TracksIter);
1074  delete oldPointer;
1075  oldPointer=nullptr;
1076  TracksEnd=tracksAtJetCandidate->end();
1077  break;
1078  }
1079  ++TracksIter;
1080 
1081  }
1082  }
1083 
1084  if (numberOfTracksBefore-numberOfTracksToDelete!=(int)myJetCandidate->vxTrackAtVertex()->size()) {
1085  ATH_MSG_DEBUG( " MISMATCH in JetFitterRoutines: the jetcandidate had: " << numberOfTracksBefore << " tracks " <<
1086  " and " << numberOfTracksToDelete << " to delete = " << myJetCandidate->vxTrackAtVertex()->size() << " tracks left! " );
1087  }
1088 
1089 
1090  //now as a last step you need to delete the vertex you're not using anymore...
1091 
1092  std::vector<VxVertexOnJetAxis*> copyOfVerticesAtJetCandidate=myJetCandidate->getVerticesOnJetAxis();
1093 
1094  const std::vector<VxVertexOnJetAxis*>::iterator VerticesBegin=copyOfVerticesAtJetCandidate.begin();
1095  std::vector<VxVertexOnJetAxis*>::iterator VerticesEnd=copyOfVerticesAtJetCandidate.end();
1096 
1097  bool found=false;
1098 
1099  //have a variable to store the vector of tracks which then need to be removed (copy)...
1100  // std::vector<VxTrackAtVertex*> tracksToRemove=vertexToDelete->getTracksAtVertex();
1101 
1102  for (std::vector<VxVertexOnJetAxis*>::iterator VerticesIter=VerticesBegin;VerticesIter!=VerticesEnd;) {
1103  if ((*VerticesIter)==vertexToDelete) {
1104  delete *VerticesIter;
1105  VerticesIter=copyOfVerticesAtJetCandidate.erase(VerticesIter);
1106  VerticesEnd=copyOfVerticesAtJetCandidate.end();
1107  found=true;
1108  break;
1109  }
1110  ++VerticesIter;
1111 
1112  }
1113 
1114  if (!found) {
1115  ATH_MSG_WARNING( "Could not find the vertex to delete... Very strange... Check!!! " );
1116  }
1117 
1118  //update myJetCandidate with the new vector of tracks
1119  myJetCandidate->setVerticesOnJetAxis(copyOfVerticesAtJetCandidate);
1120 
1121  //update the numbering scheme
1122  m_initializationHelper->updateTrackNumbering(myJetCandidate);
1123 
1124  }
1125 
1126  void JetFitterRoutines::copyRecoPositionsToLinearizationPositions(VxJetCandidate & myJetCandidate) const
1127  {
1128 
1129  const VertexPositions & newLinVertexPositions=myJetCandidate.getRecVertexPositions();
1130  Amg::VectorX linPositions=newLinVertexPositions.position();
1131 
1132  const std::vector<VxVertexOnJetAxis*> & associatedVertices=myJetCandidate.getVerticesOnJetAxis();
1133  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxBegin=associatedVertices.begin();
1134  const std::vector<VxVertexOnJetAxis*>::const_iterator VtxEnd=associatedVertices.end();
1135 
1136  for (std::vector<VxVertexOnJetAxis*>::const_iterator VtxIter=VtxBegin;VtxIter!=VtxEnd;++VtxIter) {
1137  VxVertexOnJetAxis* myVertex=(*VtxIter);
1138  if (myVertex!=nullptr) {
1139 
1140  double distOnAxis=linPositions[numRow(myVertex->getNumVertex())];
1141 
1142  auto const sinLinJetTheta = std::sin(linPositions[Trk::jet_theta]);
1143  auto const cosLinJetTheta = std::cos(linPositions[Trk::jet_theta]);
1144 
1145  auto const absLinJetTheta = std::abs(linPositions[Trk::jet_theta]);
1146  auto const abssinLinJetTheta = std::abs(sinLinJetTheta);
1147  auto const abscosLinJetTheta = std::abs(cosLinJetTheta);
1148 
1149  double R=distOnAxis*sinLinJetTheta;
1150  double Z=distOnAxis*cosLinJetTheta;
1151  if (std::abs(R)>m_maxR)
1152  {
1153  if (absLinJetTheta>1e-8)
1154  {
1155  ATH_MSG_DEBUG (" Closest distance of track to jet axis is outside ID envelope, R=" << R << ", setting to R= " << m_maxR);
1156  distOnAxis=m_maxR / abssinLinJetTheta;
1157  }
1158  }
1159 
1160  Z=distOnAxis*cosLinJetTheta;
1161  if (std::abs(Z)>m_maxZ)
1162  {
1163  if (abscosLinJetTheta>1e-8)
1164  {
1165  ATH_MSG_DEBUG( " Closest distance of track to jet axis is outside ID envelope, Z=" << Z << ", setting to Z= " << m_maxZ );
1166  distOnAxis=m_maxZ / cosLinJetTheta;
1167  }
1168  }
1169 
1170  linPositions[numRow(myVertex->getNumVertex())]=distOnAxis;
1171  }
1172  }
1173 
1174  VertexPositions & linVertexPositions=myJetCandidate.getLinearizationVertexPositions();
1175  linVertexPositions.setPosition(linPositions);
1176  //now set the linearization position for the next step to the actual fitted vertex
1177  //OLD CODE BEFORE CHECKS
1178  //myJetCandidate->setLinearizationVertexPositions(myJetCandidate->getRecVertexPositions());
1179  }
1180 
1181 }//end namespace
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
Trk::VertexPositions::position
const Amg::VectorX & position() const
return position of vertex
Definition: VertexPositions.cxx:95
covarianceTool.ndf
ndf
Definition: covarianceTool.py:678
RecVertex.h
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
JetFitterRoutines.h
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
Trk::VxJetCandidate::getVerticesOnJetAxis
const std::vector< VxVertexOnJetAxis * > & getVerticesOnJetAxis(void) const
Definition: VxJetCandidate.cxx:543
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition: EventPrimitives.h:30
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:27
get_generator_info.result
result
Definition: get_generator_info.py:21
Trk::VxJetFitterDebugInfo::setSignFlipNumFitIterations
void setSignFlipNumFitIterations(int num)
Definition: VxJetFitterDebugInfo.h:88
Trk::VxTrackAtVertex
The VxTrackAtVertex is a common class for all present TrkVertexFitters The VxTrackAtVertex is designe...
Definition: VxTrackAtVertex.h:77
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
AtlasFieldCacheCondObj
Definition: AtlasFieldCacheCondObj.h:19
Trk::jet_theta
@ jet_theta
Definition: JetVtxParamDefs.h:28
Monitored::Z
@ Z
Definition: HistogramFillerUtils.h:24
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
Trk::VxJetCandidate::getLinearizationVertexPositions
const Trk::VertexPositions & getLinearizationVertexPositions() const
Definition: VxJetCandidate.cxx:531
initialize
void initialize()
Definition: run_EoverP.cxx:894
JetVtxParamDefs.h
VxVertexOnJetAxis.h
VxJetCandidate.h
Trk::VxVertexOnJetAxis
VxVertexOnJetAxis inherits from Vertex.
Definition: VxVertexOnJetAxis.h:79
Trk::jet_xv
@ jet_xv
Definition: JetVtxParamDefs.h:27
NeutralParameters.h
VxJetFitterDebugInfo.h
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
TrkDistanceFinderNeutralCharged.h
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
dbg::ptr
void * ptr(T *p)
Definition: SGImplSvc.cxx:74
JetFitterHelper.h
Trk::VxJetCandidate::getRecVertexPositions
const Trk::RecVertexPositions & getRecVertexPositions() const
Definition: VxJetCandidate.cxx:519
VxClusteringTable.h
JetFitterInitializationHelper.h
Trk::VxJetFitterDebugInfo::setNumFitIterations
void setNumFitIterations(int num)
Definition: VxJetFitterDebugInfo.h:80
Trk::VxJetCandidate::getPrimaryVertex
const VxVertexOnJetAxis * getPrimaryVertex(void) const
Definition: VxJetCandidate.cxx:551
GeoPrimitives.h
Trk::VertexPositions::setPosition
void setPosition(const Amg::VectorX &)
Definition: VertexPositions.cxx:109
PairOfVxVertexOnJetAxis.h
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Trk::jet_yv
@ jet_yv
Definition: JetVtxParamDefs.h:27
Trk::VxClusteringTable
Definition: VxClusteringTable.h:71
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:523
VxTrackAtVertex.h
Utilities.h
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::VxVertexOnJetAxis::fitQuality
const Trk::FitQuality & fitQuality() const
Fit quality access method.
Definition: VxVertexOnJetAxis.cxx:86
Trk::FitQuality
Class to represent and store fit qualities from track reconstruction in terms of and number of degre...
Definition: FitQuality.h:97
AnalysisUtils::Delta::R
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
Definition: AnalysisMisc.h:49
Trk::RecVertexPositions::fitQuality
const Trk::FitQuality & fitQuality() const
Fit quality access method.
Definition: RecVertexPositions.cxx:167
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
Trk::VxVertexOnJetAxis::setCompatibilityToPrimaryVtx
void setCompatibilityToPrimaryVtx(float)
set compatibility to the primary vertex
Definition: VxVertexOnJetAxis.cxx:121
Trk::JetFitterRoutines::JetFitterRoutines
JetFitterRoutines(const std::string &t, const std::string &n, const IInterface *p)
Constructor.
Definition: JetFitterRoutines.cxx:74
Trk::jet_zv
@ jet_zv
position x,y,z of primary vertex
Definition: JetVtxParamDefs.h:27
KalmanVertexOnJetAxisSmoother.h
Trk::VxJetCandidate
Definition: VxJetCandidate.h:72
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::RecVertexPositions
Definition: RecVertexPositions.h:34
Trk::NeutralTrack
Definition: NeutralTrack.h:10
Trk::RecVertexPositions::covariancePosition
Amg::MatrixX const & covariancePosition() const
return the covDeltaV matrix of the vertex fit
Definition: RecVertexPositions.cxx:171
a
TList * a
Definition: liststreamerinfos.cxx:10
Trk::VxJetCandidate::getClusteringTable
Trk::VxClusteringTable *& getClusteringTable(void)
Definition: VxJetCandidate.cxx:569
CondAlgsOpts.found
int found
Definition: CondAlgsOpts.py:101
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
TrkDistanceFinderNeutralNeutral.h
MagField::AtlasFieldCache
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h)
Definition: AtlasFieldCache.h:43
Trk::FitQuality::chiSquared
double chiSquared() const
returns the of the overall track fit
Definition: FitQuality.h:56
Trk::jet_phi
@ jet_phi
Definition: JetVtxParamDefs.h:28
Trk::FitQuality::numberDoF
int numberDoF() const
returns the number of degrees of freedom of the overall track or vertex fit as integer
Definition: FitQuality.h:60
Trk::VertexPositions
VertexPositions class to represent and store a vertex.
Definition: VertexPositions.h:36
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
Trk::VxJetFitterDebugInfo
Definition: VxJetFitterDebugInfo.h:44
python.AutoConfigFlags.msg
msg
Definition: AutoConfigFlags.py:7
RecVertexPositions.h
python.Dumpers.FitQuality
FitQuality
Definition: Dumpers.py:63
KalmanVertexOnJetAxisUpdator.h
Generated on Sun Dec 22 2024 21:12:30 for ATLAS Offline Software by   doxygen 1.8.18