diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/CMakeLists.txt b/Tracking/TrkVertexFitter/TrkVertexFitters/CMakeLists.txt index 9cbb8c70af470c13cfed07c6cf388659f02e64b4..cf3e86f316b625932d3528d11e6f188d4dc16d43 100644 --- a/Tracking/TrkVertexFitter/TrkVertexFitters/CMakeLists.txt +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/CMakeLists.txt @@ -33,3 +33,18 @@ atlas_add_component( TrkVertexFitters src/components/*.cxx LINK_LIBRARIES AthenaBaseComps xAODTracking GaudiKernel TrkParameters TrkParametersBase TrkParticleBase VxVertex TrkVertexFitterInterfaces TrkSurfaces TrkLinks TrkTrack VxMultiVertex TrkExInterfaces TrkVertexFittersLib ) + +# Install files from the package: +atlas_install_joboptions( share/*.py ) + + +atlas_add_test( AdaptiveVertexFitter_test + SCRIPT athena.py TrkVertexFitters/AdaptiveVertexFitter_test.py + PROPERTIES TIMEOUT 300 + EXTRA_PATTERNS " INFO |WARNING |found service|Adding private|^ +[+]|HepPDT Version|^indet|^pixel|^phi|^eta|^GEOPIXEL|in material|no dictionary|^subdet|^part|^barrel_endcap|^layer|^sct|^bec|^side0|^strip0|^lay_disk|TrackingGeometrySvc|^PixelID" ) + + +atlas_add_test( AdaptiveMultiVertexFitter_test + SCRIPT athena.py TrkVertexFitters/AdaptiveMultiVertexFitter_test.py + PROPERTIES TIMEOUT 300 + EXTRA_PATTERNS " INFO |WARNING |found service|Adding private|^ +[+]|HepPDT Version|^indet|^pixel|^phi|^eta|^GEOPIXEL|in material|no dictionary|^subdet|^part|^barrel_endcap|^layer|^sct|^bec|^side0|^strip0|^lay_disk|TrackingGeometrySvc|^PixelID" ) diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveMultiVertexFitter_test.py b/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveMultiVertexFitter_test.py new file mode 100644 index 0000000000000000000000000000000000000000..b5527ad49615027245dfc3476a08fce6049ba7dd --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveMultiVertexFitter_test.py @@ -0,0 +1,64 @@ +# +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration. +# +# File: TrkVertexFitters/share/AdaptiveMultiVertexFitter_test.py +# Author: scott snyder <snyder@bnl.gov> +# Data: Jul, 2019 +# Brief: Unit test for AdaptiveMultiVertexFitter. Incomplete! +# + + +from __future__ import print_function + + +from AthenaCommon.DetFlags import DetFlags +DetFlags.detdescr.ID_setOn() + +RunNumber = 284500 + +import sys +import string +import ROOT +import math +from AtlasGeoModel import SetGeometryVersion +from AtlasGeoModel import GeoModelInit +from AtlasGeoModel import SetupRecoGeometry + +from GeoModelSvc.GeoModelSvcConf import GeoModelSvc +ServiceMgr += GeoModelSvc() +theApp.CreateSvc += [ "GeoModelSvc"] +from AtlasGeoModel import InDetGM + +from IOVDbSvc.IOVDbSvcConf import IOVDbSvc +IOVDbSvc().GlobalTag = 'OFLCOND-RUN12-SDR-35' + +import MagFieldServices.SetupField + +from TrkDetDescrSvc.AtlasTrackingGeometrySvc import AtlasTrackingGeometrySvc + +from AthenaCommon.AlgSequence import AlgSequence +topSequence = AlgSequence() + +theApp.EvtMax=2 + +from xAODEventInfoCnv.xAODEventInfoCnvConf import xAODMaker__EventInfoCnvAlg +eialg = xAODMaker__EventInfoCnvAlg (DoBeginRun = False) +topSequence += eialg + + + +# Suppress useless GeoModelSvc messages. +from AthenaCommon import Constants +GeoModelSvc().OutputLevel=Constants.WARNING + +from TrkExTools.AtlasExtrapolator import AtlasExtrapolator + + +from TrkVertexFitters.TrkVertexFittersConf import \ + Trk__AdaptiveMultiVertexFitterTestAlg, Trk__AdaptiveMultiVertexFitter +fitter = Trk__AdaptiveMultiVertexFitter ('AdaptiveMultiVertexFitter', + OutputLevel = INFO) +testalg1 = Trk__AdaptiveMultiVertexFitterTestAlg ('testalg1', + OutputLevel = VERBOSE, + Tool = fitter) +topSequence += testalg1 diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveMultiVertexFitter_test.ref b/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveMultiVertexFitter_test.ref new file mode 100644 index 0000000000000000000000000000000000000000..973203c71168ca6a9b3935010d4d19c29d331d59 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveMultiVertexFitter_test.ref @@ -0,0 +1,2 @@ +testalg1 VERBOSE execute +testalg1 VERBOSE execute diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveVertexFitter_test.py b/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveVertexFitter_test.py new file mode 100644 index 0000000000000000000000000000000000000000..083ac972875aa9c12f01f6e6c0af0255b5f1a0d6 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveVertexFitter_test.py @@ -0,0 +1,67 @@ +# +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration. +# +# File: TrkVertexFitters/share/AdaptiveVertexFitter_test.py +# Author: scott snyder <snyder@bnl.gov> +# Data: Jul, 2019 +# Brief: Unit test for AdaptiveVertexFitter. Incomplete! +# + + +from __future__ import print_function + + +from AthenaCommon.DetFlags import DetFlags +DetFlags.detdescr.ID_setOn() + +RunNumber = 284500 + +import sys +import string +import ROOT +import math +from AtlasGeoModel import SetGeometryVersion +from AtlasGeoModel import GeoModelInit +from AtlasGeoModel import SetupRecoGeometry + +from GeoModelSvc.GeoModelSvcConf import GeoModelSvc +ServiceMgr += GeoModelSvc() +theApp.CreateSvc += [ "GeoModelSvc"] +from AtlasGeoModel import InDetGM + +from IOVDbSvc.IOVDbSvcConf import IOVDbSvc +IOVDbSvc().GlobalTag = 'OFLCOND-RUN12-SDR-35' + +import MagFieldServices.SetupField + +from TrkDetDescrSvc.AtlasTrackingGeometrySvc import AtlasTrackingGeometrySvc + +from AthenaCommon.AlgSequence import AlgSequence +topSequence = AlgSequence() + +theApp.EvtMax=2 + +from xAODEventInfoCnv.xAODEventInfoCnvConf import xAODMaker__EventInfoCnvAlg +eialg = xAODMaker__EventInfoCnvAlg (DoBeginRun = False) +topSequence += eialg + + + +# Suppress useless GeoModelSvc messages. +from AthenaCommon import Constants +GeoModelSvc().OutputLevel=Constants.WARNING + +#from AthenaCommon.GlobalFlags import globalflags +#from InDetRecExample.InDetJobProperties import InDetFlags +#include ('InDetRecExample/InDetRecLoadTools.py') +from TrkExTools.AtlasExtrapolator import AtlasExtrapolator + + +from TrkVertexFitters.TrkVertexFittersConf import \ + Trk__AdaptiveVertexFitterTestAlg, Trk__AdaptiveVertexFitter +fitter = Trk__AdaptiveVertexFitter ('AdaptiveVertexFitter', + OutputLevel = INFO) +testalg1 = Trk__AdaptiveVertexFitterTestAlg ('testalg1', + OutputLevel = VERBOSE, + Tool = fitter) +topSequence += testalg1 diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveVertexFitter_test.ref b/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveVertexFitter_test.ref new file mode 100644 index 0000000000000000000000000000000000000000..973203c71168ca6a9b3935010d4d19c29d331d59 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/share/AdaptiveVertexFitter_test.ref @@ -0,0 +1,2 @@ +testalg1 VERBOSE execute +testalg1 VERBOSE execute diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitter.cxx b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitter.cxx index ec90d1b95be7e3f8deb4108c858c3217754a23fa..fb3eec7bed9c45a1f6c2c25e6eea3ca8bd67c8d6 100755 --- a/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitter.cxx +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitter.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration */ @@ -371,7 +371,7 @@ namespace Trk // std::cout << " step: " << num_steps << " considering shift of: " << *iter << " vtx z: " << (*iter)->recVertex().position().z() <<std::endl; - Amg::Vector3D vrtpos(3,0); + Amg::Vector3D vrtpos; vrtpos[0]=oldpositions[*iter][0]-(*iter)->position()[0]; vrtpos[1]=oldpositions[*iter][1]-(*iter)->position()[1]; diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitterTestAlg.cxx b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitterTestAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..223ca0c249fac60383b4c3f686ad6d6685b473d4 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitterTestAlg.cxx @@ -0,0 +1,364 @@ +/* + * Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration. + */ +/** + * @file TrkVertexFitters/src/AdaptiveMultiVertexFitterTestAlg.cxx + * @author scott snyder <snyder@bnl.gov> + * @date Jul, 2019 + * @brief Algorithm for testing AdaptiveMultiVertexFitter. + */ + + +#undef NDEBUG +#include "AdaptiveMultiVertexFitterTestAlg.h" +#include "xAODTracking/Vertex.h" +#include "xAODTracking/TrackParticle.h" +#include "VxMultiVertex/MvfFitInfo.h" +#include "VxMultiVertex/MVFVxTrackAtVertex.h" +#include "VxMultiVertex/TrackToVtxLink.h" +#include "TrkTrack/Track.h" +#include "TrkParticleBase/TrackParticleBase.h" +#include "EventPrimitives/EventPrimitivesHelpers.h" +#include "TestTools/FLOATassert.h" +#include "GaudiKernel/SystemOfUnits.h" +#include <cassert> + + +#include "CLHEP/Vector/LorentzVector.h" + +using Gaudi::Units::mm; +using Gaudi::Units::MeV; +using Gaudi::Units::GeV; + + +namespace { + + +std::unique_ptr<AmgSymMatrix(5)> cov5() +{ + auto m = std::make_unique<AmgSymMatrix(5)>(); + m->setIdentity(); + return m; +} + + +std::unique_ptr<AmgSymMatrix(5)> cov5a() +{ + auto m = std::make_unique<AmgSymMatrix(5)>(); + m->setZero(); + for (int i=0; i < 5; i++) { + (*m)(i,i) = 1e-2; + } + (*m)(1,1)=1; + return m; +} + + +typedef std::vector<std::unique_ptr<Trk::Perigee> > PerigeeUVec_t; +PerigeeUVec_t makePerigees1() +{ + Amg::Vector3D pos0 { 0, 0, 0 }; + + Amg::Vector3D pos1a { 2*mm, 1*mm, -10*mm }; + Amg::Vector3D mom1a { 400*MeV, 600*MeV, 200*MeV }; + Amg::Vector3D pos1b { 1*mm, 2*mm, -3*mm }; + Amg::Vector3D mom1b { 600*MeV, 400*MeV, -200*MeV }; + Amg::Vector3D pos1c { 1.2*mm, 1.3*mm, -7*mm }; + Amg::Vector3D mom1c { 300*MeV, 1000*MeV, 100*MeV }; + + PerigeeUVec_t ret; + + ret.emplace_back (std::make_unique<Trk::Perigee>(pos1a, mom1a, 1, pos1a, cov5().release())); + ret.emplace_back (std::make_unique<Trk::Perigee>(pos1b, mom1b, -1, pos1b, cov5().release())); + ret.emplace_back (std::make_unique<Trk::Perigee>(pos1c, mom1c, -1, pos1c, cov5().release())); + + return ret; +} + + +typedef std::vector<std::unique_ptr<Trk::Perigee> > PerigeeUVec_t; +PerigeeUVec_t makePerigees2() +{ + Amg::Vector3D pos1a { 10*mm, 0*mm, -5*mm }; + Amg::Vector3D mom1a { 1000*MeV, 0*MeV, 0*MeV }; + Amg::Vector3D pos1b { 10.5*mm, 0.5*mm, -5.5*mm }; + Amg::Vector3D mom1b { 800*MeV, 200*MeV, 200*MeV }; + Amg::Vector3D pos1c { 9.5*mm, -0.5*mm, -4.5*mm }; + Amg::Vector3D mom1c { 700*MeV, -300*MeV, -200*MeV }; + + PerigeeUVec_t ret; + + ret.emplace_back (std::make_unique<Trk::Perigee>(pos1a, mom1a, 1, pos1a, cov5a().release())); + ret.emplace_back (std::make_unique<Trk::Perigee>(pos1b, mom1b, -1, pos1a, cov5a().release())); + ret.emplace_back (std::make_unique<Trk::Perigee>(pos1c, mom1c, -1, pos1a, cov5a().release())); + + return ret; +} + + +typedef std::vector<std::unique_ptr<Trk::Track> > TrackUVec_t; +TrackUVec_t makeTracks (const PerigeeUVec_t& perigees) +{ + TrackUVec_t tracks; + + for (const std::unique_ptr<Trk::Perigee>& p : perigees) { + Trk::TrackInfo info (Trk::TrackInfo::Unknown, Trk::undefined); + auto fqual = std::make_unique<Trk::FitQuality> (0, 0); + auto tsos = std::make_unique<DataVector<const Trk::TrackStateOnSurface> >(); + std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern(0); + typePattern.set(Trk::TrackStateOnSurface::Perigee); + tsos->push_back (std::make_unique<Trk::TrackStateOnSurface> + (nullptr, p->clone(), nullptr, nullptr, typePattern)); + tracks.push_back (std::make_unique<Trk::Track> (info, + tsos.release(), + fqual.release())); + } + + return tracks; +} + + +void dumpCovariance (const AmgSymMatrix(5)& m) +{ + for (int i=0; i < 5; i++) { + for (int j=0; j < 5; j++) { + std::cout << m(i,j) << ", "; + } + } +} + + +[[maybe_unused]] +void dumpVertex (const xAOD::Vertex& v) +{ + std::cout << "vvv\n"; + std::cout << v.x() << ", " << v.y() << ", " << v.z() << "\n"; + if (v.isAvailable<short> ("vertexType")) { + std::cout << "vertexType " << v.vertexType() << "\n"; + } + std::cout << "chi2/ndof " << v.chiSquared() << ", " << v.numberDoF() << "\n"; + + std::cout << "cov "; + for (float f : v.covariance()) { + std::cout << f << ", "; + } + std::cout << "\n"; + + if (v.isAvailable<std::vector<ElementLink<xAOD::TrackParticleContainer> > > ("trackParticleLinks")) { + std::cout << "tplinks "; + for (const ElementLink< xAOD::TrackParticleContainer >& l : v.trackParticleLinks()) { + std::cout << l.dataID() << "/" << l.index() << " "; + } + std::cout << "\n"; + } + + if (v.isAvailable<std::vector<float> > ("trackWeights")) { + std::cout << "wt "; + for (float f : v.trackWeights()) { + std::cout << f << " "; + } + std::cout << "\n"; + } + + if (v.isAvailable<std::vector<ElementLink<xAOD::NeutralParticleContainer> > > ("neutralParticleLinks")) { + std::cout << "nplinks "; + for (const ElementLink< xAOD::NeutralParticleContainer >& l : v.neutralParticleLinks()) { + std::cout << l.dataID() << "/" << l.index() << " "; + } + std::cout << "\n"; + } + + if (v.isAvailable<std::vector<float> > ("neutralWeights")) { + std::cout << "wt "; + for (float f : v.neutralWeights()) { + std::cout << f << " "; + } + std::cout << "\n"; + } + + std::cout << v.vxTrackAtVertexAvailable() << "\n"; + for (const Trk::VxTrackAtVertex& vv : v.vxTrackAtVertex()) { + vv.dump (std::cout); + std::cout << "cov "; + if (vv.perigeeAtVertex()) { + dumpCovariance (*vv.perigeeAtVertex()->covariance()); + } + else { + std::cout << "(null)"; + } + std::cout << "\n"; + } +} + + +void assertVec3D (const char* which, + const Amg::Vector3D& a, + const Amg::Vector3D& b, + double thresh = 2e-5) +{ + if ( ! Athena_test::isEqual (a.x(), b.x(), thresh) || + ! Athena_test::isEqual (a.y(), b.y(), thresh) || + ! Athena_test::isEqual (a.z(), b.z(), thresh) ) + { + std::cerr << "TrkVKalVrtFitterTestAlg::assertVec3D mismatch " << which + << ": [" + << a.x() << ", " + << a.y() << ", " + << a.z() << "] / [" + << b.x() << ", " + << b.y() << ", " + << b.z() << "]\n"; + std::abort(); + } +} + + +void compareVertex (const xAOD::Vertex& a, const xAOD::Vertex& b) +{ + assertVec3D ("vertex pos", a.position(), b.position(), 5e-4); + assert (Athena_test::isEqual (a.chiSquared(), b.chiSquared(), 5e-3) ); + assert (Athena_test::isEqual (a.numberDoF(), b.numberDoF(), 5e-3) ); + + assert (a.covariance().size() == b.covariance().size()); + for (unsigned int i = 0; i < a.covariance().size(); i++) { + if (isinf(a.covariance()[i]) && isinf(b.covariance()[i])) continue; + assert (Athena_test::isEqual (a.covariance()[i], b.covariance()[i], 2e-2) ); + } +} + + +void setInitialPerigee (xAOD::Vertex& v, unsigned i, const Trk::Perigee* p) +{ + std::vector< Trk::VxTrackAtVertex >& vec = v.vxTrackAtVertex(); + if (vec.size() <= i) vec.resize(i+1); + vec[i].setInitialPerigee (p); +} + + +void setInitialPerigees (xAOD::Vertex& v, const TrackUVec_t& tracks) +{ + int i = 0; + for (const std::unique_ptr<Trk::Track>& t : tracks) { + setInitialPerigee (v, i, t->perigeeParameters()); + ++i; + } +} + + +struct VertexInfo +{ + xAOD::Vertex v; + TrackUVec_t tracks; + PerigeeUVec_t perigees; + std::unique_ptr<Trk::MvfFitInfo> fi; + std::vector<Trk::VxTrackAtVertex*> vtracks; + std::vector<Trk::TrackToVtxLink> links; +}; + + +void initVertex (VertexInfo& vi, + const Amg::Vector3D& pos, + PerigeeUVec_t&& perigees) +{ + xAOD::Vertex::Decorator< Trk::MvfFitInfo* > MvfFitInfo("MvfFitInfo"); + xAOD::Vertex::Decorator< bool > isInitialized("isInitialized"); + xAOD::Vertex::Decorator< std::vector<Trk::VxTrackAtVertex*> > VTAV("VTAV"); + + vi.perigees = std::move (perigees); + + vi.v.makePrivateStore(); + vi.v.setPosition (pos); + + vi.tracks = makeTracks (vi.perigees); + setInitialPerigees (vi.v, vi.tracks); + + auto cv = std::make_unique<xAOD::Vertex>(); + cv->makePrivateStore(); + cv->setPosition (pos); + Amg::MatrixX looseConstraintCovariance(3, 3); + looseConstraintCovariance.setIdentity(); + looseConstraintCovariance = looseConstraintCovariance * 1e+8; + cv->setCovariancePosition (looseConstraintCovariance); + cv->setFitQuality (0, -3); + vi.fi = std::make_unique<Trk::MvfFitInfo> (cv.release(), + new Amg::Vector3D(pos), + new Amg::Vector3D(pos)); + MvfFitInfo(vi.v) = vi.fi.get(); + isInitialized(vi.v) = false; + + vi.links.reserve (vi.perigees.size()); + for (std::unique_ptr<Trk::Perigee>& p : vi.perigees) { + vi.links.push_back (new std::vector<xAOD::Vertex*> {&vi.v}); + auto tav = std::make_unique<Trk::MVFVxTrackAtVertex>(1.5, p->clone(), p.get()); + tav->setLinkToVertices (&vi.links.back()); + vi.vtracks.push_back (tav.release()); + } + VTAV(vi.v) = vi.vtracks; +} + + +} // anonymous namespace + + +namespace Trk { + + +/** + * @brief Standard Gaudi initialize method. + */ +StatusCode AdaptiveMultiVertexFitterTestAlg::initialize() +{ + ATH_CHECK( m_fitter.retrieve() ); + return StatusCode::SUCCESS; +} + + +/** + * @brief Standard Gaudi execute method. + */ +StatusCode AdaptiveMultiVertexFitterTestAlg::execute() +{ + ATH_MSG_VERBOSE ("execute"); + + ATH_CHECK( test1() ); + + return StatusCode::SUCCESS; +} + + +StatusCode AdaptiveMultiVertexFitterTestAlg::test1() +{ + VertexInfo v1; + initVertex (v1, {1.5*mm, 1.7*mm, -6*mm}, makePerigees1()); + + + VertexInfo v2; + initVertex (v2, {9.8*mm, 0.2*mm, -4.8*mm}, makePerigees2()); + + std::vector<xAOD::Vertex*> verts {&v1.v, &v2.v}; + m_fitter->fit (verts); + + xAOD::Vertex exp_v1; + exp_v1.makePrivateStore(); + exp_v1.setPosition ({8.67978, 10.3569, -6.33368}); + exp_v1.setFitQuality (0.0141012, 0.797808); + exp_v1.setCovariance (std::vector<float> + {751.563, 142.738, 950.953, + -7.62904, 47.5311, 67.4557}); + + xAOD::Vertex exp_v2; + exp_v2.makePrivateStore(); + exp_v2.setPosition ({7.97404, 0.106101, -4.97333}); + exp_v2.setFitQuality (1.74029, 1.6738); + exp_v2.setCovariance (std::vector<float> + {0.348487, -0.0194744, 0.0290909, + -0.00362647, 0.000706365, 0.445415}); + + compareVertex (v1.v, exp_v1); + compareVertex (v2.v, exp_v2); + + return StatusCode::SUCCESS; +} + + +} // namespace Trk diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitterTestAlg.h b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitterTestAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..b8d77a2e620e0b78f19edda1b4f461372d8c030e --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveMultiVertexFitterTestAlg.h @@ -0,0 +1,52 @@ +// This file's extension implies that it's C, but it's really -*- C++ -*-. +/* + * Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration. + */ +/** + * @file TrkVertexFitters/AdaptiveMultiVertexFitterTestAlg.h + * @author scott snyder <snyder@bnl.gov> + * @date Jul, 2019 + * @brief Algorithm for testing AdaptiveMultiVertexFitter. + */ + + +#ifndef TRKVERTEXFITTERS_ADAPTIVEMULTIVERTEXFITTERTESTALG_H +#define TRKVERTEXFITTERS_ADAPTIVEMULTIVERTEXFITTERTESTALG_H + + +#include "AthenaBaseComps/AthAlgorithm.h" +#include "TrkVertexFitters/AdaptiveMultiVertexFitter.h" +#include "GaudiKernel/ToolHandle.h" + + + +namespace Trk { + + +class AdaptiveMultiVertexFitterTestAlg + : public AthAlgorithm +{ +public: + using AthAlgorithm::AthAlgorithm; + + + /// Standard Gaudi initialize method. + virtual StatusCode initialize() override; + + /// Execute the algorithm. + virtual StatusCode execute() override; + + +private: + StatusCode test1(); + + + ToolHandle<Trk::AdaptiveMultiVertexFitter> m_fitter + { this, "Tool", "Trk::AdaptiveMultiVertexFitter", "Tool to test." }; +}; + + +} // namespace Trk + + +#endif // not TRKVERTEXFITTERS_ADAPTIVEMULTIVERTEXFITTERTESTALG_H diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitter.cxx b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitter.cxx index fb80726f59c9201a802af589babe005fbcbb6165..787fd74322c048be8cefcf477fa514252f7521b1 100755 --- a/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitter.cxx +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitter.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration */ @@ -157,25 +157,29 @@ namespace Trk } - Amg::Vector3D* SeedPoint(0); + Amg::Vector3D SeedPoint; //now find the best point for seeding and linearization of the tracks if (IsStartingPoint) { - SeedPoint=new Amg::Vector3D(startingPoint); + SeedPoint=startingPoint; } else { if (perigeeList.size()>1) { if (IsConstraint) { - SeedPoint=new Amg::Vector3D(m_SeedFinder->findSeed(perigeeList,&constraint)); + SeedPoint=m_SeedFinder->findSeed(perigeeList,&constraint); } else { - SeedPoint=new Amg::Vector3D(m_SeedFinder->findSeed(perigeeList)); + SeedPoint=m_SeedFinder->findSeed(perigeeList); } - } else { - SeedPoint=new Amg::Vector3D(constraint.position().x(),constraint.position().y(),constraint.position().z()); + } + else if (IsConstraint) { + SeedPoint=constraint.position(); + } + else { + SeedPoint.setZero(); } } //in case m_onlyzseed is on, just use only the z component given by the seed finder if (m_onlyzseed&&IsConstraint) { - *SeedPoint=Amg::Vector3D(constraint.position().x(),constraint.position().y(),SeedPoint->z()); + SeedPoint=constraint.position(); } @@ -189,7 +193,7 @@ namespace Trk std::vector<const Trk::TrackParameters*>::const_iterator perigeesEnd=perigeeList.end(); ATH_MSG_DEBUG("Inside fitter with track perigee parameters."); - ATH_MSG_DEBUG("Seed point: " << *SeedPoint); + ATH_MSG_DEBUG("Seed point: " << SeedPoint); int myDebugNTrk(0); for (std::vector<const Trk::TrackParameters*>::const_iterator perigeesIter=perigeesBegin;perigeesIter!=perigeesEnd;++perigeesIter) { @@ -203,8 +207,8 @@ namespace Trk VxTrackAtVertex* LinTrackToAdd = new VxTrackAtVertex(0., 0, NULL, (*perigeesIter), NULL); //TODO: Must think now about when to delete this memory! -David S. - //m_LinearizedTrackFactory->linearize(*LinTrackToAdd,*SeedPoint); why linearize it? maybe you don't need it at all!!!!! <19-05-2006> - bool success=m_ImpactPoint3dEstimator->addIP3dAtaPlane(*LinTrackToAdd,*SeedPoint); + //m_LinearizedTrackFactory->linearize(*LinTrackToAdd,SeedPoint); why linearize it? maybe you don't need it at all!!!!! <19-05-2006> + bool success=m_ImpactPoint3dEstimator->addIP3dAtaPlane(*LinTrackToAdd,SeedPoint); if (!success) { msg(MSG::WARNING) << "Adding compatibility to vertex information failed. Newton distance finder didn't converge..." << endmsg; @@ -220,7 +224,7 @@ namespace Trk std::vector<const Trk::NeutralParameters*>::const_iterator neutralPerigeesEnd =neutralPerigeeList.end(); ATH_MSG_DEBUG("Inside fitter with neutral perigee parameters."); - ATH_MSG_DEBUG("Seed point: " << *SeedPoint); + ATH_MSG_DEBUG("Seed point: " << SeedPoint); int myDebugNNeutral(0); for (std::vector<const Trk::NeutralParameters*>::const_iterator neutralPerigeesIter=neutralPerigeesBegin;neutralPerigeesIter!=neutralPerigeesEnd;++neutralPerigeesIter) { @@ -230,7 +234,7 @@ namespace Trk VxTrackAtVertex* LinTrackToAdd = new VxTrackAtVertex(0., 0, NULL, NULL, (*neutralPerigeesIter) ); //TODO: Must think now about when to delete this memory! -David S. - bool success = m_ImpactPoint3dEstimator->addIP3dAtaPlane(*LinTrackToAdd,*SeedPoint); + bool success = m_ImpactPoint3dEstimator->addIP3dAtaPlane(*LinTrackToAdd,SeedPoint); if (!success) { msg(MSG::WARNING) << "Adding compatibility to vertex information failed. Newton distance finder didn't converge..." << endmsg; @@ -253,28 +257,25 @@ namespace Trk - xAOD::Vertex* ConstraintVertex(0); + xAOD::Vertex ConstraintVertex; + ConstraintVertex.makePrivateStore(); //use the previous prior vertex info for the initial vertex if there if (IsConstraint ) { - ConstraintVertex=new xAOD::Vertex(); - ConstraintVertex->makePrivateStore(); - ConstraintVertex->setPosition( constraint.position() ); - ConstraintVertex->setCovariancePosition( constraint.covariancePosition() ); - ConstraintVertex->setFitQuality( 0., 0. ); + ConstraintVertex.setPosition( constraint.position() ); + ConstraintVertex.setCovariancePosition( constraint.covariancePosition() ); + ConstraintVertex.setFitQuality( 0., 0. ); } else { AmgSymMatrix(3) startingCovMatrix; startingCovMatrix.setIdentity(); startingCovMatrix = startingCovMatrix / m_initialError; //now initialize with starting position and covariance matrix - ConstraintVertex=new xAOD::Vertex(); - ConstraintVertex->makePrivateStore(); - ConstraintVertex->setPosition( *SeedPoint ); - ConstraintVertex->setCovariancePosition( startingCovMatrix ); - ConstraintVertex->setFitQuality( 0., -3. ); + ConstraintVertex.setPosition( SeedPoint ); + ConstraintVertex.setCovariancePosition( startingCovMatrix ); + ConstraintVertex.setFitQuality( 0., -3. ); } //now put the linearizedtracks into VxTrackAtVertex - return dothefit(*ConstraintVertex,*SeedPoint,theLinTracks); + return dothefit(ConstraintVertex,SeedPoint,theLinTracks); } @@ -599,9 +600,6 @@ namespace Trk << "the ndf of the vertex is at fit end: " << ActualVertex->numberDoF() << endmsg; } - delete &ConstraintVertex; - delete &SeedVertex; - // TODO: get rid of following line //std::cout << "Number of steps: " << num_steps << ". Number of relinearizations: " << num_relinearizations << "." << std::endl << std::endl; diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitterTestAlg.cxx b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitterTestAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b1a69f7731d573859c57cbbae002d6eded7588b6 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitterTestAlg.cxx @@ -0,0 +1,755 @@ +/* + * Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration. + */ +/** + * @file TrkVertexFitters/src/AdaptiveVertexFitterTestAlg.cxx + * @author scott snyder <snyder@bnl.gov> + * @date Jul, 2019 + * @brief Algorithm for testing AdaptiveVertexFitter. + */ + + +#undef NDEBUG +#include "AdaptiveVertexFitterTestAlg.h" +#include "xAODTracking/Vertex.h" +#include "xAODTracking/TrackParticle.h" +#include "TrkTrack/Track.h" +#include "TrkParticleBase/TrackParticleBase.h" +#include "EventPrimitives/EventPrimitivesHelpers.h" +#include "TestTools/FLOATassert.h" +#include "GaudiKernel/SystemOfUnits.h" +#include <cassert> + + +#include "CLHEP/Vector/LorentzVector.h" + +using Gaudi::Units::mm; +using Gaudi::Units::MeV; +using Gaudi::Units::GeV; + + +namespace { + + +template <class T> +std::vector<const T*> asVec (const std::vector<std::unique_ptr<T> >& v) +{ + std::vector<const T*> ret; + for (const std::unique_ptr<T>& p : v) { + ret.push_back (p.get()); + } + return ret; +} + + +std::vector<const Trk::TrackParameters*> +asVec (const std::vector<std::unique_ptr<Trk::Perigee> >& v) +{ + std::vector<const Trk::TrackParameters*> ret; + for (const std::unique_ptr<Trk::Perigee>& p : v) { + ret.push_back (p.get()); + } + return ret; +} + + +std::vector<const Trk::NeutralParameters*> +asVec (const std::vector<std::unique_ptr<Trk::NeutralPerigee> >& v) +{ + std::vector<const Trk::NeutralParameters*> ret; + for (const std::unique_ptr<Trk::NeutralPerigee>& p : v) { + ret.push_back (p.get()); + } + return ret; +} + + +std::unique_ptr<AmgSymMatrix(5)> cov5() +{ + auto m = std::make_unique<AmgSymMatrix(5)>(); + m->setIdentity(); + return m; +} + + +typedef std::vector<std::unique_ptr<Trk::Perigee> > PerigeeUVec_t; +PerigeeUVec_t makePerigees1() +{ + Amg::Vector3D pos0 { 0, 0, 0 }; + + Amg::Vector3D pos1a { 2*mm, 1*mm, -10*mm }; + Amg::Vector3D mom1a { 400*MeV, 600*MeV, 200*MeV }; + Amg::Vector3D pos1b { 1*mm, 2*mm, -3*mm }; + Amg::Vector3D mom1b { 600*MeV, 400*MeV, -200*MeV }; + Amg::Vector3D pos1c { 1.2*mm, 1.3*mm, -7*mm }; + Amg::Vector3D mom1c { 300*MeV, 1000*MeV, 100*MeV }; + + PerigeeUVec_t ret; + + ret.emplace_back (std::make_unique<Trk::Perigee>(pos1a, mom1a, 1, pos1a, cov5().release())); + ret.emplace_back (std::make_unique<Trk::Perigee>(pos1b, mom1b, -1, pos1b, cov5().release())); + ret.emplace_back (std::make_unique<Trk::Perigee>(pos1c, mom1c, -1, pos1c, cov5().release())); + + return ret; +} + + +typedef std::vector<std::unique_ptr<Trk::NeutralPerigee> > NeutralUVec_t; +NeutralUVec_t makeNeutrals1() +{ + Amg::Vector3D pos0 { 0, 0, 0 }; + + Amg::Vector3D pos1a { 3*mm, 0.5*mm, -7*mm }; + Amg::Vector3D mom1a { 1000*MeV, 900*MeV, 2000*MeV }; + Amg::Vector3D pos1b { -1*mm, 2.5*mm, 1*mm }; + Amg::Vector3D mom1b { 800*MeV, 1000*MeV, 300*MeV }; + Amg::Vector3D pos1c { -1.5*mm, 1*mm, -3*mm }; + Amg::Vector3D mom1c { 500*MeV, 4000*MeV, 800*MeV }; + + NeutralUVec_t ret; + + ret.emplace_back (std::make_unique<Trk::NeutralPerigee>(pos1a, mom1a, 1, pos1a, cov5().release())); + ret.emplace_back (std::make_unique<Trk::NeutralPerigee>(pos1b, mom1b, 1, pos1b, cov5().release())); + ret.emplace_back (std::make_unique<Trk::NeutralPerigee>(pos1c, mom1c, 1, pos1c, cov5().release())); + + return ret; +} + + +typedef std::vector<std::unique_ptr<Trk::Track> > TrackUVec_t; +TrackUVec_t makeTracks (PerigeeUVec_t&& perigees) +{ + TrackUVec_t tracks; + + for (std::unique_ptr<Trk::Perigee>& p : perigees) { + Trk::TrackInfo info (Trk::TrackInfo::Unknown, Trk::undefined); + auto fqual = std::make_unique<Trk::FitQuality> (0, 0); + auto tsos = std::make_unique<DataVector<const Trk::TrackStateOnSurface> >(); + std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern(0); + typePattern.set(Trk::TrackStateOnSurface::Perigee); + tsos->push_back (std::make_unique<Trk::TrackStateOnSurface> + (nullptr, p.release(), nullptr, nullptr, typePattern)); + tracks.push_back (std::make_unique<Trk::Track> (info, + tsos.release(), + fqual.release())); + } + + return tracks; +} + + +typedef std::vector<std::unique_ptr<Trk::TrackParticleBase> > TrackParticleUVec_t; +TrackParticleUVec_t makeTrackParticles (PerigeeUVec_t&& perigees) +{ + TrackParticleUVec_t tracks; + + std::vector<const Trk::TrackParameters*> pvec; + for (std::unique_ptr<Trk::Perigee>& p : perigees) { + tracks.push_back (std::make_unique<Trk::TrackParticleBase> (nullptr, + Trk::NoVtx, + nullptr, + nullptr, + pvec, + p.release(), + nullptr)); + } + + return tracks; +} + + +// Copied from TrackParticleCreatorTool. +void setDefiningParameters( xAOD::TrackParticle& tp, + const Trk::Perigee& perigee ) +{ + tp.setDefiningParameters(perigee.parameters()[Trk::d0], + perigee.parameters()[Trk::z0], + perigee.parameters()[Trk::phi0], + perigee.parameters()[Trk::theta], + perigee.parameters()[Trk::qOverP]); + const AmgSymMatrix(5)* covMatrix = perigee.covariance(); + // see https://its.cern.ch/jira/browse/ATLASRECTS-645 for justification to comment out the following line + // assert(covMatrix && covMatrix->rows()==5&& covMatrix->cols()==5); + std::vector<float> covMatrixVec; + if( !covMatrix ) ;//ATH_MSG_WARNING("Setting Defining parameters without error matrix"); + else Amg::compress(*covMatrix,covMatrixVec); + tp.setDefiningParametersCovMatrixVec(covMatrixVec); + const Amg::Vector3D& surfaceCenter = perigee.associatedSurface().center(); + tp.setParametersOrigin(surfaceCenter.x(), surfaceCenter.y(), surfaceCenter.z() ); +} +void setDefiningParameters( xAOD::NeutralParticle& tp, + const Trk::NeutralPerigee& perigee ) +{ + tp.setDefiningParameters(perigee.parameters()[Trk::d0], + perigee.parameters()[Trk::z0], + perigee.parameters()[Trk::phi0], + perigee.parameters()[Trk::theta], + perigee.parameters()[Trk::qOverP]); + const AmgSymMatrix(5)* covMatrix = perigee.covariance(); + // see https://its.cern.ch/jira/browse/ATLASRECTS-645 for justification to comment out the following line + // assert(covMatrix && covMatrix->rows()==5&& covMatrix->cols()==5); + std::vector<float> covMatrixVec; + if( !covMatrix ) ;//ATH_MSG_WARNING("Setting Defining parameters without error matrix"); + else Amg::compress(*covMatrix,covMatrixVec); + tp.setDefiningParametersCovMatrixVec(covMatrixVec); + const Amg::Vector3D& surfaceCenter = perigee.associatedSurface().center(); + tp.setParametersOrigin(surfaceCenter.x(), surfaceCenter.y(), surfaceCenter.z() ); +} + + +typedef std::vector<std::unique_ptr<xAOD::TrackParticle> > xAODTPUVec_t; +xAODTPUVec_t makexAODTP (PerigeeUVec_t&& perigees) +{ + xAODTPUVec_t tracks; + + for (std::unique_ptr<Trk::Perigee>& p : perigees) { + auto tp = std::make_unique<xAOD::TrackParticle>(); + tp->makePrivateStore(); + setDefiningParameters (*tp, *p); + tracks.push_back (std::move (tp)); + } + + return tracks; +} + + +typedef std::vector<std::unique_ptr<xAOD::NeutralParticle> > xAODNPUVec_t; +xAODNPUVec_t makexAODNP (NeutralUVec_t&& perigees) +{ + xAODNPUVec_t tracks; + + for (std::unique_ptr<Trk::NeutralPerigee>& p : perigees) { + auto tp = std::make_unique<xAOD::NeutralParticle>(); + tp->makePrivateStore(); + setDefiningParameters (*tp, *p); + tracks.push_back (std::move (tp)); + } + + return tracks; +} + + +void dumpCovariance (const AmgSymMatrix(5)& m) +{ + for (int i=0; i < 5; i++) { + for (int j=0; j < 5; j++) { + std::cout << m(i,j) << ", "; + } + } +} + + +[[maybe_unused]] +void dumpVertex (const xAOD::Vertex& v) +{ + std::cout << "vvv\n"; + std::cout << v.x() << ", " << v.y() << ", " << v.z() << "\n"; + if (v.isAvailable<short> ("vertexType")) { + std::cout << "vertexType " << v.vertexType() << "\n"; + } + std::cout << "chi2/ndof " << v.chiSquared() << ", " << v.numberDoF() << "\n"; + + std::cout << "cov "; + for (float f : v.covariance()) { + std::cout << f << ", "; + } + std::cout << "\n"; + + if (v.isAvailable<std::vector<ElementLink<xAOD::TrackParticleContainer> > > ("trackParticleLinks")) { + std::cout << "tplinks "; + for (const ElementLink< xAOD::TrackParticleContainer >& l : v.trackParticleLinks()) { + std::cout << l.dataID() << "/" << l.index() << " "; + } + std::cout << "\n"; + } + + if (v.isAvailable<std::vector<float> > ("trackWeights")) { + std::cout << "wt "; + for (float f : v.trackWeights()) { + std::cout << f << " "; + } + std::cout << "\n"; + } + + if (v.isAvailable<std::vector<ElementLink<xAOD::NeutralParticleContainer> > > ("neutralParticleLinks")) { + std::cout << "nplinks "; + for (const ElementLink< xAOD::NeutralParticleContainer >& l : v.neutralParticleLinks()) { + std::cout << l.dataID() << "/" << l.index() << " "; + } + std::cout << "\n"; + } + + if (v.isAvailable<std::vector<float> > ("neutralWeights")) { + std::cout << "wt "; + for (float f : v.neutralWeights()) { + std::cout << f << " "; + } + std::cout << "\n"; + } + + std::cout << v.vxTrackAtVertexAvailable() << "\n"; + for (const Trk::VxTrackAtVertex& vv : v.vxTrackAtVertex()) { + vv.dump (std::cout); + std::cout << "cov "; + if (vv.perigeeAtVertex()) { + dumpCovariance (*vv.perigeeAtVertex()->covariance()); + } + else { + std::cout << "(null)"; + } + std::cout << "\n"; + } +} + + +void assertVec3D (const char* which, + const Amg::Vector3D& a, + const Amg::Vector3D& b) +{ + if ( ! Athena_test::isEqual (a.x(), b.x(), 2e-5) || + ! Athena_test::isEqual (a.y(), b.y(), 2e-5) || + ! Athena_test::isEqual (a.z(), b.z(), 2e-5) ) + { + std::cerr << "TrkVKalVrtFitterTestAlg::assertVec3D mismatch " << which + << ": [" + << a.x() << ", " + << a.y() << ", " + << a.z() << "] / [" + << b.x() << ", " + << b.y() << ", " + << b.z() << "]\n"; + std::abort(); + } +} + + +void comparePerigee (const Trk::TrackParameters* a, + const Trk::TrackParameters* b) +{ + if (!a && !b) return; + if (!a || !b) std::abort(); + assertVec3D ("perigee pos", a->position(), b->position()); + assertVec3D ("perigee mom", a->momentum(), b->momentum()); + assert (a->charge() == b->charge()); + assertVec3D ("perigee surf", + a->associatedSurface().center(), + b->associatedSurface().center()); +} + + +void compareVertex (const xAOD::Vertex& a, const xAOD::Vertex& b) +{ + assertVec3D ("vertex pos", a.position(), b.position()); + assert (Athena_test::isEqual (a.chiSquared(), b.chiSquared(), 5e-5) ); + assert (Athena_test::isEqual (a.numberDoF(), b.numberDoF(), 1e-5) ); + + assert (a.covariance().size() == b.covariance().size()); + for (unsigned int i = 0; i < a.covariance().size(); i++) { + if (isinf(a.covariance()[i]) && isinf(b.covariance()[i])) continue; + assert (Athena_test::isEqual (a.covariance()[i], b.covariance()[i], 2e-2) ); + } + + assert (a.vxTrackAtVertexAvailable() == b.vxTrackAtVertexAvailable()); + if (a.vxTrackAtVertexAvailable()) { + const std::vector< Trk::VxTrackAtVertex >& avec = a.vxTrackAtVertex(); + const std::vector< Trk::VxTrackAtVertex >& bvec = b.vxTrackAtVertex(); + assert (avec.size() == bvec.size()); + for (unsigned int i = 0; i < avec.size(); i++) { + comparePerigee (avec[i].initialPerigee(), bvec[i].initialPerigee()); + comparePerigee (avec[i].perigeeAtVertex(), bvec[i].perigeeAtVertex()); + assert (Athena_test::isEqual (avec[i].trackQuality().chiSquared(), + bvec[i].trackQuality().chiSquared(), + 3e-2)); + assert (avec[i].trackQuality().numberDoF() == + bvec[i].trackQuality().numberDoF()); + } + } +} + + +void setInitialPerigee (xAOD::Vertex& v, unsigned i, const Trk::Perigee* p) +{ + std::vector< Trk::VxTrackAtVertex >& vec = v.vxTrackAtVertex(); + if (vec.size() <= i) vec.resize(i+1); + vec[i].setInitialPerigee (p); +} + + +void setInitialPerigees (xAOD::Vertex& v, const TrackUVec_t& tracks) +{ + int i = 0; + for (const std::unique_ptr<Trk::Track>& t : tracks) { + setInitialPerigee (v, i, t->perigeeParameters()); + ++i; + } +} + + +void clearInitialPerigees (xAOD::Vertex& v) +{ + for (Trk::VxTrackAtVertex& v : v.vxTrackAtVertex()) { + v.setInitialPerigee (static_cast<const Trk::Perigee*>(nullptr)); + } +} + + +void setRefittedPerigee (xAOD::Vertex& v, unsigned i, + float charge, + const Amg::Vector3D& pos, + const Amg::Vector3D& mom, + const std::vector<float>& c) +{ + std::vector< Trk::VxTrackAtVertex >& vec = v.vxTrackAtVertex(); + if (vec.size() <= i) vec.resize(i+1); + + std::unique_ptr<AmgSymMatrix(5)> cov = cov5(); + for (int i=0; i < 5; i++) { + for (int j=0; j < 5; j++) { + unsigned ipos = i*5 + j; + (*cov)(i,j) = ipos < c.size() ? c[ipos] : 0; + } + } + + const Amg::Vector3D& vpos = v.position(); + auto p = std::make_unique<Trk::Perigee> (pos, mom, charge, vpos, + cov.release()); + vec[i].setPerigeeAtVertex (p.release()); +} + + +void setFitQuality (xAOD::Vertex& v, unsigned i, float chi2, int ndof) +{ + std::vector< Trk::VxTrackAtVertex >& vec = v.vxTrackAtVertex(); + if (vec.size() <= i) vec.resize(i+1); + vec[i].setTrackQuality (Trk::FitQuality (chi2, ndof)); +} + + +} // anonymous namespace + + +namespace Trk { + + +/** + * @brief Standard Gaudi initialize method. + */ +StatusCode AdaptiveVertexFitterTestAlg::initialize() +{ + ATH_CHECK( m_fitter.retrieve() ); + return StatusCode::SUCCESS; +} + + +/** + * @brief Standard Gaudi execute method. + */ +StatusCode AdaptiveVertexFitterTestAlg::execute() +{ + ATH_MSG_VERBOSE ("execute"); + + ATH_CHECK( test1() ); + ATH_CHECK( test2() ); + ATH_CHECK( test3() ); + ATH_CHECK( test4() ); + ATH_CHECK( test5() ); + ATH_CHECK( test6() ); + + return StatusCode::SUCCESS; +} + + +// Charged, no constraint. +StatusCode AdaptiveVertexFitterTestAlg::test1() +{ + xAOD::Vertex exp_v0; + exp_v0.makePrivateStore(); + exp_v0.setPosition ({8.43429, 9.21288, -6.60683}); + exp_v0.setFitQuality (0.235948, 1.83137); + exp_v0.setCovariance (std::vector<float> + {35.7266, 30.7883, 47.0715, 3.14669, 4.52157, 6.74599 }); + setRefittedPerigee (exp_v0, 0, 1, + {7.8113830, 9.6359772, -7.1130181}, + {405.1715584, 596.5212045, 199.9961520}, + {1239.56, -45.9167, 227.862, -0.121956, -33.6206, -45.9167, 128.234, -8.63731, -11.1976, 1.30089, 227.862, -8.63731, 42.8345, -0.0234628, -6.46787, -0.121956, -11.1976, -0.0234628, 1.00001, 0.00399819, -33.6206, 1.30089, -6.46787, 0.00399819, 1} + ); + setFitQuality (exp_v0, 0, 0, 2); + setRefittedPerigee (exp_v0, 1, -1, + {9.4647993, 7.6951354, -5.8295872}, + {596.5909112, 405.0689507, -199.9960745}, + {1150.81, -107.955, 215.525, -0.120009, -32.3376, -107.955, 132.06, -20.6979, -10.9741, 3.17765, 215.525, -20.6979, 41.3112, -0.0235635, -6.34888, -0.120009, -10.9741, -0.0235635, 1.00001, 0.00408136, -32.3376, 3.17765, -6.34888, 0.00408136, 1} + ); + setFitQuality (exp_v0, 1, 0.175, 2); + setRefittedPerigee (exp_v0, 2, -1, + {3.9422703, 10.5334868, -6.0773998}, + {294.4712678, 1001.6416761, 100.0031255}, + {872.333, -76.5329, 170.831, -0.0897349, -27.9026, -76.5329, 102.388, -15.3986, -9.71236, 2.58803, 170.831, -15.3986, 34.3815, -0.0185874, -5.77814, -0.0897349, -9.71236, -0.0185874, 1.00001, 0.00312552, -27.9026, 2.58803, -5.77814, 0.00312552, 1} + ); + setFitQuality (exp_v0, 2, 0.027, 2); + + TrackUVec_t tracks = makeTracks (makePerigees1()); + std::unique_ptr<xAOD::Vertex> v1 (m_fitter->fit (asVec (tracks))); + setInitialPerigees (exp_v0, tracks); + compareVertex (*v1, exp_v0); + + PerigeeUVec_t perigees = makePerigees1(); + std::unique_ptr<xAOD::Vertex> v2 (m_fitter->fit (asVec (perigees))); + compareVertex (*v2, exp_v0); + + return StatusCode::SUCCESS; +} + + +// Neutral, no constraint. +StatusCode AdaptiveVertexFitterTestAlg::test2() +{ + xAOD::Vertex exp_v0; + exp_v0.makePrivateStore(); + exp_v0.setPosition ({4.27611, 4.35683, 1.62467}); + exp_v0.setFitQuality (1.69878, 0.79376); + exp_v0.setCovariance (std::vector<float> + {28.318, 26.913, 37.8531, 17.3323, 17.3041, 23.6152}); + setFitQuality (exp_v0, 0, 0.936, 1); + setFitQuality (exp_v0, 1, 0.000, 2); + setFitQuality (exp_v0, 2, 0.000, 2); + + NeutralUVec_t neutrals = makeNeutrals1(); + std::unique_ptr<xAOD::Vertex> v1 (m_fitter->fit (std::vector<const Trk::TrackParameters*>(), + asVec (neutrals))); + compareVertex (*v1, exp_v0); + + return StatusCode::SUCCESS; +} + + +// Charged + Vector3D constraint. +StatusCode AdaptiveVertexFitterTestAlg::test3() +{ + xAOD::Vertex exp_v0; + exp_v0.makePrivateStore(); + exp_v0.setPosition ({8.82019, 9.90325, -6.4356}); + exp_v0.setFitQuality (0.0430536, 1.85853); + exp_v0.setCovariance (std::vector<float> + {198.998, 111.025, 256.182, 7.89219, 18.6016, 25.5699}); + setRefittedPerigee (exp_v0, 0, 1, + {8.2545569, 10.2879887, -6.8943690}, + {405.5620128, 596.2559088, 199.9958625}, + {1640.25, -51.1682, 281.924, -0.152864, -38.9104, -51.1682, 147.868, -8.97462, -12.0465, 1.26021, 281.924, -8.97462, 49.4107, -0.0273351, -6.9577, -0.152864, -12.0465, -0.0273351, 1.00002, 0.00429935, -38.9104, 1.26021, -6.9577, 0.00429935, 1} + ); + setFitQuality (exp_v0, 0, 0.000, 2); + setRefittedPerigee (exp_v0, 1, -1, + {10.0507183, 8.0932102, -6.0260413}, + {596.3526270, 405.4198090, -199.9958068}, + {1487.84, -156.737, 262.111, -0.147597, -36.9838, -156.737, 155.885, -28.189, -11.7322, 4.05867, 262.111, -28.189, 47.1302, -0.0271073, -6.79166, -0.147597, -11.7322, -0.0271073, 1.00002, 0.00436106, -36.9838, 4.05867, -6.79166, 0.00436106, 1 } + ); + setFitQuality (exp_v0, 1, 0.017, 2); + setRefittedPerigee (exp_v0, 2, -1, + {4.1590180, 11.2713536, -6.0037370}, + {294.0294524, 1001.7714343, 100.0033730}, + {1167.78, -99.1794, 213.359, -0.112662, -32.5361, -99.1794, 119.804, -18.5507, -10.4868, 2.89953, 213.359, -18.5507, 39.9175, -0.0216083, -6.23884, -0.112662, -10.4868, -0.0216083, 1.00001, 0.00337217, -32.5361, 2.89953, -6.23884, 0.00337217, 1 } + ); + setFitQuality (exp_v0, 2, 0.020, 2); + + Amg::Vector3D pnt1 (5, 6, -3); + + TrackUVec_t tracks = makeTracks (makePerigees1()); + std::unique_ptr<xAOD::Vertex> v1 (m_fitter->fit (asVec (tracks), pnt1)); + setInitialPerigees (exp_v0, tracks); + compareVertex (*v1, exp_v0); + + PerigeeUVec_t perigees = makePerigees1(); + std::unique_ptr<xAOD::Vertex> v2 (m_fitter->fit (asVec (perigees), pnt1)); + compareVertex (*v2, exp_v0); + + TrackParticleUVec_t tps = makeTrackParticles (makePerigees1()); + std::unique_ptr<xAOD::Vertex> v3 (m_fitter->fit (asVec (tps), pnt1)); + compareVertex (*v3, exp_v0); + + xAODTPUVec_t xtps = makexAODTP (makePerigees1()); + std::unique_ptr<xAOD::Vertex> v4 (m_fitter->fit (asVec (xtps), pnt1)); + clearInitialPerigees (exp_v0); + compareVertex (*v4, exp_v0); + + return StatusCode::SUCCESS; +} + + +// Charged + Vertex constraint. +StatusCode AdaptiveVertexFitterTestAlg::test4() +{ + xAOD::Vertex exp_v0; + exp_v0.makePrivateStore(); + exp_v0.setPosition ({5.00265, 6.02993, -3.2627}); + exp_v0.setFitQuality (0.904414, 4.84469); + exp_v0.setCovariance (std::vector<float> + {1.24242, 0.00898474, 1.25033, + -0.0018625, 0.00503271, 1.16608}); + setRefittedPerigee (exp_v0, 0, 1, + {5.2563912, 5.8589792, -8.3777169}, + {402.9097324, 598.0506755, 199.9978428}, + { 148.334, 3.90191, 44.5024, -0.019887, -10.635, 3.90191, 40.8102, 1.24424, -6.30148, -0.32025, 44.5024, 1.24424, 14.2092, -0.00679499, -3.63453, -0.019887, -6.30148, -0.00679499, 1, 0.00224026, -10.635, -0.32025, -3.63453, 0.00224026, 1 } + ); + setFitQuality (exp_v0, 0, 0.682, 2); + setRefittedPerigee (exp_v0, 1, -1, + {5.6281930, 5.1009479, -4.5451087}, + {598.1437842, 402.7714914, -199.9978479}, + {124.605, -11.9636, 38.8859, -0.0179296, -9.62427, -11.9636, 38.1524, -3.99608, -5.99156, 1.06341, 38.8859, -3.99608, 12.9843, -0.00645023, -3.46157, -0.0179296, -5.99156, -0.00645023, 1, 0.00223385, -9.62427, 1.06341, -3.46157, 0.00223385, 1 } + ); + setFitQuality (exp_v0, 1, 0.061, 2); + setRefittedPerigee (exp_v0, 2, -1, + {2.8050119, 6.6815112, -6.4621037}, + {296.7777220, 1000.9608467, 100.0018360}, + { 123.052, -8.25619, 37.7378, -0.0183804, -9.50958, -8.25619, 33.7962, -2.7124, -5.67336, 0.738344, 37.7378, -2.7124, 12.4073, -0.0065315, -3.37786, -0.0183804, -5.67336, -0.0065315, 1, 0.0018425, -9.50958, 0.738344, -3.37786, 0.0018425, 1 } + ); + setFitQuality (exp_v0, 2, 0.352, 2); + + Amg::Vector3D pnt1 (5, 6, -3); + xAOD::Vertex pnt2; + pnt2.makePrivateStore(); + pnt2.setPosition (pnt1); + AmgSymMatrix(3) pnt2covar; + pnt2covar.setIdentity(); + pnt2.setCovariancePosition (pnt2covar); + + TrackUVec_t tracks = makeTracks (makePerigees1()); + std::unique_ptr<xAOD::Vertex> v1 (m_fitter->fit (asVec (tracks), pnt2)); + setInitialPerigees (exp_v0, tracks); + compareVertex (*v1, exp_v0); + + PerigeeUVec_t perigees = makePerigees1(); + std::unique_ptr<xAOD::Vertex> v2 (m_fitter->fit (asVec (perigees), pnt2)); + compareVertex (*v2, exp_v0); + + TrackParticleUVec_t tps = makeTrackParticles (makePerigees1()); + std::unique_ptr<xAOD::Vertex> v3 (m_fitter->fit (asVec (tps), pnt2)); + compareVertex (*v3, exp_v0); + + xAODTPUVec_t xtps = makexAODTP (makePerigees1()); + std::unique_ptr<xAOD::Vertex> v4 (m_fitter->fit (asVec (xtps), pnt2)); + clearInitialPerigees (exp_v0); + compareVertex (*v4, exp_v0); + + return StatusCode::SUCCESS; +} + + +// Charged + Neutral + Vector3D constraint +StatusCode AdaptiveVertexFitterTestAlg::test5() +{ + xAOD::Vertex exp_v0; + exp_v0.makePrivateStore(); + exp_v0.setPosition ({-0.664619, 1.88247, -4.28452}); + exp_v0.setFitQuality (1.36786, 6.60783); + exp_v0.setCovariance (std::vector<float> + {25.8826, 26.6122, 91.0458, 6.34189, 14.6174, 13.9884}); + setRefittedPerigee (exp_v0, 0, 1, + {1.5885284, 0.3823804, -10.2058300}, + {399.6301520, 600.2463141, 200.0002601}, + { 1.58278, -0.618193, -0.821391, 0.00010879, -0.171884, -0.618193, 2.32566, 0.908938, 0.800925, 0.346876, -0.821391, 0.908938, 1.20816, -0.000294073, 0.461137, 0.00010879, 0.800925, -0.000294073, 1, -0.000269915, -0.171884, 0.346876, 0.461137, -0.000269915, 1 } + ); + setFitQuality (exp_v0, 0, 0.000, 2); + setRefittedPerigee (exp_v0, 1, -1, + {-0.2079830, 1.1957289, -2.5975002}, + {600.4814296, 399.2766328, -200.0005578}, + { 3.53014, 0.466334, -2.04043, 0.000621337, -0.653413, 0.466334, 3.53405, -0.415368, 1.56192, -0.206492, -2.04043, -0.415368, 1.81448, -0.000856625, 0.901728, 0.000621337, 1.56192, -0.000856625, 1, -0.000578862, -0.653413, -0.206492, 0.901728, -0.000578862, 1 } + ); + setFitQuality (exp_v0, 1, 0.017, 2); + setRefittedPerigee (exp_v0, 2, -1, + {1.2063615, 1.3212056, -6.9978818}, + {299.9873170, 1000.0038041, 100.0000072}, + { 1.00049, 0.00413671, 0.0221412, -2.7918e-08, -0.000147081, 0.00413671, 1.03551, 0.18734, -0.0223173, -0.00248881, 0.0221412, 0.18734, 1.00242, -1.91988e-06, -0.0133167, -2.7918e-08, -0.0223173, -1.91988e-06, 1, 7.24261e-06, -0.000147081, -0.00248881, -0.0133167, 7.24261e-06, 1 } + ); + setFitQuality (exp_v0, 2, 0.020, 2); + setFitQuality (exp_v0, 3, 0.136, 2); + setFitQuality (exp_v0, 4, 0.000, 2); + setFitQuality (exp_v0, 5, 0.000, 2); + + Amg::Vector3D pnt1 (5, 6, -3); + + TrackUVec_t tracks = makeTracks (makePerigees1()); + setInitialPerigees (exp_v0, tracks); + + PerigeeUVec_t perigees = makePerigees1(); + NeutralUVec_t neutrals = makeNeutrals1(); + std::unique_ptr<xAOD::Vertex> v1 (m_fitter->fit (asVec (perigees), + asVec (neutrals), + pnt1)); + compareVertex (*v1, exp_v0); + + xAODTPUVec_t xtps = makexAODTP (makePerigees1()); + xAODNPUVec_t xaodnp = makexAODNP (makeNeutrals1()); + std::unique_ptr<xAOD::Vertex> v2 (m_fitter->fit (asVec (xtps), + asVec (xaodnp), + pnt1)); + clearInitialPerigees (exp_v0); + compareVertex (*v2, exp_v0); + + return StatusCode::SUCCESS; +} + + +// Charged+Neutral + Vertex constraint +StatusCode AdaptiveVertexFitterTestAlg::test6() +{ + xAOD::Vertex exp_v0; + exp_v0.makePrivateStore(); + exp_v0.setPosition ({4.85215, 5.94893, -3.13472}); + exp_v0.setFitQuality (2.38431, 8.54327); + exp_v0.setCovariance (std::vector<float> + {1.183, 0.0323074, 1.21271, + 0.00903037, 0.0167373, 1.12584}); + setRefittedPerigee (exp_v0, 0, 1, + {5.1719083, 5.7335618, -8.4196568}, + {402.8346276, 598.1012479, 199.9979000}, + { 135.368, 4.54292, 41.4349, -0.0182748, -10.0936, 4.54292, 38.8427, 1.48244, -6.13913, -0.389733, 41.4349, 1.48244, 13.5349, -0.00640884, -3.54057, -0.0182748, -6.13913, -0.00640884, 1, 0.00218082, -10.0936, -0.389733, -3.54057, 0.00218082, 1 } + ); + setFitQuality (exp_v0, 0, 0.682, 2); + + setRefittedPerigee (exp_v0, 1, -1, + {5.4869777, 5.0058724, -4.4978936}, + {598.2006961, 402.6869274, -199.9979142}, + { 111.922, -11.1715, 35.7532, -0.0162277, -9.04482, -11.1715, 35.9577, -3.83303, -5.80842, 1.04702, 35.7532, -3.83303, 12.2632, -0.0060222, -3.35579, -0.0162277, -5.80842, -0.0060222, 1, 0.00216502, -9.04482, 1.04702, -3.35579, 0.00216502, 1 } + ); + setFitQuality (exp_v0, 1, 0.061, 2); + + setRefittedPerigee (exp_v0, 2, -1, + {2.7708142, 6.5661849, -6.4736255}, + {296.8467752, 1000.9403742, 100.0017963}, + { 114.181, -7.37341, 35.593, -0.0172367, -9.10529, -7.37341, 32.345, -2.46816, -5.55155, 0.684107, 35.593, -2.46816, 11.9239, -0.00626001, -3.30551, -0.0172367, -5.55155, -0.00626001, 1, 0.00180269, -9.10529, 0.684107, -3.30551, 0.00180269, 1 } + ); + setFitQuality (exp_v0, 2, 0.352, 2); + setFitQuality (exp_v0, 3, 1.119, 1); + setFitQuality (exp_v0, 4, 0.000, 2); + setFitQuality (exp_v0, 5, 0.000, 2); + + Amg::Vector3D pnt1 (5, 6, -3); + xAOD::Vertex pnt2; + pnt2.makePrivateStore(); + pnt2.setPosition (pnt1); + AmgSymMatrix(3) pnt2covar; + pnt2covar.setIdentity(); + pnt2.setCovariancePosition (pnt2covar); + + TrackUVec_t tracks = makeTracks (makePerigees1()); + setInitialPerigees (exp_v0, tracks); + + PerigeeUVec_t perigees = makePerigees1(); + NeutralUVec_t neutrals = makeNeutrals1(); + std::unique_ptr<xAOD::Vertex> v1 (m_fitter->fit (asVec (perigees), + asVec (neutrals), + pnt2)); + compareVertex (*v1, exp_v0); + + xAODTPUVec_t xtps = makexAODTP (makePerigees1()); + xAODNPUVec_t xaodnp = makexAODNP (makeNeutrals1()); + std::unique_ptr<xAOD::Vertex> v2 (m_fitter->fit (asVec (xtps), + asVec (xaodnp), + pnt2)); + clearInitialPerigees (exp_v0); + compareVertex (*v2, exp_v0); + + return StatusCode::SUCCESS; +} + + +} // namespace Trk diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitterTestAlg.h b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitterTestAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..359d24b3728c87d7bfbeffe2a4d9f73d19f8d618 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/src/AdaptiveVertexFitterTestAlg.h @@ -0,0 +1,57 @@ +// This file's extension implies that it's C, but it's really -*- C++ -*-. +/* + * Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration. + */ +/** + * @file TrkVertexFitters/src/AdaptiveVertexFitterTestAlg.h + * @author scott snyder <snyder@bnl.gov> + * @date Jul, 2019 + * @brief Algorithm for testing AdaptiveVertexFitter. + */ + + +#ifndef TRKVERTEXFITTERS_ADAPTIVEVERTEXFITTERTESTALG_H +#define TRKVERTEXFITTERS_ADAPTIVEVERTEXFITTERTESTALG_H + + +#include "AthenaBaseComps/AthAlgorithm.h" +#include "TrkVertexFitterInterfaces/IVertexFitter.h" +#include "GaudiKernel/ToolHandle.h" + + + +namespace Trk { + + +class AdaptiveVertexFitterTestAlg + : public AthAlgorithm +{ +public: + using AthAlgorithm::AthAlgorithm; + + + /// Standard Gaudi initialize method. + virtual StatusCode initialize() override; + + /// Execute the algorithm. + virtual StatusCode execute() override; + + +private: + StatusCode test1(); + StatusCode test2(); + StatusCode test3(); + StatusCode test4(); + StatusCode test5(); + StatusCode test6(); + + + ToolHandle<Trk::IVertexFitter> m_fitter + { this, "Tool", "Trk::AdaptiveVertexFitter", "Tool to test." }; +}; + + +} // namespace Trk + + +#endif // not TRKVERTEXFITTERS_ADAPTIVEVERTEXFITTERTESTALG_H diff --git a/Tracking/TrkVertexFitter/TrkVertexFitters/src/components/TrkVertexFitters_entries.cxx b/Tracking/TrkVertexFitter/TrkVertexFitters/src/components/TrkVertexFitters_entries.cxx index 3306f166b690e34be1205ecacefe1b450d39869c..901907d71ee188aa989b1c9e614602b24d832429 100644 --- a/Tracking/TrkVertexFitter/TrkVertexFitters/src/components/TrkVertexFitters_entries.cxx +++ b/Tracking/TrkVertexFitter/TrkVertexFitters/src/components/TrkVertexFitters_entries.cxx @@ -3,6 +3,8 @@ #include "TrkVertexFitters/AdaptiveVertexFitter.h" #include "TrkVertexFitters/SequentialVertexSmoother.h" #include "TrkVertexFitters/DummyVertexSmoother.h" +#include "../AdaptiveVertexFitterTestAlg.h" +#include "../AdaptiveMultiVertexFitterTestAlg.h" using namespace Trk; @@ -11,4 +13,6 @@ DECLARE_COMPONENT( AdaptiveMultiVertexFitter ) DECLARE_COMPONENT( SequentialVertexSmoother ) DECLARE_COMPONENT( AdaptiveVertexFitter ) DECLARE_COMPONENT( DummyVertexSmoother ) +DECLARE_COMPONENT( AdaptiveVertexFitterTestAlg ) +DECLARE_COMPONENT( AdaptiveMultiVertexFitterTestAlg )