diff --git a/CUDADataFormats/Track/BuildFile.xml b/CUDADataFormats/Track/BuildFile.xml
new file mode 100644
index 0000000000000..e3f9a0910bbd8
--- /dev/null
+++ b/CUDADataFormats/Track/BuildFile.xml
@@ -0,0 +1,9 @@
+
+
+
+
+
+
+
+
+
diff --git a/CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h b/CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h
new file mode 100644
index 0000000000000..3ee5af80353dd
--- /dev/null
+++ b/CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h
@@ -0,0 +1,9 @@
+#ifndef CUDADataFormats_Track_PixelTrackHeterogeneous_h
+#define CUDADataFormats_Track_PixelTrackHeterogeneous_h
+
+#include "CUDADataFormats/Common/interface/HeterogeneousSoA.h"
+#include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h"
+
+using PixelTrackHeterogeneous = HeterogeneousSoA;
+
+#endif // #ifndef CUDADataFormats_Track_PixelTrackHeterogeneous_h
\ No newline at end of file
diff --git a/CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h b/CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h
new file mode 100644
index 0000000000000..bd39f3c4d3bfe
--- /dev/null
+++ b/CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h
@@ -0,0 +1,73 @@
+#ifndef CUDADataFormats_Track_TrackHeterogeneousT_H
+#define CUDADataFormats_Track_TrackHeterogeneousT_H
+
+#include "CUDADataFormats/Track/interface/TrajectoryStateSoAT.h"
+#include "HeterogeneousCore/CUDAUtilities/interface/HistoContainer.h"
+
+#include "CUDADataFormats/Common/interface/HeterogeneousSoA.h"
+
+namespace pixelTrack {
+ enum class Quality : uint8_t { bad = 0, dup, loose, strict, tight, highPurity };
+}
+
+template
+class TrackSoAHeterogeneousT {
+public:
+ static constexpr int32_t stride() { return S; }
+
+ using Quality = pixelTrack::Quality;
+ using hindex_type = uint32_t;
+ using HitContainer = cms::cuda::OneToManyAssoc;
+
+ // Always check quality is at least loose!
+ // CUDA does not support enums in __lgc ...
+private:
+ eigenSoA::ScalarSoA quality_;
+
+public:
+ constexpr Quality quality(int32_t i) const { return (Quality)(quality_(i)); }
+ constexpr Quality &quality(int32_t i) { return (Quality &)(quality_(i)); }
+ constexpr Quality const *qualityData() const { return (Quality const *)(quality_.data()); }
+ constexpr Quality *qualityData() { return (Quality *)(quality_.data()); }
+
+ // this is chi2/ndof as not necessarely all hits are used in the fit
+ eigenSoA::ScalarSoA chi2;
+
+ constexpr int nHits(int i) const { return detIndices.size(i); }
+
+ // State at the Beam spot
+ // phi,tip,1/pt,cotan(theta),zip
+ TrajectoryStateSoAT stateAtBS;
+ eigenSoA::ScalarSoA eta;
+ eigenSoA::ScalarSoA pt;
+ constexpr float charge(int32_t i) const { return std::copysign(1.f, stateAtBS.state(i)(2)); }
+ constexpr float phi(int32_t i) const { return stateAtBS.state(i)(0); }
+ constexpr float tip(int32_t i) const { return stateAtBS.state(i)(1); }
+ constexpr float zip(int32_t i) const { return stateAtBS.state(i)(4); }
+
+ // state at the detector of the outermost hit
+ // representation to be decided...
+ // not yet filled on GPU
+ // TrajectoryStateSoA stateAtOuterDet;
+
+ HitContainer hitIndices;
+ HitContainer detIndices;
+};
+
+namespace pixelTrack {
+
+#ifdef GPU_SMALL_EVENTS
+ // kept for testing and debugging
+ constexpr uint32_t maxNumber() { return 2 * 1024; }
+#else
+ // tested on MC events with 55-75 pileup events
+ constexpr uint32_t maxNumber() { return 32 * 1024; }
+#endif
+
+ using TrackSoA = TrackSoAHeterogeneousT;
+ using TrajectoryState = TrajectoryStateSoAT;
+ using HitContainer = TrackSoA::HitContainer;
+
+} // namespace pixelTrack
+
+#endif // CUDADataFormats_Track_TrackHeterogeneousT_H
diff --git a/CUDADataFormats/Track/interface/TrajectoryStateSoAT.h b/CUDADataFormats/Track/interface/TrajectoryStateSoAT.h
new file mode 100644
index 0000000000000..64fcd573a6991
--- /dev/null
+++ b/CUDADataFormats/Track/interface/TrajectoryStateSoAT.h
@@ -0,0 +1,59 @@
+#ifndef CUDADataFormats_Track_TrajectoryStateSOAT_H
+#define CUDADataFormats_Track_TrajectoryStateSOAT_H
+
+#include
+#include "HeterogeneousCore/CUDAUtilities/interface/eigenSoA.h"
+
+template
+struct TrajectoryStateSoAT {
+ using Vector5f = Eigen::Matrix;
+ using Vector15f = Eigen::Matrix;
+
+ using Vector5d = Eigen::Matrix;
+ using Matrix5d = Eigen::Matrix;
+
+ static constexpr int32_t stride() { return S; }
+
+ eigenSoA::MatrixSoA state;
+ eigenSoA::MatrixSoA covariance;
+
+ template
+ __host__ __device__ inline void copyFromCircle(
+ V3 const& cp, M3 const& ccov, V2 const& lp, M2 const& lcov, float b, int32_t i) {
+ state(i) << cp.template cast(), lp.template cast();
+ state(i)(2) *= b;
+ auto cov = covariance(i);
+ cov(0) = ccov(0, 0);
+ cov(1) = ccov(0, 1);
+ cov(2) = b * float(ccov(0, 2));
+ cov(4) = cov(3) = 0;
+ cov(5) = ccov(1, 1);
+ cov(6) = b * float(ccov(1, 2));
+ cov(8) = cov(7) = 0;
+ cov(9) = b * b * float(ccov(2, 2));
+ cov(11) = cov(10) = 0;
+ cov(12) = lcov(0, 0);
+ cov(13) = lcov(0, 1);
+ cov(14) = lcov(1, 1);
+ }
+
+ template
+ __host__ __device__ inline void copyFromDense(V5 const& v, M5 const& cov, int32_t i) {
+ state(i) = v.template cast();
+ for (int j = 0, ind = 0; j < 5; ++j)
+ for (auto k = j; k < 5; ++k)
+ covariance(i)(ind++) = cov(j, k);
+ }
+
+ template
+ __host__ __device__ inline void copyToDense(V5& v, M5& cov, int32_t i) const {
+ v = state(i).template cast();
+ for (int j = 0, ind = 0; j < 5; ++j) {
+ cov(j, j) = covariance(i)(ind++);
+ for (auto k = j + 1; k < 5; ++k)
+ cov(k, j) = cov(j, k) = covariance(i)(ind++);
+ }
+ }
+};
+
+#endif // CUDADataFormats_Track_TrajectoryStateSOAT_H
diff --git a/CUDADataFormats/Track/src/classes.h b/CUDADataFormats/Track/src/classes.h
new file mode 100644
index 0000000000000..97c116f6c88d3
--- /dev/null
+++ b/CUDADataFormats/Track/src/classes.h
@@ -0,0 +1,9 @@
+#ifndef CUDADataFormats_Track_src_classes_h
+#define CUDADataFormats_Track_src_classes_h
+
+#include "CUDADataFormats/Common/interface/Product.h"
+#include "CUDADataFormats/Common/interface/HostProduct.h"
+#include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h"
+#include "DataFormats/Common/interface/Wrapper.h"
+
+#endif // CUDADataFormats_Track_src_classes_h
diff --git a/CUDADataFormats/Track/src/classes_def.xml b/CUDADataFormats/Track/src/classes_def.xml
new file mode 100644
index 0000000000000..9c80ae91baf29
--- /dev/null
+++ b/CUDADataFormats/Track/src/classes_def.xml
@@ -0,0 +1,6 @@
+
+
+
+
+
+
diff --git a/CUDADataFormats/Track/test/BuildFile.xml b/CUDADataFormats/Track/test/BuildFile.xml
new file mode 100644
index 0000000000000..598b345d4709d
--- /dev/null
+++ b/CUDADataFormats/Track/test/BuildFile.xml
@@ -0,0 +1,13 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/CUDADataFormats/Track/test/TrajectoryStateSOA_t.cpp b/CUDADataFormats/Track/test/TrajectoryStateSOA_t.cpp
new file mode 100644
index 0000000000000..d6ff539a642b0
--- /dev/null
+++ b/CUDADataFormats/Track/test/TrajectoryStateSOA_t.cpp
@@ -0,0 +1 @@
+#include "TrajectoryStateSOA_t.h"
diff --git a/CUDADataFormats/Track/test/TrajectoryStateSOA_t.cu b/CUDADataFormats/Track/test/TrajectoryStateSOA_t.cu
new file mode 100644
index 0000000000000..d6ff539a642b0
--- /dev/null
+++ b/CUDADataFormats/Track/test/TrajectoryStateSOA_t.cu
@@ -0,0 +1 @@
+#include "TrajectoryStateSOA_t.h"
diff --git a/CUDADataFormats/Track/test/TrajectoryStateSOA_t.h b/CUDADataFormats/Track/test/TrajectoryStateSOA_t.h
new file mode 100644
index 0000000000000..97b88873c2613
--- /dev/null
+++ b/CUDADataFormats/Track/test/TrajectoryStateSOA_t.h
@@ -0,0 +1,75 @@
+#include "CUDADataFormats/Track/interface/TrajectoryStateSoAT.h"
+
+using Vector5d = Eigen::Matrix;
+using Matrix5d = Eigen::Matrix;
+
+__host__ __device__ Matrix5d loadCov(Vector5d const& e) {
+ Matrix5d cov;
+ for (int i = 0; i < 5; ++i)
+ cov(i, i) = e(i) * e(i);
+ for (int i = 0; i < 5; ++i) {
+ for (int j = 0; j < i; ++j) {
+ double v = 0.3 * std::sqrt(cov(i, i) * cov(j, j)); // this makes the matrix pos defined
+ cov(i, j) = (i + j) % 2 ? -0.4 * v : 0.1 * v;
+ cov(j, i) = cov(i, j);
+ }
+ }
+ return cov;
+}
+
+using TS = TrajectoryStateSoAT<128>;
+
+__global__ void testTSSoA(TS* pts, int n) {
+ assert(n <= 128);
+
+ Vector5d par0;
+ par0 << 0.2, 0.1, 3.5, 0.8, 0.1;
+ Vector5d e0;
+ e0 << 0.01, 0.01, 0.035, -0.03, -0.01;
+ auto cov0 = loadCov(e0);
+
+ TS& ts = *pts;
+
+ int first = threadIdx.x + blockIdx.x * blockDim.x;
+
+ for (int i = first; i < n; i += blockDim.x * gridDim.x) {
+ ts.copyFromDense(par0, cov0, i);
+ Vector5d par1;
+ Matrix5d cov1;
+ ts.copyToDense(par1, cov1, i);
+ Vector5d delV = par1 - par0;
+ Matrix5d delM = cov1 - cov0;
+ for (int j = 0; j < 5; ++j) {
+ assert(std::abs(delV(j)) < 1.e-5);
+ for (auto k = j; k < 5; ++k) {
+ assert(cov0(k, j) == cov0(j, k));
+ assert(cov1(k, j) == cov1(j, k));
+ assert(std::abs(delM(k, j)) < 1.e-5);
+ }
+ }
+ }
+}
+
+#ifdef __CUDACC__
+#include "HeterogeneousCore/CUDAUtilities/interface/requireDevices.h"
+#include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h"
+#endif
+
+int main() {
+#ifdef __CUDACC__
+ cms::cudatest::requireDevices();
+#endif
+
+ TS ts;
+
+#ifdef __CUDACC__
+ TS* ts_d;
+ cudaCheck(cudaMalloc(&ts_d, sizeof(TS)));
+ testTSSoA<<<1, 64>>>(ts_d, 128);
+ cudaCheck(cudaGetLastError());
+ cudaCheck(cudaMemcpy(&ts, ts_d, sizeof(TS), cudaMemcpyDefault));
+ cudaCheck(cudaDeviceSynchronize());
+#else
+ testTSSoA(&ts, 128);
+#endif
+}
diff --git a/DQM/TrackingMonitorClient/python/pixelTrackingEffFromHitPattern_cff.py b/DQM/TrackingMonitorClient/python/pixelTrackingEffFromHitPattern_cff.py
index 15ceaf93ed20a..cff85e56d94f7 100644
--- a/DQM/TrackingMonitorClient/python/pixelTrackingEffFromHitPattern_cff.py
+++ b/DQM/TrackingMonitorClient/python/pixelTrackingEffFromHitPattern_cff.py
@@ -21,7 +21,10 @@ def _layers(suffix, quant, histoPostfix):
]
pixelTrackingEffFromHitPattern = DQMEDHarvester("DQMGenericClient",
- subDirs = cms.untracked.vstring("Tracking/PixelTrackParameters/HitEffFromHitPattern*"),
+ subDirs = cms.untracked.vstring("Tracking/PixelTrackParameters/pixelTracks/HitEffFromHitPattern*",
+ "Tracking/PixelTrackParameters/dzPV0p1/HitEffFromHitPattern*",
+ "Tracking/PixelTrackParameters/pt_0to1/HitEffFromHitPattern*",
+ "Tracking/PixelTrackParameters/pt_1/HitEffFromHitPattern*"),
efficiency = cms.vstring(
_layers("PU", "GoodNumVertices", "") +
_layers("BX", "BX", "VsBX") +
diff --git a/DQM/TrackingMonitorSource/python/pixelTracksMonitoring_cff.py b/DQM/TrackingMonitorSource/python/pixelTracksMonitoring_cff.py
index a075f671f05ce..d5deba78b46c8 100644
--- a/DQM/TrackingMonitorSource/python/pixelTracksMonitoring_cff.py
+++ b/DQM/TrackingMonitorSource/python/pixelTracksMonitoring_cff.py
@@ -1,23 +1,77 @@
import FWCore.ParameterSet.Config as cms
import DQM.TrackingMonitor.TrackerCollisionTrackingMonitor_cfi
-pixelTracksMonitoring = DQM.TrackingMonitor.TrackerCollisionTrackingMonitor_cfi.TrackerCollisionTrackMon.clone()
-pixelTracksMonitoring.FolderName = 'Tracking/PixelTrackParameters'
-pixelTracksMonitoring.TrackProducer = 'pixelTracks'
-pixelTracksMonitoring.allTrackProducer = 'pixelTracks'
-pixelTracksMonitoring.beamSpot = 'offlineBeamSpot'
-pixelTracksMonitoring.primaryVertex = 'pixelVertices'
-pixelTracksMonitoring.pvNDOF = 1
-pixelTracksMonitoring.doAllPlots = True
-pixelTracksMonitoring.doLumiAnalysis = True
-pixelTracksMonitoring.doProfilesVsLS = True
-pixelTracksMonitoring.doDCAPlots = True
-pixelTracksMonitoring.doProfilesVsLS = True
-pixelTracksMonitoring.doPlotsVsGoodPVtx = True
-pixelTracksMonitoring.doEffFromHitPatternVsPU = False
-pixelTracksMonitoring.doEffFromHitPatternVsBX = False
-pixelTracksMonitoring.doEffFromHitPatternVsLUMI = False
-pixelTracksMonitoring.doPlotsVsGoodPVtx = True
-pixelTracksMonitoring.doPlotsVsLUMI = True
-pixelTracksMonitoring.doPlotsVsBX = True
+pixelTracksMonitor = DQM.TrackingMonitor.TrackerCollisionTrackingMonitor_cfi.TrackerCollisionTrackMon.clone()
+pixelTracksMonitor.FolderName = 'Tracking/PixelTrackParameters/pixelTracks'
+pixelTracksMonitor.TrackProducer = 'pixelTracks'
+pixelTracksMonitor.allTrackProducer = 'pixelTracks'
+pixelTracksMonitor.beamSpot = 'offlineBeamSpot'
+pixelTracksMonitor.primaryVertex = 'pixelVertices'
+pixelTracksMonitor.pvNDOF = 1
+pixelTracksMonitor.doAllPlots = True
+pixelTracksMonitor.doLumiAnalysis = True
+pixelTracksMonitor.doProfilesVsLS = True
+pixelTracksMonitor.doDCAPlots = True
+pixelTracksMonitor.doProfilesVsLS = True
+pixelTracksMonitor.doPlotsVsGoodPVtx = True
+pixelTracksMonitor.doEffFromHitPatternVsPU = False
+pixelTracksMonitor.doEffFromHitPatternVsBX = False
+pixelTracksMonitor.doEffFromHitPatternVsLUMI = False
+pixelTracksMonitor.doPlotsVsGoodPVtx = True
+pixelTracksMonitor.doPlotsVsLUMI = True
+pixelTracksMonitor.doPlotsVsBX = True
+_trackSelector = cms.EDFilter('TrackSelector',
+ src = cms.InputTag('pixelTracks'),
+ cut = cms.string("")
+)
+
+pixelTracksPt0to1 = _trackSelector.clone(cut = "pt >= 0 & pt < 1 ")
+pixelTracksPt1 = _trackSelector.clone(cut = "pt >= 1 ")
+from DQM.TrackingMonitorSource.TrackCollections2monitor_cff import highPurityPV0p1 as _highPurityPV0p1
+pixelTracksPV0p1 = _highPurityPV0p1.clone(
+ src = "pixelTracks",
+ quality = "",
+ vertexTag = "goodPixelVertices"
+)
+
+pixelTracksMonitorPt0to1 = pixelTracksMonitor.clone(
+ TrackProducer = "pixelTracksPt0to1",
+ FolderName = "Tracking/PixelTrackParameters/pt_0to1"
+)
+pixelTracksMonitorPt1 = pixelTracksMonitor.clone(
+ TrackProducer = "pixelTracksPt1",
+ FolderName = "Tracking/PixelTrackParameters/pt_1"
+)
+pixelTracksMonitorPV0p1 = pixelTracksMonitor.clone(
+ TrackProducer = "pixelTracksPV0p1",
+ FolderName = "Tracking/PixelTrackParameters/dzPV0p1"
+)
+
+
+from CommonTools.ParticleFlow.goodOfflinePrimaryVertices_cfi import goodOfflinePrimaryVertices as _goodOfflinePrimaryVertices
+goodPixelVertices = _goodOfflinePrimaryVertices.clone(
+ src = "pixelVertices",
+)
+
+from DQM.TrackingMonitor.primaryVertexResolution_cfi import primaryVertexResolution as _primaryVertexResolution
+pixelVertexResolution = _primaryVertexResolution.clone(
+ vertexSrc = "goodPixelVertices",
+ rootFolder = "OfflinePixelPV/Resolution",
+)
+
+pixelTracksMonitoringTask = cms.Task(
+ goodPixelVertices,
+ pixelTracksPt0to1,
+ pixelTracksPt1,
+ pixelTracksPV0p1,
+)
+
+pixelTracksMonitoring = cms.Sequence(
+ pixelTracksMonitor +
+ pixelTracksMonitorPt0to1 +
+ pixelTracksMonitorPt1 +
+ pixelTracksMonitorPV0p1 +
+ pixelVertexResolution,
+ pixelTracksMonitoringTask
+)
diff --git a/RecoPixelVertexing/Configuration/python/customizePixelTracksForProfiling.py b/RecoPixelVertexing/Configuration/python/customizePixelTracksForProfiling.py
deleted file mode 100644
index 4713b64e5e48a..0000000000000
--- a/RecoPixelVertexing/Configuration/python/customizePixelTracksForProfiling.py
+++ /dev/null
@@ -1,15 +0,0 @@
-import FWCore.ParameterSet.Config as cms
-
-def customizePixelTracksForProfiling(process):
- process.out = cms.OutputModule("AsciiOutputModule",
- outputCommands = cms.untracked.vstring(
- "keep *_pixelTracks_*_*",
- ),
- verbosity = cms.untracked.uint32(0),
- )
-
- process.outPath = cms.EndPath(process.out)
-
- process.schedule = cms.Schedule(process.raw2digi_step, process.reconstruction_step, process.outPath)
-
- return process
diff --git a/RecoPixelVertexing/Configuration/python/customizePixelTracksSoAonCPU.py b/RecoPixelVertexing/Configuration/python/customizePixelTracksSoAonCPU.py
new file mode 100644
index 0000000000000..24cc16e02b463
--- /dev/null
+++ b/RecoPixelVertexing/Configuration/python/customizePixelTracksSoAonCPU.py
@@ -0,0 +1,53 @@
+import FWCore.ParameterSet.Config as cms
+
+def customizePixelTracksSoAonCPU(process):
+
+ process.CUDAService = cms.Service('CUDAService',
+ enabled = cms.untracked.bool(False)
+ )
+
+ # ensure the same results when running on GPU (which supports only the 'HLT' payload) and CPU
+ process.siPixelClustersPreSplitting.cpu.payloadType = cms.string('HLT')
+
+ from RecoLocalTracker.SiPixelRecHits.siPixelRecHitSoAFromLegacy_cfi import siPixelRecHitSoAFromLegacy
+ process.siPixelRecHitsPreSplitting = siPixelRecHitSoAFromLegacy.clone(
+ convertToLegacy = True
+ )
+
+ from RecoPixelVertexing.PixelTriplets.caHitNtupletCUDA_cfi import caHitNtupletCUDA
+ process.pixelTrackSoA = caHitNtupletCUDA.clone(
+ onGPU = False,
+ pixelRecHitSrc = 'siPixelRecHitsPreSplitting'
+ )
+
+ from RecoPixelVertexing.PixelTrackFitting.pixelTrackProducerFromSoA_cfi import pixelTrackProducerFromSoA
+ process.pixelTracks = pixelTrackProducerFromSoA.clone(
+ pixelRecHitLegacySrc = 'siPixelRecHitsPreSplitting'
+ )
+
+ process.reconstruction_step += process.siPixelRecHitsPreSplitting + process.pixelTrackSoA
+
+ return process
+
+
+def customizePixelTracksForTriplets(process):
+
+ from HLTrigger.Configuration.common import producers_by_type
+ for producer in producers_by_type(process, 'CAHitNtupletCUDA'):
+ producer.includeJumpingForwardDoublets = True
+ producer.minHitsPerNtuplet = 3
+
+ return process
+
+
+def customizePixelTracksSoAonCPUForProfiling(process):
+
+ process = customizePixelTracksSoAonCPU(process)
+
+ process.siPixelRecHitSoAFromLegacy.convertToLegacy = False
+
+ process.TkSoA = cms.Path(process.offlineBeamSpot + process.siPixelDigis + process.siPixelClustersPreSplitting + process.siPixelRecHitSoAFromLegacy + process.pixelTrackSoA)
+
+ process.schedule = cms.Schedule(process.TkSoA)
+
+ return process
diff --git a/RecoPixelVertexing/PixelTrackFitting/BuildFile.xml b/RecoPixelVertexing/PixelTrackFitting/BuildFile.xml
index e6fc938dc25a7..a589aad036996 100644
--- a/RecoPixelVertexing/PixelTrackFitting/BuildFile.xml
+++ b/RecoPixelVertexing/PixelTrackFitting/BuildFile.xml
@@ -1,3 +1,5 @@
+
+
diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h b/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h
new file mode 100644
index 0000000000000..86fe6a278777c
--- /dev/null
+++ b/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h
@@ -0,0 +1,606 @@
+#ifndef RecoPixelVertexing_PixelTrackFitting_interface_BrokenLine_h
+#define RecoPixelVertexing_PixelTrackFitting_interface_BrokenLine_h
+
+#include
+
+#include "RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h"
+
+namespace brokenline {
+
+ //!< Karimäki's parameters: (phi, d, k=1/R)
+ /*!< covariance matrix: \n
+ |cov(phi,phi)|cov( d ,phi)|cov( k ,phi)| \n
+ |cov(phi, d )|cov( d , d )|cov( k , d )| \n
+ |cov(phi, k )|cov( d , k )|cov( k , k )| \n
+ as defined in Karimäki V., 1990, Effective circle fitting for particle trajectories,
+ Nucl. Instr. and Meth. A305 (1991) 187.
+ */
+ using karimaki_circle_fit = riemannFit::CircleFit;
+
+ /*!
+ \brief data needed for the Broken Line fit procedure.
+ */
+ template
+ struct PreparedBrokenLineData {
+ int qCharge; //!< particle charge
+ riemannFit::Matrix2xNd radii; //!< xy data in the system in which the pre-fitted center is the origin
+ riemannFit::VectorNd sTransverse; //!< total distance traveled in the transverse plane
+ // starting from the pre-fitted closest approach
+ riemannFit::VectorNd sTotal; //!< total distance traveled (three-dimensional)
+ riemannFit::VectorNd zInSZplane; //!< orthogonal coordinate to the pre-fitted line in the sz plane
+ riemannFit::VectorNd varBeta; //!< kink angles in the SZ plane
+ };
+
+ /*!
+ \brief Computes the Coulomb multiple scattering variance of the planar angle.
+
+ \param length length of the track in the material.
+ \param bField magnetic field in Gev/cm/c.
+ \param radius radius of curvature (needed to evaluate p).
+ \param layer denotes which of the four layers of the detector is the endpoint of the
+ * multiple scattered track. For example, if Layer=3, then the particle has
+ * just gone through the material between the second and the third layer.
+
+ \todo add another Layer variable to identify also the start point of the track,
+ * so if there are missing hits or multiple hits, the part of the detector that
+ * the particle has traversed can be exactly identified.
+
+ \warning the formula used here assumes beta=1, and so neglects the dependence
+ * of theta_0 on the mass of the particle at fixed momentum.
+
+ \return the variance of the planar angle ((theta_0)^2 /3).
+ */
+ __host__ __device__ inline double multScatt(
+ const double& length, const double bField, const double radius, int layer, double slope) {
+ // limit R to 20GeV...
+ auto pt2 = std::min(20., bField * radius);
+ pt2 *= pt2;
+ constexpr double inv_X0 = 0.06 / 16.; //!< inverse of radiation length of the material in cm
+ //if(Layer==1) XXI_0=0.06/16.;
+ // else XXI_0=0.06/16.;
+ //XX_0*=1;
+
+ //! number between 1/3 (uniform material) and 1 (thin scatterer) to be manually tuned
+ constexpr double geometry_factor = 0.7;
+ constexpr double fact = geometry_factor * riemannFit::sqr(13.6 / 1000.);
+ return fact / (pt2 * (1. + riemannFit::sqr(slope))) * (std::abs(length) * inv_X0) *
+ riemannFit::sqr(1. + 0.038 * log(std::abs(length) * inv_X0));
+ }
+
+ /*!
+ \brief Computes the 2D rotation matrix that transforms the line y=slope*x into the line y=0.
+
+ \param slope tangent of the angle of rotation.
+
+ \return 2D rotation matrix.
+ */
+ __host__ __device__ inline riemannFit::Matrix2d rotationMatrix(double slope) {
+ riemannFit::Matrix2d rot;
+ rot(0, 0) = 1. / sqrt(1. + riemannFit::sqr(slope));
+ rot(0, 1) = slope * rot(0, 0);
+ rot(1, 0) = -rot(0, 1);
+ rot(1, 1) = rot(0, 0);
+ return rot;
+ }
+
+ /*!
+ \brief Changes the Karimäki parameters (and consequently their covariance matrix) under a
+ * translation of the coordinate system, such that the old origin has coordinates (x0,y0)
+ * in the new coordinate system. The formulas are taken from Karimäki V., 1990, Effective
+ * circle fitting for particle trajectories, Nucl. Instr. and Meth. A305 (1991) 187.
+
+ \param circle circle fit in the old coordinate system. circle.par(0) is phi, circle.par(1) is d and circle.par(2) is rho.
+ \param x0 x coordinate of the translation vector.
+ \param y0 y coordinate of the translation vector.
+ \param jacobian passed by reference in order to save stack.
+ */
+ __host__ __device__ inline void translateKarimaki(karimaki_circle_fit& circle,
+ double x0,
+ double y0,
+ riemannFit::Matrix3d& jacobian) {
+ // Avoid multiple access to the circle.par vector.
+ using scalar = std::remove_reference::type;
+ scalar phi = circle.par(0);
+ scalar dee = circle.par(1);
+ scalar rho = circle.par(2);
+
+ // Avoid repeated trig. computations
+ scalar sinPhi = sin(phi);
+ scalar cosPhi = cos(phi);
+
+ // Intermediate computations for the circle parameters
+ scalar deltaPara = x0 * cosPhi + y0 * sinPhi;
+ scalar deltaOrth = x0 * sinPhi - y0 * cosPhi + dee;
+ scalar tempSmallU = 1 + rho * dee;
+ scalar tempC = -rho * y0 + tempSmallU * cosPhi;
+ scalar tempB = rho * x0 + tempSmallU * sinPhi;
+ scalar tempA = 2. * deltaOrth + rho * (riemannFit::sqr(deltaOrth) + riemannFit::sqr(deltaPara));
+ scalar tempU = sqrt(1. + rho * tempA);
+
+ // Intermediate computations for the error matrix transform
+ scalar xi = 1. / (riemannFit::sqr(tempB) + riemannFit::sqr(tempC));
+ scalar tempV = 1. + rho * deltaOrth;
+ scalar lambda = (0.5 * tempA) / (riemannFit::sqr(1. + tempU) * tempU);
+ scalar mu = 1. / (tempU * (1. + tempU)) + rho * lambda;
+ scalar zeta = riemannFit::sqr(deltaOrth) + riemannFit::sqr(deltaPara);
+ jacobian << xi * tempSmallU * tempV, -xi * riemannFit::sqr(rho) * deltaOrth, xi * deltaPara,
+ 2. * mu * tempSmallU * deltaPara, 2. * mu * tempV, mu * zeta - lambda * tempA, 0, 0, 1.;
+
+ // translated circle parameters
+ // phi
+ circle.par(0) = atan2(tempB, tempC);
+ // d
+ circle.par(1) = tempA / (1 + tempU);
+ // rho after translation. It is invariant, so noop
+ // circle.par(2)= rho;
+
+ // translated error matrix
+ circle.cov = jacobian * circle.cov * jacobian.transpose();
+ }
+
+ /*!
+ \brief Computes the data needed for the Broken Line fit procedure that are mainly common for the circle and the line fit.
+
+ \param hits hits coordinates.
+ \param fast_fit pre-fit result in the form (X0,Y0,R,tan(theta)).
+ \param bField magnetic field in Gev/cm/c.
+ \param results PreparedBrokenLineData to be filled (see description of PreparedBrokenLineData).
+ */
+ template
+ __host__ __device__ inline void prepareBrokenLineData(const M3xN& hits,
+ const V4& fast_fit,
+ const double bField,
+ PreparedBrokenLineData& results) {
+ riemannFit::Vector2d dVec;
+ riemannFit::Vector2d eVec;
+
+ dVec = hits.block(0, 1, 2, 1) - hits.block(0, 0, 2, 1);
+ eVec = hits.block(0, n - 1, 2, 1) - hits.block(0, n - 2, 2, 1);
+ results.qCharge = riemannFit::cross2D(dVec, eVec) > 0 ? -1 : 1;
+
+ const double slope = -results.qCharge / fast_fit(3);
+
+ riemannFit::Matrix2d rotMat = rotationMatrix(slope);
+
+ // calculate radii and s
+ results.radii = hits.block(0, 0, 2, n) - fast_fit.head(2) * riemannFit::MatrixXd::Constant(1, n, 1);
+ eVec = -fast_fit(2) * fast_fit.head(2) / fast_fit.head(2).norm();
+ for (u_int i = 0; i < n; i++) {
+ dVec = results.radii.block(0, i, 2, 1);
+ results.sTransverse(i) = results.qCharge * fast_fit(2) *
+ atan2(riemannFit::cross2D(dVec, eVec), dVec.dot(eVec)); // calculates the arc length
+ }
+ riemannFit::VectorNd zVec = hits.block(2, 0, 1, n).transpose();
+
+ //calculate sTotal and zVec
+ riemannFit::Matrix2xNd pointsSZ = riemannFit::Matrix2xNd::Zero();
+ for (u_int i = 0; i < n; i++) {
+ pointsSZ(0, i) = results.sTransverse(i);
+ pointsSZ(1, i) = zVec(i);
+ pointsSZ.block(0, i, 2, 1) = rotMat * pointsSZ.block(0, i, 2, 1);
+ }
+ results.sTotal = pointsSZ.block(0, 0, 1, n).transpose();
+ results.zInSZplane = pointsSZ.block(1, 0, 1, n).transpose();
+
+ //calculate varBeta
+ results.varBeta(0) = results.varBeta(n - 1) = 0;
+ for (u_int i = 1; i < n - 1; i++) {
+ results.varBeta(i) = multScatt(results.sTotal(i + 1) - results.sTotal(i), bField, fast_fit(2), i + 2, slope) +
+ multScatt(results.sTotal(i) - results.sTotal(i - 1), bField, fast_fit(2), i + 1, slope);
+ }
+ }
+
+ /*!
+ \brief Computes the n-by-n band matrix obtained minimizing the Broken Line's cost function w.r.t u.
+ * This is the whole matrix in the case of the line fit and the main n-by-n block in the case
+ * of the circle fit.
+
+ \param weights weights of the first part of the cost function, the one with the measurements
+ * and not the angles (\sum_{i=1}^n w*(y_i-u_i)^2).
+ \param sTotal total distance traveled by the particle from the pre-fitted closest approach.
+ \param varBeta kink angles' variance.
+
+ \return the n-by-n matrix of the linear system
+ */
+ template
+ __host__ __device__ inline riemannFit::MatrixNd matrixC_u(const riemannFit::VectorNd& weights,
+ const riemannFit::VectorNd& sTotal,
+ const riemannFit::VectorNd& varBeta) {
+ riemannFit::MatrixNd c_uMat = riemannFit::MatrixNd::Zero();
+ for (u_int i = 0; i < n; i++) {
+ c_uMat(i, i) = weights(i);
+ if (i > 1)
+ c_uMat(i, i) += 1. / (varBeta(i - 1) * riemannFit::sqr(sTotal(i) - sTotal(i - 1)));
+ if (i > 0 && i < n - 1)
+ c_uMat(i, i) +=
+ (1. / varBeta(i)) * riemannFit::sqr((sTotal(i + 1) - sTotal(i - 1)) /
+ ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
+ if (i < n - 2)
+ c_uMat(i, i) += 1. / (varBeta(i + 1) * riemannFit::sqr(sTotal(i + 1) - sTotal(i)));
+
+ if (i > 0 && i < n - 1)
+ c_uMat(i, i + 1) =
+ 1. / (varBeta(i) * (sTotal(i + 1) - sTotal(i))) *
+ (-(sTotal(i + 1) - sTotal(i - 1)) / ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
+ if (i < n - 2)
+ c_uMat(i, i + 1) +=
+ 1. / (varBeta(i + 1) * (sTotal(i + 1) - sTotal(i))) *
+ (-(sTotal(i + 2) - sTotal(i)) / ((sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i))));
+
+ if (i < n - 2)
+ c_uMat(i, i + 2) = 1. / (varBeta(i + 1) * (sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i)));
+
+ c_uMat(i, i) *= 0.5;
+ }
+ return c_uMat + c_uMat.transpose();
+ }
+
+ /*!
+ \brief A very fast helix fit.
+
+ \param hits the measured hits.
+
+ \return (X0,Y0,R,tan(theta)).
+
+ \warning sign of theta is (intentionally, for now) mistaken for negative charges.
+ */
+
+ template
+ __host__ __device__ inline void fastFit(const M3xN& hits, V4& result) {
+ constexpr uint32_t n = M3xN::ColsAtCompileTime;
+
+ const riemannFit::Vector2d a = hits.block(0, n / 2, 2, 1) - hits.block(0, 0, 2, 1);
+ const riemannFit::Vector2d b = hits.block(0, n - 1, 2, 1) - hits.block(0, n / 2, 2, 1);
+ const riemannFit::Vector2d c = hits.block(0, 0, 2, 1) - hits.block(0, n - 1, 2, 1);
+
+ auto tmp = 0.5 / riemannFit::cross2D(c, a);
+ result(0) = hits(0, 0) - (a(1) * c.squaredNorm() + c(1) * a.squaredNorm()) * tmp;
+ result(1) = hits(1, 0) + (a(0) * c.squaredNorm() + c(0) * a.squaredNorm()) * tmp;
+ // check Wikipedia for these formulas
+
+ result(2) = sqrt(a.squaredNorm() * b.squaredNorm() * c.squaredNorm()) / (2. * std::abs(riemannFit::cross2D(b, a)));
+ // Using Math Olympiad's formula R=abc/(4A)
+
+ const riemannFit::Vector2d d = hits.block(0, 0, 2, 1) - result.head(2);
+ const riemannFit::Vector2d e = hits.block(0, n - 1, 2, 1) - result.head(2);
+
+ result(3) = result(2) * atan2(riemannFit::cross2D(d, e), d.dot(e)) / (hits(2, n - 1) - hits(2, 0));
+ // ds/dz slope between last and first point
+ }
+
+ /*!
+ \brief Performs the Broken Line fit in the curved track case (that is, the fit
+ * parameters are the interceptions u and the curvature correction \Delta\kappa).
+
+ \param hits hits coordinates.
+ \param hits_cov hits covariance matrix.
+ \param fast_fit pre-fit result in the form (X0,Y0,R,tan(theta)).
+ \param bField magnetic field in Gev/cm/c.
+ \param data PreparedBrokenLineData.
+ \param circle_results struct to be filled with the results in this form:
+ -par parameter of the line in this form: (phi, d, k); \n
+ -cov covariance matrix of the fitted parameter; \n
+ -chi2 value of the cost function in the minimum.
+
+ \details The function implements the steps 2 and 3 of the Broken Line fit
+ * with the curvature correction.\n
+ * The step 2 is the least square fit, done by imposing the minimum constraint on
+ * the cost function and solving the consequent linear system. It determines the
+ * fitted parameters u and \Delta\kappa and their covariance matrix.
+ * The step 3 is the correction of the fast pre-fitted parameters for the innermost
+ * part of the track. It is first done in a comfortable coordinate system (the one
+ * in which the first hit is the origin) and then the parameters and their
+ * covariance matrix are transformed to the original coordinate system.
+ */
+ template
+ __host__ __device__ inline void circleFit(const M3xN& hits,
+ const M6xN& hits_ge,
+ const V4& fast_fit,
+ const double bField,
+ PreparedBrokenLineData& data,
+ karimaki_circle_fit& circle_results) {
+ circle_results.qCharge = data.qCharge;
+ auto& radii = data.radii;
+ const auto& sTransverse = data.sTransverse;
+ const auto& sTotal = data.sTotal;
+ auto& zInSZplane = data.zInSZplane;
+ auto& varBeta = data.varBeta;
+ const double slope = -circle_results.qCharge / fast_fit(3);
+ varBeta *= 1. + riemannFit::sqr(slope); // the kink angles are projected!
+
+ for (u_int i = 0; i < n; i++) {
+ zInSZplane(i) = radii.block(0, i, 2, 1).norm() - fast_fit(2);
+ }
+
+ riemannFit::Matrix2d vMat; // covariance matrix
+ riemannFit::VectorNd weightsVec; // weights
+ riemannFit::Matrix2d rotMat; // rotation matrix point by point
+ for (u_int i = 0; i < n; i++) {
+ vMat(0, 0) = hits_ge.col(i)[0]; // x errors
+ vMat(0, 1) = vMat(1, 0) = hits_ge.col(i)[1]; // cov_xy
+ vMat(1, 1) = hits_ge.col(i)[2]; // y errors
+ rotMat = rotationMatrix(-radii(0, i) / radii(1, i));
+ weightsVec(i) =
+ 1. / ((rotMat * vMat * rotMat.transpose())(1, 1)); // compute the orthogonal weight point by point
+ }
+
+ riemannFit::VectorNplusONEd r_uVec;
+ r_uVec(n) = 0;
+ for (u_int i = 0; i < n; i++) {
+ r_uVec(i) = weightsVec(i) * zInSZplane(i);
+ }
+
+ riemannFit::MatrixNplusONEd c_uMat;
+ c_uMat.block(0, 0, n, n) = matrixC_u(weightsVec, sTransverse, varBeta);
+ c_uMat(n, n) = 0;
+ //add the border to the c_uMat matrix
+ for (u_int i = 0; i < n; i++) {
+ c_uMat(i, n) = 0;
+ if (i > 0 && i < n - 1) {
+ c_uMat(i, n) +=
+ -(sTransverse(i + 1) - sTransverse(i - 1)) * (sTransverse(i + 1) - sTransverse(i - 1)) /
+ (2. * varBeta(i) * (sTransverse(i + 1) - sTransverse(i)) * (sTransverse(i) - sTransverse(i - 1)));
+ }
+ if (i > 1) {
+ c_uMat(i, n) +=
+ (sTransverse(i) - sTransverse(i - 2)) / (2. * varBeta(i - 1) * (sTransverse(i) - sTransverse(i - 1)));
+ }
+ if (i < n - 2) {
+ c_uMat(i, n) +=
+ (sTransverse(i + 2) - sTransverse(i)) / (2. * varBeta(i + 1) * (sTransverse(i + 1) - sTransverse(i)));
+ }
+ c_uMat(n, i) = c_uMat(i, n);
+ if (i > 0 && i < n - 1)
+ c_uMat(n, n) += riemannFit::sqr(sTransverse(i + 1) - sTransverse(i - 1)) / (4. * varBeta(i));
+ }
+
+#ifdef CPP_DUMP
+ std::cout << "CU5\n" << c_uMat << std::endl;
+#endif
+ riemannFit::MatrixNplusONEd iMat;
+ math::cholesky::invert(c_uMat, iMat);
+#ifdef CPP_DUMP
+ std::cout << "I5\n" << iMat << std::endl;
+#endif
+
+ riemannFit::VectorNplusONEd uVec = iMat * r_uVec; // obtain the fitted parameters by solving the linear system
+
+ // compute (phi, d_ca, k) in the system in which the midpoint of the first two corrected hits is the origin...
+
+ radii.block(0, 0, 2, 1) /= radii.block(0, 0, 2, 1).norm();
+ radii.block(0, 1, 2, 1) /= radii.block(0, 1, 2, 1).norm();
+
+ riemannFit::Vector2d dVec = hits.block(0, 0, 2, 1) + (-zInSZplane(0) + uVec(0)) * radii.block(0, 0, 2, 1);
+ riemannFit::Vector2d eVec = hits.block(0, 1, 2, 1) + (-zInSZplane(1) + uVec(1)) * radii.block(0, 1, 2, 1);
+
+ circle_results.par << atan2((eVec - dVec)(1), (eVec - dVec)(0)),
+ -circle_results.qCharge *
+ (fast_fit(2) - sqrt(riemannFit::sqr(fast_fit(2)) - 0.25 * (eVec - dVec).squaredNorm())),
+ circle_results.qCharge * (1. / fast_fit(2) + uVec(n));
+
+ assert(circle_results.qCharge * circle_results.par(1) <= 0);
+
+ riemannFit::Vector2d eMinusd = eVec - dVec;
+ double tmp1 = eMinusd.squaredNorm();
+ double tmp2 = sqrt(riemannFit::sqr(2 * fast_fit(2)) - tmp1);
+
+ riemannFit::Matrix3d jacobian;
+ jacobian << (radii(1, 0) * eMinusd(0) - eMinusd(1) * radii(0, 0)) / tmp1,
+ (radii(1, 1) * eMinusd(0) - eMinusd(1) * radii(0, 1)) / tmp1, 0,
+ (circle_results.qCharge / 2) * (eMinusd(0) * radii(0, 0) + eMinusd(1) * radii(1, 0)) / tmp2,
+ (circle_results.qCharge / 2) * (eMinusd(0) * radii(0, 1) + eMinusd(1) * radii(1, 1)) / tmp2, 0, 0, 0,
+ circle_results.qCharge;
+
+ circle_results.cov << iMat(0, 0), iMat(0, 1), iMat(0, n), iMat(1, 0), iMat(1, 1), iMat(1, n), iMat(n, 0),
+ iMat(n, 1), iMat(n, n);
+
+ circle_results.cov = jacobian * circle_results.cov * jacobian.transpose();
+
+ //...Translate in the system in which the first corrected hit is the origin, adding the m.s. correction...
+
+ auto eMinusDVec = eVec - dVec;
+ translateKarimaki(circle_results, 0.5 * eMinusDVec(0), 0.5 * eMinusDVec(1), jacobian);
+ circle_results.cov(0, 0) +=
+ (1 + riemannFit::sqr(slope)) * multScatt(sTotal(1) - sTotal(0), bField, fast_fit(2), 2, slope);
+
+ //...And translate back to the original system
+
+ translateKarimaki(circle_results, dVec(0), dVec(1), jacobian);
+
+ // compute chi2
+ circle_results.chi2 = 0;
+ for (u_int i = 0; i < n; i++) {
+ circle_results.chi2 += weightsVec(i) * riemannFit::sqr(zInSZplane(i) - uVec(i));
+ if (i > 0 && i < n - 1)
+ circle_results.chi2 +=
+ riemannFit::sqr(uVec(i - 1) / (sTransverse(i) - sTransverse(i - 1)) -
+ uVec(i) * (sTransverse(i + 1) - sTransverse(i - 1)) /
+ ((sTransverse(i + 1) - sTransverse(i)) * (sTransverse(i) - sTransverse(i - 1))) +
+ uVec(i + 1) / (sTransverse(i + 1) - sTransverse(i)) +
+ (sTransverse(i + 1) - sTransverse(i - 1)) * uVec(n) / 2) /
+ varBeta(i);
+ }
+
+ // assert(circle_results.chi2>=0);
+ }
+
+ /*!
+ \brief Performs the Broken Line fit in the straight track case (that is, the fit parameters are only the interceptions u).
+
+ \param hits hits coordinates.
+ \param fast_fit pre-fit result in the form (X0,Y0,R,tan(theta)).
+ \param bField magnetic field in Gev/cm/c.
+ \param data PreparedBrokenLineData.
+ \param line_results struct to be filled with the results in this form:
+ -par parameter of the line in this form: (cot(theta), Zip); \n
+ -cov covariance matrix of the fitted parameter; \n
+ -chi2 value of the cost function in the minimum.
+
+ \details The function implements the steps 2 and 3 of the Broken Line fit without
+ * the curvature correction.\n
+ * The step 2 is the least square fit, done by imposing the minimum constraint
+ * on the cost function and solving the consequent linear system. It determines
+ * the fitted parameters u and their covariance matrix.
+ * The step 3 is the correction of the fast pre-fitted parameters for the innermost
+ * part of the track. It is first done in a comfortable coordinate system (the one
+ * in which the first hit is the origin) and then the parameters and their covariance
+ * matrix are transformed to the original coordinate system.
+ */
+ template
+ __host__ __device__ inline void lineFit(const M6xN& hits_ge,
+ const V4& fast_fit,
+ const double bField,
+ const PreparedBrokenLineData& data,
+ riemannFit::LineFit& line_results) {
+ const auto& radii = data.radii;
+ const auto& sTotal = data.sTotal;
+ const auto& zInSZplane = data.zInSZplane;
+ const auto& varBeta = data.varBeta;
+
+ const double slope = -data.qCharge / fast_fit(3);
+ riemannFit::Matrix2d rotMat = rotationMatrix(slope);
+
+ riemannFit::Matrix3d vMat = riemannFit::Matrix3d::Zero(); // covariance matrix XYZ
+ riemannFit::Matrix2x3d jacobXYZtosZ =
+ riemannFit::Matrix2x3d::Zero(); // jacobian for computation of the error on s (xyz -> sz)
+ riemannFit::VectorNd weights = riemannFit::VectorNd::Zero();
+ for (u_int i = 0; i < n; i++) {
+ vMat(0, 0) = hits_ge.col(i)[0]; // x errors
+ vMat(0, 1) = vMat(1, 0) = hits_ge.col(i)[1]; // cov_xy
+ vMat(0, 2) = vMat(2, 0) = hits_ge.col(i)[3]; // cov_xz
+ vMat(1, 1) = hits_ge.col(i)[2]; // y errors
+ vMat(2, 1) = vMat(1, 2) = hits_ge.col(i)[4]; // cov_yz
+ vMat(2, 2) = hits_ge.col(i)[5]; // z errors
+ auto tmp = 1. / radii.block(0, i, 2, 1).norm();
+ jacobXYZtosZ(0, 0) = radii(1, i) * tmp;
+ jacobXYZtosZ(0, 1) = -radii(0, i) * tmp;
+ jacobXYZtosZ(1, 2) = 1.;
+ weights(i) = 1. / ((rotMat * jacobXYZtosZ * vMat * jacobXYZtosZ.transpose() * rotMat.transpose())(
+ 1, 1)); // compute the orthogonal weight point by point
+ }
+
+ riemannFit::VectorNd r_u;
+ for (u_int i = 0; i < n; i++) {
+ r_u(i) = weights(i) * zInSZplane(i);
+ }
+#ifdef CPP_DUMP
+ std::cout << "CU4\n" << matrixC_u(w, sTotal, varBeta) << std::endl;
+#endif
+ riemannFit::MatrixNd iMat;
+ math::cholesky::invert(matrixC_u(weights, sTotal, varBeta), iMat);
+#ifdef CPP_DUMP
+ std::cout << "I4\n" << iMat << std::endl;
+#endif
+
+ riemannFit::VectorNd uVec = iMat * r_u; // obtain the fitted parameters by solving the linear system
+
+ // line parameters in the system in which the first hit is the origin and with axis along SZ
+ line_results.par << (uVec(1) - uVec(0)) / (sTotal(1) - sTotal(0)), uVec(0);
+ auto idiff = 1. / (sTotal(1) - sTotal(0));
+ line_results.cov << (iMat(0, 0) - 2 * iMat(0, 1) + iMat(1, 1)) * riemannFit::sqr(idiff) +
+ multScatt(sTotal(1) - sTotal(0), bField, fast_fit(2), 2, slope),
+ (iMat(0, 1) - iMat(0, 0)) * idiff, (iMat(0, 1) - iMat(0, 0)) * idiff, iMat(0, 0);
+
+ // translate to the original SZ system
+ riemannFit::Matrix2d jacobian;
+ jacobian(0, 0) = 1.;
+ jacobian(0, 1) = 0;
+ jacobian(1, 0) = -sTotal(0);
+ jacobian(1, 1) = 1.;
+ line_results.par(1) += -line_results.par(0) * sTotal(0);
+ line_results.cov = jacobian * line_results.cov * jacobian.transpose();
+
+ // rotate to the original sz system
+ auto tmp = rotMat(0, 0) - line_results.par(0) * rotMat(0, 1);
+ jacobian(1, 1) = 1. / tmp;
+ jacobian(0, 0) = jacobian(1, 1) * jacobian(1, 1);
+ jacobian(0, 1) = 0;
+ jacobian(1, 0) = line_results.par(1) * rotMat(0, 1) * jacobian(0, 0);
+ line_results.par(1) = line_results.par(1) * jacobian(1, 1);
+ line_results.par(0) = (rotMat(0, 1) + line_results.par(0) * rotMat(0, 0)) * jacobian(1, 1);
+ line_results.cov = jacobian * line_results.cov * jacobian.transpose();
+
+ // compute chi2
+ line_results.chi2 = 0;
+ for (u_int i = 0; i < n; i++) {
+ line_results.chi2 += weights(i) * riemannFit::sqr(zInSZplane(i) - uVec(i));
+ if (i > 0 && i < n - 1)
+ line_results.chi2 += riemannFit::sqr(uVec(i - 1) / (sTotal(i) - sTotal(i - 1)) -
+ uVec(i) * (sTotal(i + 1) - sTotal(i - 1)) /
+ ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))) +
+ uVec(i + 1) / (sTotal(i + 1) - sTotal(i))) /
+ varBeta(i);
+ }
+ }
+
+ /*!
+ \brief Helix fit by three step:
+ -fast pre-fit (see Fast_fit() for further info); \n
+ -circle fit of the hits projected in the transverse plane by Broken Line algorithm (see BL_Circle_fit() for further info); \n
+ -line fit of the hits projected on the (pre-fitted) cilinder surface by Broken Line algorithm (see BL_Line_fit() for further info); \n
+ Points must be passed ordered (from inner to outer layer).
+
+ \param hits Matrix3xNd hits coordinates in this form: \n
+ |x1|x2|x3|...|xn| \n
+ |y1|y2|y3|...|yn| \n
+ |z1|z2|z3|...|zn|
+ \param hits_cov Matrix3Nd covariance matrix in this form (()->cov()): \n
+ |(x1,x1)|(x2,x1)|(x3,x1)|(x4,x1)|.|(y1,x1)|(y2,x1)|(y3,x1)|(y4,x1)|.|(z1,x1)|(z2,x1)|(z3,x1)|(z4,x1)| \n
+ |(x1,x2)|(x2,x2)|(x3,x2)|(x4,x2)|.|(y1,x2)|(y2,x2)|(y3,x2)|(y4,x2)|.|(z1,x2)|(z2,x2)|(z3,x2)|(z4,x2)| \n
+ |(x1,x3)|(x2,x3)|(x3,x3)|(x4,x3)|.|(y1,x3)|(y2,x3)|(y3,x3)|(y4,x3)|.|(z1,x3)|(z2,x3)|(z3,x3)|(z4,x3)| \n
+ |(x1,x4)|(x2,x4)|(x3,x4)|(x4,x4)|.|(y1,x4)|(y2,x4)|(y3,x4)|(y4,x4)|.|(z1,x4)|(z2,x4)|(z3,x4)|(z4,x4)| \n
+ . . . . . . . . . . . . . . . \n
+ |(x1,y1)|(x2,y1)|(x3,y1)|(x4,y1)|.|(y1,y1)|(y2,y1)|(y3,x1)|(y4,y1)|.|(z1,y1)|(z2,y1)|(z3,y1)|(z4,y1)| \n
+ |(x1,y2)|(x2,y2)|(x3,y2)|(x4,y2)|.|(y1,y2)|(y2,y2)|(y3,x2)|(y4,y2)|.|(z1,y2)|(z2,y2)|(z3,y2)|(z4,y2)| \n
+ |(x1,y3)|(x2,y3)|(x3,y3)|(x4,y3)|.|(y1,y3)|(y2,y3)|(y3,x3)|(y4,y3)|.|(z1,y3)|(z2,y3)|(z3,y3)|(z4,y3)| \n
+ |(x1,y4)|(x2,y4)|(x3,y4)|(x4,y4)|.|(y1,y4)|(y2,y4)|(y3,x4)|(y4,y4)|.|(z1,y4)|(z2,y4)|(z3,y4)|(z4,y4)| \n
+ . . . . . . . . . . . . . . . \n
+ |(x1,z1)|(x2,z1)|(x3,z1)|(x4,z1)|.|(y1,z1)|(y2,z1)|(y3,z1)|(y4,z1)|.|(z1,z1)|(z2,z1)|(z3,z1)|(z4,z1)| \n
+ |(x1,z2)|(x2,z2)|(x3,z2)|(x4,z2)|.|(y1,z2)|(y2,z2)|(y3,z2)|(y4,z2)|.|(z1,z2)|(z2,z2)|(z3,z2)|(z4,z2)| \n
+ |(x1,z3)|(x2,z3)|(x3,z3)|(x4,z3)|.|(y1,z3)|(y2,z3)|(y3,z3)|(y4,z3)|.|(z1,z3)|(z2,z3)|(z3,z3)|(z4,z3)| \n
+ |(x1,z4)|(x2,z4)|(x3,z4)|(x4,z4)|.|(y1,z4)|(y2,z4)|(y3,z4)|(y4,z4)|.|(z1,z4)|(z2,z4)|(z3,z4)|(z4,z4)|
+ \param bField magnetic field in the center of the detector in Gev/cm/c, in order to perform the p_t calculation.
+
+ \warning see BL_Circle_fit(), BL_Line_fit() and Fast_fit() warnings.
+
+ \bug see BL_Circle_fit(), BL_Line_fit() and Fast_fit() bugs.
+
+ \return (phi,Tip,p_t,cot(theta)),Zip), their covariance matrix and the chi2's of the circle and line fits.
+ */
+ template
+ inline riemannFit::HelixFit helixFit(const riemannFit::Matrix3xNd& hits,
+ const Eigen::Matrix& hits_ge,
+ const double bField) {
+ riemannFit::HelixFit helix;
+ riemannFit::Vector4d fast_fit;
+ fastFit(hits, fast_fit);
+
+ PreparedBrokenLineData data;
+ karimaki_circle_fit circle;
+ riemannFit::LineFit line;
+ riemannFit::Matrix3d jacobian;
+
+ prepareBrokenLineData(hits, fast_fit, bField, data);
+ lineFit(hits_ge, fast_fit, bField, data, line);
+ circleFit(hits, hits_ge, fast_fit, bField, data, circle);
+
+ // the circle fit gives k, but here we want p_t, so let's change the parameter and the covariance matrix
+ jacobian << 1., 0, 0, 0, 1., 0, 0, 0,
+ -std::abs(circle.par(2)) * bField / (riemannFit::sqr(circle.par(2)) * circle.par(2));
+ circle.par(2) = bField / std::abs(circle.par(2));
+ circle.cov = jacobian * circle.cov * jacobian.transpose();
+
+ helix.par << circle.par, line.par;
+ helix.cov = riemannFit::MatrixXd::Zero(5, 5);
+ helix.cov.block(0, 0, 3, 3) = circle.cov;
+ helix.cov.block(3, 3, 2, 2) = line.cov;
+ helix.qCharge = circle.qCharge;
+ helix.chi2_circle = circle.chi2;
+ helix.chi2_line = line.chi2;
+
+ return helix;
+ }
+
+} // namespace brokenline
+
+#endif // RecoPixelVertexing_PixelTrackFitting_interface_BrokenLine_h
diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/FitResult.h b/RecoPixelVertexing/PixelTrackFitting/interface/FitResult.h
new file mode 100644
index 0000000000000..01497719d2998
--- /dev/null
+++ b/RecoPixelVertexing/PixelTrackFitting/interface/FitResult.h
@@ -0,0 +1,65 @@
+#ifndef RecoPixelVertexing_PixelTrackFitting_interface_FitResult_h
+#define RecoPixelVertexing_PixelTrackFitting_interface_FitResult_h
+
+#include
+#include
+
+#include
+#include
+#include
+
+namespace riemannFit {
+
+ using Vector2d = Eigen::Vector2d;
+ using Vector3d = Eigen::Vector3d;
+ using Vector4d = Eigen::Vector4d;
+ using Vector5d = Eigen::Matrix;
+ using Matrix2d = Eigen::Matrix2d;
+ using Matrix3d = Eigen::Matrix3d;
+ using Matrix4d = Eigen::Matrix4d;
+ using Matrix5d = Eigen::Matrix;
+ using Matrix6d = Eigen::Matrix;
+
+ template
+ using Matrix3xNd = Eigen::Matrix; // used for inputs hits
+
+ struct CircleFit {
+ Vector3d par; //!< parameter: (X0,Y0,R)
+ Matrix3d cov;
+ /*!< covariance matrix: \n
+ |cov(X0,X0)|cov(Y0,X0)|cov( R,X0)| \n
+ |cov(X0,Y0)|cov(Y0,Y0)|cov( R,Y0)| \n
+ |cov(X0, R)|cov(Y0, R)|cov( R, R)|
+ */
+ int32_t qCharge; //!< particle charge
+ float chi2;
+ };
+
+ struct LineFit {
+ Vector2d par; //!<(cotan(theta),Zip)
+ Matrix2d cov;
+ /*!<
+ |cov(c_t,c_t)|cov(Zip,c_t)| \n
+ |cov(c_t,Zip)|cov(Zip,Zip)|
+ */
+ double chi2;
+ };
+
+ struct HelixFit {
+ Vector5d par; //!<(phi,Tip,pt,cotan(theta)),Zip)
+ Matrix5d cov;
+ /*!< ()->cov() \n
+ |(phi,phi)|(Tip,phi)|(p_t,phi)|(c_t,phi)|(Zip,phi)| \n
+ |(phi,Tip)|(Tip,Tip)|(p_t,Tip)|(c_t,Tip)|(Zip,Tip)| \n
+ |(phi,p_t)|(Tip,p_t)|(p_t,p_t)|(c_t,p_t)|(Zip,p_t)| \n
+ |(phi,c_t)|(Tip,c_t)|(p_t,c_t)|(c_t,c_t)|(Zip,c_t)| \n
+ |(phi,Zip)|(Tip,Zip)|(p_t,Zip)|(c_t,Zip)|(Zip,Zip)|
+ */
+ float chi2_circle;
+ float chi2_line;
+ // Vector4d fast_fit;
+ int32_t qCharge; //!< particle charge
+ }; // __attribute__((aligned(16)));
+
+} // namespace riemannFit
+#endif
diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h b/RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h
new file mode 100644
index 0000000000000..2fe74f53a7bd2
--- /dev/null
+++ b/RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h
@@ -0,0 +1,243 @@
+#ifndef RecoPixelVertexing_PixelTrackFitting_interface_FitUtils_h
+#define RecoPixelVertexing_PixelTrackFitting_interface_FitUtils_h
+
+#include "DataFormats/Math/interface/choleskyInversion.h"
+#include "HeterogeneousCore/CUDAUtilities/interface/cuda_assert.h"
+#include "RecoPixelVertexing/PixelTrackFitting/interface/FitResult.h"
+
+namespace riemannFit {
+
+ constexpr double epsilon = 1.e-4; //!< used in numerical derivative (J2 in Circle_fit())
+
+ using VectorXd = Eigen::VectorXd;
+ using MatrixXd = Eigen::MatrixXd;
+ template
+ using MatrixNd = Eigen::Matrix;
+ template
+ using MatrixNplusONEd = Eigen::Matrix;
+ template
+ using ArrayNd = Eigen::Array;
+ template
+ using Matrix2Nd = Eigen::Matrix;
+ template
+ using Matrix3Nd = Eigen::Matrix;
+ template
+ using Matrix2xNd = Eigen::Matrix;
+ template
+ using Array2xNd = Eigen::Array;
+ template
+ using MatrixNx3d = Eigen::Matrix;
+ template
+ using MatrixNx5d = Eigen::Matrix;
+ template
+ using VectorNd = Eigen::Matrix;
+ template
+ using VectorNplusONEd = Eigen::Matrix;
+ template
+ using Vector2Nd = Eigen::Matrix;
+ template
+ using Vector3Nd = Eigen::Matrix;
+ template
+ using RowVectorNd = Eigen::Matrix;
+ template
+ using RowVector2Nd = Eigen::Matrix;
+
+ using Matrix2x3d = Eigen::Matrix;
+
+ using Matrix3f = Eigen::Matrix3f;
+ using Vector3f = Eigen::Vector3f;
+ using Vector4f = Eigen::Vector4f;
+ using Vector6f = Eigen::Matrix;
+
+ template
+ __host__ __device__ void printIt(C* m, const char* prefix = "") {
+#ifdef RFIT_DEBUG
+ for (uint r = 0; r < m->rows(); ++r) {
+ for (uint c = 0; c < m->cols(); ++c) {
+ printf("%s Matrix(%d,%d) = %g\n", prefix, r, c, (*m)(r, c));
+ }
+ }
+#endif
+ }
+
+ /*!
+ \brief raise to square.
+ */
+ template
+ constexpr T sqr(const T a) {
+ return a * a;
+ }
+
+ /*!
+ \brief Compute cross product of two 2D vector (assuming z component 0),
+ returning z component of the result.
+ \param a first 2D vector in the product.
+ \param b second 2D vector in the product.
+ \return z component of the cross product.
+ */
+
+ __host__ __device__ inline double cross2D(const Vector2d& a, const Vector2d& b) {
+ return a.x() * b.y() - a.y() * b.x();
+ }
+
+ /*!
+ * load error in CMSSW format to our formalism
+ *
+ */
+ template
+ __host__ __device__ void loadCovariance2D(M6xNf const& ge, M2Nd& hits_cov) {
+ // Index numerology:
+ // i: index of the hits/point (0,..,3)
+ // j: index of space component (x,y,z)
+ // l: index of space components (x,y,z)
+ // ge is always in sync with the index i and is formatted as:
+ // ge[] ==> [xx, xy, yy, xz, yz, zz]
+ // in (j,l) notation, we have:
+ // ge[] ==> [(0,0), (0,1), (1,1), (0,2), (1,2), (2,2)]
+ // so the index ge_idx corresponds to the matrix elements:
+ // | 0 1 3 |
+ // | 1 2 4 |
+ // | 3 4 5 |
+ constexpr uint32_t hits_in_fit = M6xNf::ColsAtCompileTime;
+ for (uint32_t i = 0; i < hits_in_fit; ++i) {
+ {
+ constexpr uint32_t ge_idx = 0, j = 0, l = 0;
+ hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
+ }
+ {
+ constexpr uint32_t ge_idx = 2, j = 1, l = 1;
+ hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
+ }
+ {
+ constexpr uint32_t ge_idx = 1, j = 1, l = 0;
+ hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) = hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) =
+ ge.col(i)[ge_idx];
+ }
+ }
+ }
+
+ template
+ __host__ __device__ void loadCovariance(M6xNf const& ge, M3xNd& hits_cov) {
+ // Index numerology:
+ // i: index of the hits/point (0,..,3)
+ // j: index of space component (x,y,z)
+ // l: index of space components (x,y,z)
+ // ge is always in sync with the index i and is formatted as:
+ // ge[] ==> [xx, xy, yy, xz, yz, zz]
+ // in (j,l) notation, we have:
+ // ge[] ==> [(0,0), (0,1), (1,1), (0,2), (1,2), (2,2)]
+ // so the index ge_idx corresponds to the matrix elements:
+ // | 0 1 3 |
+ // | 1 2 4 |
+ // | 3 4 5 |
+ constexpr uint32_t hits_in_fit = M6xNf::ColsAtCompileTime;
+ for (uint32_t i = 0; i < hits_in_fit; ++i) {
+ {
+ constexpr uint32_t ge_idx = 0, j = 0, l = 0;
+ hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
+ }
+ {
+ constexpr uint32_t ge_idx = 2, j = 1, l = 1;
+ hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
+ }
+ {
+ constexpr uint32_t ge_idx = 5, j = 2, l = 2;
+ hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
+ }
+ {
+ constexpr uint32_t ge_idx = 1, j = 1, l = 0;
+ hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) = hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) =
+ ge.col(i)[ge_idx];
+ }
+ {
+ constexpr uint32_t ge_idx = 3, j = 2, l = 0;
+ hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) = hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) =
+ ge.col(i)[ge_idx];
+ }
+ {
+ constexpr uint32_t ge_idx = 4, j = 2, l = 1;
+ hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) = hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) =
+ ge.col(i)[ge_idx];
+ }
+ }
+ }
+
+ /*!
+ \brief Transform circle parameter from (X0,Y0,R) to (phi,Tip,p_t) and
+ consequently covariance matrix.
+ \param circle_uvr parameter (X0,Y0,R), covariance matrix to
+ be transformed and particle charge.
+ \param B magnetic field in Gev/cm/c unit.
+ \param error flag for errors computation.
+ */
+ __host__ __device__ inline void par_uvrtopak(CircleFit& circle, const double B, const bool error) {
+ Vector3d par_pak;
+ const double temp0 = circle.par.head(2).squaredNorm();
+ const double temp1 = sqrt(temp0);
+ par_pak << atan2(circle.qCharge * circle.par(0), -circle.qCharge * circle.par(1)),
+ circle.qCharge * (temp1 - circle.par(2)), circle.par(2) * B;
+ if (error) {
+ const double temp2 = sqr(circle.par(0)) * 1. / temp0;
+ const double temp3 = 1. / temp1 * circle.qCharge;
+ Matrix3d j4Mat;
+ j4Mat << -circle.par(1) * temp2 * 1. / sqr(circle.par(0)), temp2 * 1. / circle.par(0), 0., circle.par(0) * temp3,
+ circle.par(1) * temp3, -circle.qCharge, 0., 0., B;
+ circle.cov = j4Mat * circle.cov * j4Mat.transpose();
+ }
+ circle.par = par_pak;
+ }
+
+ /*!
+ \brief Transform circle parameter from (X0,Y0,R) to (phi,Tip,q/R) and
+ consequently covariance matrix.
+ \param circle_uvr parameter (X0,Y0,R), covariance matrix to
+ be transformed and particle charge.
+ */
+ __host__ __device__ inline void fromCircleToPerigee(CircleFit& circle) {
+ Vector3d par_pak;
+ const double temp0 = circle.par.head(2).squaredNorm();
+ const double temp1 = sqrt(temp0);
+ par_pak << atan2(circle.qCharge * circle.par(0), -circle.qCharge * circle.par(1)),
+ circle.qCharge * (temp1 - circle.par(2)), circle.qCharge / circle.par(2);
+
+ const double temp2 = sqr(circle.par(0)) * 1. / temp0;
+ const double temp3 = 1. / temp1 * circle.qCharge;
+ Matrix3d j4Mat;
+ j4Mat << -circle.par(1) * temp2 * 1. / sqr(circle.par(0)), temp2 * 1. / circle.par(0), 0., circle.par(0) * temp3,
+ circle.par(1) * temp3, -circle.qCharge, 0., 0., -circle.qCharge / (circle.par(2) * circle.par(2));
+ circle.cov = j4Mat * circle.cov * j4Mat.transpose();
+
+ circle.par = par_pak;
+ }
+
+ // transformation between the "perigee" to cmssw localcoord frame
+ // the plane of the latter is the perigee plane...
+ // from //!<(phi,Tip,q/pt,cotan(theta)),Zip)
+ // to q/p,dx/dz,dy/dz,x,z
+ template
+ __host__ __device__ inline void transformToPerigeePlane(VI5 const& ip, MI5 const& icov, VO5& op, MO5& ocov) {
+ auto sinTheta2 = 1. / (1. + ip(3) * ip(3));
+ auto sinTheta = std::sqrt(sinTheta2);
+ auto cosTheta = ip(3) * sinTheta;
+
+ op(0) = sinTheta * ip(2);
+ op(1) = 0.;
+ op(2) = -ip(3);
+ op(3) = ip(1);
+ op(4) = -ip(4);
+
+ Matrix5d jMat = Matrix5d::Zero();
+
+ jMat(0, 2) = sinTheta;
+ jMat(0, 3) = -sinTheta2 * cosTheta * ip(2);
+ jMat(1, 0) = 1.;
+ jMat(2, 3) = -1.;
+ jMat(3, 1) = 1.;
+ jMat(4, 4) = -1;
+
+ ocov = jMat * icov * jMat.transpose();
+ }
+
+} // namespace riemannFit
+
+#endif // RecoPixelVertexing_PixelTrackFitting_interface_FitUtils_h
diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/PixelNtupletsFitter.h b/RecoPixelVertexing/PixelTrackFitting/interface/PixelNtupletsFitter.h
new file mode 100644
index 0000000000000..9fb8843589669
--- /dev/null
+++ b/RecoPixelVertexing/PixelTrackFitting/interface/PixelNtupletsFitter.h
@@ -0,0 +1,27 @@
+#ifndef RecoPixelVertexing_PixelTrackFitting_interface_PixelNtupletsFitter_h
+#define RecoPixelVertexing_PixelTrackFitting_interface_PixelNtupletsFitter_h
+
+#include
+
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelFitterBase.h"
+#include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
+
+class PixelNtupletsFitter final : public PixelFitterBase {
+public:
+ explicit PixelNtupletsFitter(float nominalB, const MagneticField* field, bool useRiemannFit);
+ ~PixelNtupletsFitter() override = default;
+ std::unique_ptr run(const std::vector& hits,
+ const TrackingRegion& region,
+ const edm::EventSetup& setup) const override;
+
+private:
+ float nominalB_;
+ const MagneticField* field_;
+ bool useRiemannFit_;
+};
+
+#endif // RecoPixelVertexing_PixelTrackFitting_interface_PixelNtupletsFitter_h
diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h
new file mode 100644
index 0000000000000..52cf4b637fb37
--- /dev/null
+++ b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h
@@ -0,0 +1,1008 @@
+#ifndef RecoPixelVertexing_PixelTrackFitting_interface_RiemannFit_h
+#define RecoPixelVertexing_PixelTrackFitting_interface_RiemannFit_h
+
+#include "RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h"
+
+namespace riemannFit {
+
+ /*! Compute the Radiation length in the uniform hypothesis
+ *
+ * The Pixel detector, barrel and forward, is considered as an homogeneous
+ * cylinder of material, whose radiation lengths has been derived from the TDR
+ * plot that shows that 16cm correspond to 0.06 radiation lengths. Therefore
+ * one radiation length corresponds to 16cm/0.06 =~ 267 cm. All radiation
+ * lengths are computed using this unique number, in both regions, barrel and
+ * endcap.
+ *
+ * NB: no angle corrections nor projections are computed inside this routine.
+ * It is therefore the responsibility of the caller to supply the proper
+ * lengths in input. These lengths are the path traveled by the particle along
+ * its trajectory, namely the so called S of the helix in 3D space.
+ *
+ * \param length_values vector of incremental distances that will be translated
+ * into radiation length equivalent. Each radiation length i is computed
+ * incrementally with respect to the previous length i-1. The first length has
+ * no reference point (i.e. it has the dca).
+ *
+ * \return incremental radiation lengths that correspond to each segment.
+ */
+
+ template
+ __host__ __device__ inline void computeRadLenUniformMaterial(const VNd1& length_values, VNd2& rad_lengths) {
+ // Radiation length of the pixel detector in the uniform assumption, with
+ // 0.06 rad_len at 16 cm
+ constexpr double xx_0_inv = 0.06 / 16.;
+ uint n = length_values.rows();
+ rad_lengths(0) = length_values(0) * xx_0_inv;
+ for (uint j = 1; j < n; ++j) {
+ rad_lengths(j) = std::abs(length_values(j) - length_values(j - 1)) * xx_0_inv;
+ }
+ }
+
+ /*!
+ \brief Compute the covariance matrix along cartesian S-Z of points due to
+ multiple Coulomb scattering to be used in the line_fit, for the barrel
+ and forward cases.
+ The input covariance matrix is in the variables s-z, original and
+ unrotated.
+ The multiple scattering component is computed in the usual linear
+ approximation, using the 3D path which is computed as the squared root of
+ the squared sum of the s and z components passed in.
+ Internally a rotation by theta is performed and the covariance matrix
+ returned is the one in the direction orthogonal to the rotated S3D axis,
+ i.e. along the rotated Z axis.
+ The choice of the rotation is not arbitrary, but derived from the fact that
+ putting the horizontal axis along the S3D direction allows the usage of the
+ ordinary least squared fitting techiques with the trivial parametrization y
+ = mx + q, avoiding the patological case with m = +/- inf, that would
+ correspond to the case at eta = 0.
+ */
+
+ template
+ __host__ __device__ inline auto scatterCovLine(Matrix2d const* cov_sz,
+ const V4& fast_fit,
+ VNd1 const& s_arcs,
+ VNd2 const& z_values,
+ const double theta,
+ const double bField,
+ MatrixNd& ret) {
+#ifdef RFIT_DEBUG
+ riemannFit::printIt(&s_arcs, "Scatter_cov_line - s_arcs: ");
+#endif
+ constexpr uint n = N;
+ double p_t = std::min(20., fast_fit(2) * bField); // limit pt to avoid too small error!!!
+ double p_2 = p_t * p_t * (1. + 1. / sqr(fast_fit(3)));
+ VectorNd rad_lengths_S;
+ // See documentation at http://eigen.tuxfamily.org/dox/group__TutorialArrayClass.html
+ // Basically, to perform cwise operations on Matrices and Vectors, you need
+ // to transform them into Array-like objects.
+ VectorNd s_values = s_arcs.array() * s_arcs.array() + z_values.array() * z_values.array();
+ s_values = s_values.array().sqrt();
+ computeRadLenUniformMaterial(s_values, rad_lengths_S);
+ VectorNd sig2_S;
+ sig2_S = .000225 / p_2 * (1. + 0.038 * rad_lengths_S.array().log()).abs2() * rad_lengths_S.array();
+#ifdef RFIT_DEBUG
+ riemannFit::printIt(cov_sz, "Scatter_cov_line - cov_sz: ");
+#endif
+ Matrix2Nd tmp = Matrix2Nd::Zero();
+ for (uint k = 0; k < n; ++k) {
+ tmp(k, k) = cov_sz[k](0, 0);
+ tmp(k + n, k + n) = cov_sz[k](1, 1);
+ tmp(k, k + n) = tmp(k + n, k) = cov_sz[k](0, 1);
+ }
+ for (uint k = 0; k < n; ++k) {
+ for (uint l = k; l < n; ++l) {
+ for (uint i = 0; i < std::min(k, l); ++i) {
+ tmp(k + n, l + n) += std::abs(s_values(k) - s_values(i)) * std::abs(s_values(l) - s_values(i)) * sig2_S(i);
+ }
+ tmp(l + n, k + n) = tmp(k + n, l + n);
+ }
+ }
+ // We are interested only in the errors orthogonal to the rotated s-axis
+ // which, in our formalism, are in the lower square matrix.
+#ifdef RFIT_DEBUG
+ riemannFit::printIt(&tmp, "Scatter_cov_line - tmp: ");
+#endif
+ ret = tmp.block(n, n, n, n);
+ }
+
+ /*!
+ \brief Compute the covariance matrix (in radial coordinates) of points in
+ the transverse plane due to multiple Coulomb scattering.
+ \param p2D 2D points in the transverse plane.
+ \param fast_fit fast_fit Vector4d result of the previous pre-fit
+ structured in this form:(X0, Y0, R, Tan(Theta))).
+ \param B magnetic field use to compute p
+ \return scatter_cov_rad errors due to multiple scattering.
+ \warning input points must be ordered radially from the detector center
+ (from inner layer to outer ones; points on the same layer must ordered too).
+ \details Only the tangential component is computed (the radial one is
+ negligible).
+ */
+ template
+ __host__ __device__ inline MatrixNd scatter_cov_rad(const M2xN& p2D,
+ const V4& fast_fit,
+ VectorNd const& rad,
+ double B) {
+ constexpr uint n = N;
+ double p_t = std::min(20., fast_fit(2) * B); // limit pt to avoid too small error!!!
+ double p_2 = p_t * p_t * (1. + 1. / sqr(fast_fit(3)));
+ double theta = atan(fast_fit(3));
+ theta = theta < 0. ? theta + M_PI : theta;
+ VectorNd s_values;
+ VectorNd rad_lengths;
+ const Vector2d oVec(fast_fit(0), fast_fit(1));
+
+ // associated Jacobian, used in weights and errors computation
+ for (uint i = 0; i < n; ++i) { // x
+ Vector2d pVec = p2D.block(0, i, 2, 1) - oVec;
+ const double cross = cross2D(-oVec, pVec);
+ const double dot = (-oVec).dot(pVec);
+ const double tempAtan2 = atan2(cross, dot);
+ s_values(i) = std::abs(tempAtan2 * fast_fit(2));
+ }
+ computeRadLenUniformMaterial(s_values * sqrt(1. + 1. / sqr(fast_fit(3))), rad_lengths);
+ MatrixNd scatter_cov_rad = MatrixNd::Zero();
+ VectorNd sig2 = (1. + 0.038 * rad_lengths.array().log()).abs2() * rad_lengths.array();
+ sig2 *= 0.000225 / (p_2 * sqr(sin(theta)));
+ for (uint k = 0; k < n; ++k) {
+ for (uint l = k; l < n; ++l) {
+ for (uint i = 0; i < std::min(k, l); ++i) {
+ scatter_cov_rad(k, l) += (rad(k) - rad(i)) * (rad(l) - rad(i)) * sig2(i);
+ }
+ scatter_cov_rad(l, k) = scatter_cov_rad(k, l);
+ }
+ }
+#ifdef RFIT_DEBUG
+ riemannFit::printIt(&scatter_cov_rad, "Scatter_cov_rad - scatter_cov_rad: ");
+#endif
+ return scatter_cov_rad;
+ }
+
+ /*!
+ \brief Transform covariance matrix from radial (only tangential component)
+ to Cartesian coordinates (only transverse plane component).
+ \param p2D 2D points in the transverse plane.
+ \param cov_rad covariance matrix in radial coordinate.
+ \return cov_cart covariance matrix in Cartesian coordinates.
+*/
+
+ template
+ __host__ __device__ inline Matrix2Nd cov_radtocart(const M2xN& p2D,
+ const MatrixNd& cov_rad,
+ const VectorNd& rad) {
+#ifdef RFIT_DEBUG
+ printf("Address of p2D: %p\n", &p2D);
+#endif
+ printIt(&p2D, "cov_radtocart - p2D:");
+ constexpr uint n = N;
+ Matrix2Nd cov_cart = Matrix2Nd::Zero();
+ VectorNd rad_inv = rad.cwiseInverse();
+ printIt(&rad_inv, "cov_radtocart - rad_inv:");
+ for (uint i = 0; i < n; ++i) {
+ for (uint j = i; j < n; ++j) {
+ cov_cart(i, j) = cov_rad(i, j) * p2D(1, i) * rad_inv(i) * p2D(1, j) * rad_inv(j);
+ cov_cart(i + n, j + n) = cov_rad(i, j) * p2D(0, i) * rad_inv(i) * p2D(0, j) * rad_inv(j);
+ cov_cart(i, j + n) = -cov_rad(i, j) * p2D(1, i) * rad_inv(i) * p2D(0, j) * rad_inv(j);
+ cov_cart(i + n, j) = -cov_rad(i, j) * p2D(0, i) * rad_inv(i) * p2D(1, j) * rad_inv(j);
+ cov_cart(j, i) = cov_cart(i, j);
+ cov_cart(j + n, i + n) = cov_cart(i + n, j + n);
+ cov_cart(j + n, i) = cov_cart(i, j + n);
+ cov_cart(j, i + n) = cov_cart(i + n, j);
+ }
+ }
+ return cov_cart;
+ }
+
+ /*!
+ \brief Transform covariance matrix from Cartesian coordinates (only
+ transverse plane component) to radial coordinates (both radial and
+ tangential component but only diagonal terms, correlation between different
+ point are not managed).
+ \param p2D 2D points in transverse plane.
+ \param cov_cart covariance matrix in Cartesian coordinates.
+ \return cov_rad covariance matrix in raidal coordinate.
+ \warning correlation between different point are not computed.
+*/
+ template
+ __host__ __device__ inline VectorNd cov_carttorad(const M2xN& p2D,
+ const Matrix2Nd& cov_cart,
+ const VectorNd& rad) {
+ constexpr uint n = N;
+ VectorNd cov_rad;
+ const VectorNd rad_inv2 = rad.cwiseInverse().array().square();
+ for (uint i = 0; i < n; ++i) {
+ //!< in case you have (0,0) to avoid dividing by 0 radius
+ if (rad(i) < 1.e-4)
+ cov_rad(i) = cov_cart(i, i);
+ else {
+ cov_rad(i) = rad_inv2(i) * (cov_cart(i, i) * sqr(p2D(1, i)) + cov_cart(i + n, i + n) * sqr(p2D(0, i)) -
+ 2. * cov_cart(i, i + n) * p2D(0, i) * p2D(1, i));
+ }
+ }
+ return cov_rad;
+ }
+
+ /*!
+ \brief Transform covariance matrix from Cartesian coordinates (only
+ transverse plane component) to coordinates system orthogonal to the
+ pre-fitted circle in each point.
+ Further information in attached documentation.
+ \param p2D 2D points in transverse plane.
+ \param cov_cart covariance matrix in Cartesian coordinates.
+ \param fast_fit fast_fit Vector4d result of the previous pre-fit
+ structured in this form:(X0, Y0, R, tan(theta))).
+ \return cov_rad covariance matrix in the pre-fitted circle's
+ orthogonal system.
+*/
+ template
+ __host__ __device__ inline VectorNd cov_carttorad_prefit(const M2xN& p2D,
+ const Matrix2Nd& cov_cart,
+ V4& fast_fit,
+ const VectorNd& rad) {
+ constexpr uint n = N;
+ VectorNd cov_rad;
+ for (uint i = 0; i < n; ++i) {
+ //!< in case you have (0,0) to avoid dividing by 0 radius
+ if (rad(i) < 1.e-4)
+ cov_rad(i) = cov_cart(i, i); // TO FIX
+ else {
+ Vector2d a = p2D.col(i);
+ Vector2d b = p2D.col(i) - fast_fit.head(2);
+ const double x2 = a.dot(b);
+ const double y2 = cross2D(a, b);
+ const double tan_c = -y2 / x2;
+ const double tan_c2 = sqr(tan_c);
+ cov_rad(i) =
+ 1. / (1. + tan_c2) * (cov_cart(i, i) + cov_cart(i + n, i + n) * tan_c2 + 2 * cov_cart(i, i + n) * tan_c);
+ }
+ }
+ return cov_rad;
+ }
+
+ /*!
+ \brief Compute the points' weights' vector for the circle fit when multiple
+ scattering is managed.
+ Further information in attached documentation.
+ \param cov_rad_inv covariance matrix inverse in radial coordinated
+ (or, beter, pre-fitted circle's orthogonal system).
+ \return weight VectorNd points' weights' vector.
+ \bug I'm not sure this is the right way to compute the weights for non
+ diagonal cov matrix. Further investigation needed.
+*/
+
+ template
+ __host__ __device__ inline VectorNd weightCircle(const MatrixNd& cov_rad_inv) {
+ return cov_rad_inv.colwise().sum().transpose();
+ }
+
+ /*!
+ \brief Find particle q considering the sign of cross product between
+ particles velocity (estimated by the first 2 hits) and the vector radius
+ between the first hit and the center of the fitted circle.
+ \param p2D 2D points in transverse plane.
+ \param par_uvr result of the circle fit in this form: (X0,Y0,R).
+ \return q int 1 or -1.
+*/
+ template
+ __host__ __device__ inline int32_t charge(const M2xN& p2D, const Vector3d& par_uvr) {
+ return ((p2D(0, 1) - p2D(0, 0)) * (par_uvr.y() - p2D(1, 0)) - (p2D(1, 1) - p2D(1, 0)) * (par_uvr.x() - p2D(0, 0)) >
+ 0)
+ ? -1
+ : 1;
+ }
+
+ /*!
+ \brief Compute the eigenvector associated to the minimum eigenvalue.
+ \param A the Matrix you want to know eigenvector and eigenvalue.
+ \param chi2 the double were the chi2-related quantity will be stored.
+ \return the eigenvector associated to the minimum eigenvalue.
+ \warning double precision is needed for a correct assessment of chi2.
+ \details The minimus eigenvalue is related to chi2.
+ We exploit the fact that the matrix is symmetrical and small (2x2 for line
+ fit and 3x3 for circle fit), so the SelfAdjointEigenSolver from Eigen
+ library is used, with the computedDirect method (available only for 2x2
+ and 3x3 Matrix) wich computes eigendecomposition of given matrix using a
+ fast closed-form algorithm.
+ For this optimization the matrix type must be known at compiling time.
+*/
+
+ __host__ __device__ inline Vector3d min_eigen3D(const Matrix3d& A, double& chi2) {
+#ifdef RFIT_DEBUG
+ printf("min_eigen3D - enter\n");
+#endif
+ Eigen::SelfAdjointEigenSolver solver(3);
+ solver.computeDirect(A);
+ int min_index;
+ chi2 = solver.eigenvalues().minCoeff(&min_index);
+#ifdef RFIT_DEBUG
+ printf("min_eigen3D - exit\n");
+#endif
+ return solver.eigenvectors().col(min_index);
+ }
+
+ /*!
+ \brief A faster version of min_eigen3D() where double precision is not
+ needed.
+ \param A the Matrix you want to know eigenvector and eigenvalue.
+ \param chi2 the double were the chi2-related quantity will be stored
+ \return the eigenvector associated to the minimum eigenvalue.
+ \detail The computedDirect() method of SelfAdjointEigenSolver for 3x3 Matrix
+ indeed, use trigonometry function (it solves a third degree equation) which
+ speed up in single precision.
+*/
+
+ __host__ __device__ inline Vector3d min_eigen3D_fast(const Matrix3d& A) {
+ Eigen::SelfAdjointEigenSolver solver(3);
+ solver.computeDirect(A.cast());
+ int min_index;
+ solver.eigenvalues().minCoeff(&min_index);
+ return solver.eigenvectors().col(min_index).cast();
+ }
+
+ /*!
+ \brief 2D version of min_eigen3D().
+ \param aMat the Matrix you want to know eigenvector and eigenvalue.
+ \param chi2 the double were the chi2-related quantity will be stored
+ \return the eigenvector associated to the minimum eigenvalue.
+ \detail The computedDirect() method of SelfAdjointEigenSolver for 2x2 Matrix
+ do not use special math function (just sqrt) therefore it doesn't speed up
+ significantly in single precision.
+*/
+
+ __host__ __device__ inline Vector2d min_eigen2D(const Matrix2d& aMat, double& chi2) {
+ Eigen::SelfAdjointEigenSolver solver(2);
+ solver.computeDirect(aMat);
+ int min_index;
+ chi2 = solver.eigenvalues().minCoeff(&min_index);
+ return solver.eigenvectors().col(min_index);
+ }
+
+ /*!
+ \brief A very fast helix fit: it fits a circle by three points (first, middle
+ and last point) and a line by two points (first and last).
+ \param hits points to be fitted
+ \return result in this form: (X0,Y0,R,tan(theta)).
+ \warning points must be passed ordered (from internal layer to external) in
+ order to maximize accuracy and do not mistake tan(theta) sign.
+ \details This fast fit is used as pre-fit which is needed for:
+ - weights estimation and chi2 computation in line fit (fundamental);
+ - weights estimation and chi2 computation in circle fit (useful);
+ - computation of error due to multiple scattering.
+*/
+
+ template
+ __host__ __device__ inline void fastFit(const M3xN& hits, V4& result) {
+ constexpr uint32_t N = M3xN::ColsAtCompileTime;
+ constexpr auto n = N; // get the number of hits
+ printIt(&hits, "Fast_fit - hits: ");
+
+ // CIRCLE FIT
+ // Make segments between middle-to-first(b) and last-to-first(c) hits
+ const Vector2d bVec = hits.block(0, n / 2, 2, 1) - hits.block(0, 0, 2, 1);
+ const Vector2d cVec = hits.block(0, n - 1, 2, 1) - hits.block(0, 0, 2, 1);
+ printIt(&bVec, "Fast_fit - b: ");
+ printIt(&cVec, "Fast_fit - c: ");
+ // Compute their lengths
+ auto b2 = bVec.squaredNorm();
+ auto c2 = cVec.squaredNorm();
+ // The algebra has been verified (MR). The usual approach has been followed:
+ // * use an orthogonal reference frame passing from the first point.
+ // * build the segments (chords)
+ // * build orthogonal lines through mid points
+ // * make a system and solve for X0 and Y0.
+ // * add the initial point
+ bool flip = abs(bVec.x()) < abs(bVec.y());
+ auto bx = flip ? bVec.y() : bVec.x();
+ auto by = flip ? bVec.x() : bVec.y();
+ auto cx = flip ? cVec.y() : cVec.x();
+ auto cy = flip ? cVec.x() : cVec.y();
+ //!< in case b.x is 0 (2 hits with same x)
+ auto div = 2. * (cx * by - bx * cy);
+ // if aligned TO FIX
+ auto y0 = (cx * b2 - bx * c2) / div;
+ auto x0 = (0.5 * b2 - y0 * by) / bx;
+ result(0) = hits(0, 0) + (flip ? y0 : x0);
+ result(1) = hits(1, 0) + (flip ? x0 : y0);
+ result(2) = sqrt(sqr(x0) + sqr(y0));
+ printIt(&result, "Fast_fit - result: ");
+
+ // LINE FIT
+ const Vector2d dVec = hits.block(0, 0, 2, 1) - result.head(2);
+ const Vector2d eVec = hits.block(0, n - 1, 2, 1) - result.head(2);
+ printIt(&eVec, "Fast_fit - e: ");
+ printIt(&dVec, "Fast_fit - d: ");
+ // Compute the arc-length between first and last point: L = R * theta = R * atan (tan (Theta) )
+ auto dr = result(2) * atan2(cross2D(dVec, eVec), dVec.dot(eVec));
+ // Simple difference in Z between last and first hit
+ auto dz = hits(2, n - 1) - hits(2, 0);
+
+ result(3) = (dr / dz);
+
+#ifdef RFIT_DEBUG
+ printf("Fast_fit: [%f, %f, %f, %f]\n", result(0), result(1), result(2), result(3));
+#endif
+ }
+
+ /*!
+ \brief Fit a generic number of 2D points with a circle using Riemann-Chernov
+ algorithm. Covariance matrix of fitted parameter is optionally computed.
+ Multiple scattering (currently only in barrel layer) is optionally handled.
+ \param hits2D 2D points to be fitted.
+ \param hits_cov2D covariance matrix of 2D points.
+ \param fast_fit pre-fit result in this form: (X0,Y0,R,tan(theta)).
+ (tan(theta) is not used).
+ \param bField magnetic field
+ \param error flag for error computation.
+ \param scattering flag for multiple scattering
+ \return circle circle_fit:
+ -par parameter of the fitted circle in this form (X0,Y0,R); \n
+ -cov covariance matrix of the fitted parameter (not initialized if
+ error = false); \n
+ -q charge of the particle; \n
+ -chi2.
+ \warning hits must be passed ordered from inner to outer layer (double hits
+ on the same layer must be ordered too) so that multiple scattering is
+ treated properly.
+ \warning Multiple scattering for barrel is still not tested.
+ \warning Multiple scattering for endcap hits is not handled (yet). Do not
+ fit endcap hits with scattering = true !
+ \bug for small pt (<0.3 Gev/c) chi2 could be slightly underestimated.
+ \bug further investigation needed for error propagation with multiple
+ scattering.
+*/
+ template
+ __host__ __device__ inline CircleFit circleFit(const M2xN& hits2D,
+ const Matrix2Nd& hits_cov2D,
+ const V4& fast_fit,
+ const VectorNd& rad,
+ const double bField,
+ const bool error) {
+#ifdef RFIT_DEBUG
+ printf("circle_fit - enter\n");
+#endif
+ // INITIALIZATION
+ Matrix2Nd vMat = hits_cov2D;
+ constexpr uint n = N;
+ printIt(&hits2D, "circle_fit - hits2D:");
+ printIt(&hits_cov2D, "circle_fit - hits_cov2D:");
+
+#ifdef RFIT_DEBUG
+ printf("circle_fit - WEIGHT COMPUTATION\n");
+#endif
+ // WEIGHT COMPUTATION
+ VectorNd weight;
+ MatrixNd gMat;
+ double renorm;
+ {
+ MatrixNd cov_rad = cov_carttorad_prefit(hits2D, vMat, fast_fit, rad).asDiagonal();
+ MatrixNd scatterCovRadMat = scatter_cov_rad(hits2D, fast_fit, rad, bField);
+ printIt(&scatterCovRadMat, "circle_fit - scatter_cov_rad:");
+ printIt(&hits2D, "circle_fit - hits2D bis:");
+#ifdef RFIT_DEBUG
+ printf("Address of hits2D: a) %p\n", &hits2D);
+#endif
+ vMat += cov_radtocart(hits2D, scatterCovRadMat, rad);
+ printIt(&vMat, "circle_fit - V:");
+ cov_rad += scatterCovRadMat;
+ printIt(&cov_rad, "circle_fit - cov_rad:");
+ math::cholesky::invert(cov_rad, gMat);
+ // gMat = cov_rad.inverse();
+ renorm = gMat.sum();
+ gMat *= 1. / renorm;
+ weight = weightCircle(gMat);
+ }
+ printIt(&weight, "circle_fit - weight:");
+
+ // SPACE TRANSFORMATION
+#ifdef RFIT_DEBUG
+ printf("circle_fit - SPACE TRANSFORMATION\n");
+#endif
+
+ // center
+#ifdef RFIT_DEBUG
+ printf("Address of hits2D: b) %p\n", &hits2D);
+#endif
+ const Vector2d hCentroid = hits2D.rowwise().mean(); // centroid
+ printIt(&hCentroid, "circle_fit - h_:");
+ Matrix3xNd p3D;
+ p3D.block(0, 0, 2, n) = hits2D.colwise() - hCentroid;
+ printIt(&p3D, "circle_fit - p3D: a)");
+ Vector2Nd mc; // centered hits, used in error computation
+ mc << p3D.row(0).transpose(), p3D.row(1).transpose();
+ printIt(&mc, "circle_fit - mc(centered hits):");
+
+ // scale
+ const double tempQ = mc.squaredNorm();
+ const double tempS = sqrt(n * 1. / tempQ); // scaling factor
+ p3D *= tempS;
+
+ // project on paraboloid
+ p3D.row(2) = p3D.block(0, 0, 2, n).colwise().squaredNorm();
+ printIt(&p3D, "circle_fit - p3D: b)");
+
+#ifdef RFIT_DEBUG
+ printf("circle_fit - COST FUNCTION\n");
+#endif
+ // COST FUNCTION
+
+ // compute
+ Vector3d r0;
+ r0.noalias() = p3D * weight; // center of gravity
+ const Matrix3xNd xMat = p3D.colwise() - r0;
+ Matrix3d aMat = xMat * gMat * xMat.transpose();
+ printIt(&aMat, "circle_fit - A:");
+
+#ifdef RFIT_DEBUG
+ printf("circle_fit - MINIMIZE\n");
+#endif
+ // minimize
+ double chi2;
+ Vector3d vVec = min_eigen3D(aMat, chi2);
+#ifdef RFIT_DEBUG
+ printf("circle_fit - AFTER MIN_EIGEN\n");
+#endif
+ printIt(&vVec, "v BEFORE INVERSION");
+ vVec *= (vVec(2) > 0) ? 1 : -1; // TO FIX dovrebbe essere N(3)>0
+ printIt(&vVec, "v AFTER INVERSION");
+ // This hack to be able to run on GPU where the automatic assignment to a
+ // double from the vector multiplication is not working.
+#ifdef RFIT_DEBUG
+ printf("circle_fit - AFTER MIN_EIGEN 1\n");
+#endif
+ Eigen::Matrix cm;
+#ifdef RFIT_DEBUG
+ printf("circle_fit - AFTER MIN_EIGEN 2\n");
+#endif
+ cm = -vVec.transpose() * r0;
+#ifdef RFIT_DEBUG
+ printf("circle_fit - AFTER MIN_EIGEN 3\n");
+#endif
+ const double tempC = cm(0, 0);
+
+#ifdef RFIT_DEBUG
+ printf("circle_fit - COMPUTE CIRCLE PARAMETER\n");
+#endif
+ // COMPUTE CIRCLE PARAMETER
+
+ // auxiliary quantities
+ const double tempH = sqrt(1. - sqr(vVec(2)) - 4. * tempC * vVec(2));
+ const double v2x2_inv = 1. / (2. * vVec(2));
+ const double s_inv = 1. / tempS;
+ Vector3d par_uvr; // used in error propagation
+ par_uvr << -vVec(0) * v2x2_inv, -vVec(1) * v2x2_inv, tempH * v2x2_inv;
+
+ CircleFit circle;
+ circle.par << par_uvr(0) * s_inv + hCentroid(0), par_uvr(1) * s_inv + hCentroid(1), par_uvr(2) * s_inv;
+ circle.qCharge = charge(hits2D, circle.par);
+ circle.chi2 = abs(chi2) * renorm / sqr(2 * vVec(2) * par_uvr(2) * tempS);
+ printIt(&circle.par, "circle_fit - CIRCLE PARAMETERS:");
+ printIt(&circle.cov, "circle_fit - CIRCLE COVARIANCE:");
+#ifdef RFIT_DEBUG
+ printf("circle_fit - CIRCLE CHARGE: %d\n", circle.qCharge);
+#endif
+
+#ifdef RFIT_DEBUG
+ printf("circle_fit - ERROR PROPAGATION\n");
+#endif
+ // ERROR PROPAGATION
+ if (error) {
+#ifdef RFIT_DEBUG
+ printf("circle_fit - ERROR PRPAGATION ACTIVATED\n");
+#endif
+ ArrayNd vcsMat[2][2]; // cov matrix of center & scaled points
+ MatrixNd cMat[3][3]; // cov matrix of 3D transformed points
+#ifdef RFIT_DEBUG
+ printf("circle_fit - ERROR PRPAGATION ACTIVATED 2\n");
+#endif
+ {
+ Eigen::Matrix cm;
+ Eigen::Matrix cm2;
+ cm = mc.transpose() * vMat * mc;
+ const double tempC2 = cm(0, 0);
+ Matrix2Nd tempVcsMat;
+ tempVcsMat.template triangularView() =
+ (sqr(tempS) * vMat + sqr(sqr(tempS)) * 1. / (4. * tempQ * n) *
+ (2. * vMat.squaredNorm() + 4. * tempC2) * // mc.transpose() * V * mc) *
+ (mc * mc.transpose()));
+
+ printIt(&tempVcsMat, "circle_fit - Vcs:");
+ cMat[0][0] = tempVcsMat.block(0, 0, n, n).template selfadjointView();
+ vcsMat[0][1] = tempVcsMat.block(0, n, n, n);
+ cMat[1][1] = tempVcsMat.block(n, n, n, n).template selfadjointView();
+ vcsMat[1][0] = vcsMat[0][1].transpose();
+ printIt(&tempVcsMat, "circle_fit - Vcs:");
+ }
+
+ {
+ const ArrayNd t0 = (VectorXd::Constant(n, 1.) * p3D.row(0));
+ const ArrayNd t1 = (VectorXd::Constant(n, 1.) * p3D.row(1));
+ const ArrayNd t00 = p3D.row(0).transpose() * p3D.row(0);
+ const ArrayNd t01 = p3D.row(0).transpose() * p3D.row(1);
+ const ArrayNd t11 = p3D.row(1).transpose() * p3D.row(1);
+ const ArrayNd t10 = t01.transpose();
+ vcsMat[0][0] = cMat[0][0];
+ cMat[0][1] = vcsMat[0][1];
+ cMat[0][2] = 2. * (vcsMat[0][0] * t0 + vcsMat[0][1] * t1);
+ vcsMat[1][1] = cMat[1][1];
+ cMat[1][2] = 2. * (vcsMat[1][0] * t0 + vcsMat[1][1] * t1);
+ MatrixNd tmp;
+ tmp.template triangularView() =
+ (2. * (vcsMat[0][0] * vcsMat[0][0] + vcsMat[0][0] * vcsMat[0][1] + vcsMat[1][1] * vcsMat[1][0] +
+ vcsMat[1][1] * vcsMat[1][1]) +
+ 4. * (vcsMat[0][0] * t00 + vcsMat[0][1] * t01 + vcsMat[1][0] * t10 + vcsMat[1][1] * t11))
+ .matrix();
+ cMat[2][2] = tmp.template selfadjointView();
+ }
+ printIt(&cMat[0][0], "circle_fit - C[0][0]:");
+
+ Matrix3d c0Mat; // cov matrix of center of gravity (r0.x,r0.y,r0.z)
+ for (uint i = 0; i < 3; ++i) {
+ for (uint j = i; j < 3; ++j) {
+ Eigen::Matrix tmp;
+ tmp = weight.transpose() * cMat[i][j] * weight;
+ // Workaround to get things working in GPU
+ const double tempC = tmp(0, 0);
+ c0Mat(i, j) = tempC; //weight.transpose() * C[i][j] * weight;
+ c0Mat(j, i) = c0Mat(i, j);
+ }
+ }
+ printIt(&c0Mat, "circle_fit - C0:");
+
+ const MatrixNd wMat = weight * weight.transpose();
+ const MatrixNd hMat = MatrixNd::Identity().rowwise() - weight.transpose();
+ const MatrixNx3d s_v = hMat * p3D.transpose();
+ printIt(&wMat, "circle_fit - W:");
+ printIt(&hMat, "circle_fit - H:");
+ printIt(&s_v, "circle_fit - s_v:");
+
+ MatrixNd dMat[3][3]; // cov(s_v)
+ dMat[0][0] = (hMat * cMat[0][0] * hMat.transpose()).cwiseProduct(wMat);
+ dMat[0][1] = (hMat * cMat[0][1] * hMat.transpose()).cwiseProduct(wMat);
+ dMat[0][2] = (hMat * cMat[0][2] * hMat.transpose()).cwiseProduct(wMat);
+ dMat[1][1] = (hMat * cMat[1][1] * hMat.transpose()).cwiseProduct(wMat);
+ dMat[1][2] = (hMat * cMat[1][2] * hMat.transpose()).cwiseProduct(wMat);
+ dMat[2][2] = (hMat * cMat[2][2] * hMat.transpose()).cwiseProduct(wMat);
+ dMat[1][0] = dMat[0][1].transpose();
+ dMat[2][0] = dMat[0][2].transpose();
+ dMat[2][1] = dMat[1][2].transpose();
+ printIt(&dMat[0][0], "circle_fit - D_[0][0]:");
+
+ constexpr uint nu[6][2] = {{0, 0}, {0, 1}, {0, 2}, {1, 1}, {1, 2}, {2, 2}};
+
+ Matrix6d eMat; // cov matrix of the 6 independent elements of A
+ for (uint a = 0; a < 6; ++a) {
+ const uint i = nu[a][0], j = nu[a][1];
+ for (uint b = a; b < 6; ++b) {
+ const uint k = nu[b][0], l = nu[b][1];
+ VectorNd t0(n);
+ VectorNd t1(n);
+ if (l == k) {
+ t0 = 2. * dMat[j][l] * s_v.col(l);
+ if (i == j)
+ t1 = t0;
+ else
+ t1 = 2. * dMat[i][l] * s_v.col(l);
+ } else {
+ t0 = dMat[j][l] * s_v.col(k) + dMat[j][k] * s_v.col(l);
+ if (i == j)
+ t1 = t0;
+ else
+ t1 = dMat[i][l] * s_v.col(k) + dMat[i][k] * s_v.col(l);
+ }
+
+ if (i == j) {
+ Eigen::Matrix cm;
+ cm = s_v.col(i).transpose() * (t0 + t1);
+ // Workaround to get things working in GPU
+ const double tempC = cm(0, 0);
+ eMat(a, b) = 0. + tempC;
+ } else {
+ Eigen::Matrix cm;
+ cm = (s_v.col(i).transpose() * t0) + (s_v.col(j).transpose() * t1);
+ // Workaround to get things working in GPU
+ const double tempC = cm(0, 0);
+ eMat(a, b) = 0. + tempC; //(s_v.col(i).transpose() * t0) + (s_v.col(j).transpose() * t1);
+ }
+ if (b != a)
+ eMat(b, a) = eMat(a, b);
+ }
+ }
+ printIt(&eMat, "circle_fit - E:");
+
+ Eigen::Matrix j2Mat; // Jacobian of min_eigen() (numerically computed)
+ for (uint a = 0; a < 6; ++a) {
+ const uint i = nu[a][0], j = nu[a][1];
+ Matrix3d delta = Matrix3d::Zero();
+ delta(i, j) = delta(j, i) = abs(aMat(i, j) * epsilon);
+ j2Mat.col(a) = min_eigen3D_fast(aMat + delta);
+ const int sign = (j2Mat.col(a)(2) > 0) ? 1 : -1;
+ j2Mat.col(a) = (j2Mat.col(a) * sign - vVec) / delta(i, j);
+ }
+ printIt(&j2Mat, "circle_fit - J2:");
+
+ Matrix4d cvcMat; // joint cov matrix of (v0,v1,v2,c)
+ {
+ Matrix3d t0 = j2Mat * eMat * j2Mat.transpose();
+ Vector3d t1 = -t0 * r0;
+ cvcMat.block(0, 0, 3, 3) = t0;
+ cvcMat.block(0, 3, 3, 1) = t1;
+ cvcMat.block(3, 0, 1, 3) = t1.transpose();
+ Eigen::Matrix cm1;
+ Eigen::Matrix cm3;
+ cm1 = (vVec.transpose() * c0Mat * vVec);
+ // cm2 = (c0Mat.cwiseProduct(t0)).sum();
+ cm3 = (r0.transpose() * t0 * r0);
+ // Workaround to get things working in GPU
+ const double tempC = cm1(0, 0) + (c0Mat.cwiseProduct(t0)).sum() + cm3(0, 0);
+ cvcMat(3, 3) = tempC;
+ // (v.transpose() * c0Mat * v) + (c0Mat.cwiseProduct(t0)).sum() + (r0.transpose() * t0 * r0);
+ }
+ printIt(&cvcMat, "circle_fit - Cvc:");
+
+ Eigen::Matrix j3Mat; // Jacobian (v0,v1,v2,c)->(X0,Y0,R)
+ {
+ const double t = 1. / tempH;
+ j3Mat << -v2x2_inv, 0, vVec(0) * sqr(v2x2_inv) * 2., 0, 0, -v2x2_inv, vVec(1) * sqr(v2x2_inv) * 2., 0,
+ vVec(0) * v2x2_inv * t, vVec(1) * v2x2_inv * t,
+ -tempH * sqr(v2x2_inv) * 2. - (2. * tempC + vVec(2)) * v2x2_inv * t, -t;
+ }
+ printIt(&j3Mat, "circle_fit - J3:");
+
+ const RowVector2Nd Jq = mc.transpose() * tempS * 1. / n; // var(q)
+ printIt(&Jq, "circle_fit - Jq:");
+
+ Matrix3d cov_uvr = j3Mat * cvcMat * j3Mat.transpose() * sqr(s_inv) // cov(X0,Y0,R)
+ + (par_uvr * par_uvr.transpose()) * (Jq * vMat * Jq.transpose());
+
+ circle.cov = cov_uvr;
+ }
+
+ printIt(&circle.cov, "Circle cov:");
+#ifdef RFIT_DEBUG
+ printf("circle_fit - exit\n");
+#endif
+ return circle;
+ }
+
+ /*! \brief Perform an ordinary least square fit in the s-z plane to compute
+ * the parameters cotTheta and Zip.
+ *
+ * The fit is performed in the rotated S3D-Z' plane, following the formalism of
+ * Frodesen, Chapter 10, p. 259.
+ *
+ * The system has been rotated to both try to use the combined errors in s-z
+ * along Z', as errors in the Y direction and to avoid the patological case of
+ * degenerate lines with angular coefficient m = +/- inf.
+ *
+ * The rotation is using the information on the theta angle computed in the
+ * fast fit. The rotation is such that the S3D axis will be the X-direction,
+ * while the rotated Z-axis will be the Y-direction. This pretty much follows
+ * what is done in the same fit in the Broken Line approach.
+ */
+
+ template
+ __host__ __device__ inline LineFit lineFit(const M3xN& hits,
+ const M6xN& hits_ge,
+ const CircleFit& circle,
+ const V4& fast_fit,
+ const double bField,
+ const bool error) {
+ constexpr uint32_t N = M3xN::ColsAtCompileTime;
+ constexpr auto n = N;
+ double theta = -circle.qCharge * atan(fast_fit(3));
+ theta = theta < 0. ? theta + M_PI : theta;
+
+ // Prepare the Rotation Matrix to rotate the points
+ Eigen::Matrix rot;
+ rot << sin(theta), cos(theta), -cos(theta), sin(theta);
+
+ // PROJECTION ON THE CILINDER
+ //
+ // p2D will be:
+ // [s1, s2, s3, ..., sn]
+ // [z1, z2, z3, ..., zn]
+ // s values will be ordinary x-values
+ // z values will be ordinary y-values
+
+ Matrix2xNd p2D = Matrix2xNd::Zero();
+ Eigen::Matrix jxMat;
+
+#ifdef RFIT_DEBUG
+ printf("Line_fit - B: %g\n", bField);
+ printIt(&hits, "Line_fit points: ");
+ printIt(&hits_ge, "Line_fit covs: ");
+ printIt(&rot, "Line_fit rot: ");
+#endif
+ // x & associated Jacobian
+ // cfr https://indico.cern.ch/event/663159/contributions/2707659/attachments/1517175/2368189/Riemann_fit.pdf
+ // Slide 11
+ // a ==> -o i.e. the origin of the circle in XY plane, negative
+ // b ==> p i.e. distances of the points wrt the origin of the circle.
+ const Vector2d oVec(circle.par(0), circle.par(1));
+
+ // associated Jacobian, used in weights and errors computation
+ Matrix6d covMat = Matrix6d::Zero();
+ Matrix2d cov_sz[N];
+ for (uint i = 0; i < n; ++i) {
+ Vector2d pVec = hits.block(0, i, 2, 1) - oVec;
+ const double cross = cross2D(-oVec, pVec);
+ const double dot = (-oVec).dot(pVec);
+ // atan2(cross, dot) give back the angle in the transverse plane so tha the
+ // final equation reads: x_i = -q*R*theta (theta = angle returned by atan2)
+ const double tempQAtan2 = -circle.qCharge * atan2(cross, dot);
+ // p2D.coeffRef(1, i) = atan2_ * circle.par(2);
+ p2D(0, i) = tempQAtan2 * circle.par(2);
+
+ // associated Jacobian, used in weights and errors- computation
+ const double temp0 = -circle.qCharge * circle.par(2) * 1. / (sqr(dot) + sqr(cross));
+ double d_X0 = 0., d_Y0 = 0., d_R = 0.; // good approximation for big pt and eta
+ if (error) {
+ d_X0 = -temp0 * ((pVec(1) + oVec(1)) * dot - (pVec(0) - oVec(0)) * cross);
+ d_Y0 = temp0 * ((pVec(0) + oVec(0)) * dot - (oVec(1) - pVec(1)) * cross);
+ d_R = tempQAtan2;
+ }
+ const double d_x = temp0 * (oVec(1) * dot + oVec(0) * cross);
+ const double d_y = temp0 * (-oVec(0) * dot + oVec(1) * cross);
+ jxMat << d_X0, d_Y0, d_R, d_x, d_y, 0., 0., 0., 0., 0., 0., 1.;
+
+ covMat.block(0, 0, 3, 3) = circle.cov;
+ covMat(3, 3) = hits_ge.col(i)[0]; // x errors
+ covMat(4, 4) = hits_ge.col(i)[2]; // y errors
+ covMat(5, 5) = hits_ge.col(i)[5]; // z errors
+ covMat(3, 4) = covMat(4, 3) = hits_ge.col(i)[1]; // cov_xy
+ covMat(3, 5) = covMat(5, 3) = hits_ge.col(i)[3]; // cov_xz
+ covMat(4, 5) = covMat(5, 4) = hits_ge.col(i)[4]; // cov_yz
+ Matrix2d tmp = jxMat * covMat * jxMat.transpose();
+ cov_sz[i].noalias() = rot * tmp * rot.transpose();
+ }
+ // Math of d_{X0,Y0,R,x,y} all verified by hand
+ p2D.row(1) = hits.row(2);
+
+ // The following matrix will contain errors orthogonal to the rotated S
+ // component only, with the Multiple Scattering properly treated!!
+ MatrixNd cov_with_ms;
+ scatterCovLine(cov_sz, fast_fit, p2D.row(0), p2D.row(1), theta, bField, cov_with_ms);
+#ifdef RFIT_DEBUG
+ printIt(cov_sz, "line_fit - cov_sz:");
+ printIt(&cov_with_ms, "line_fit - cov_with_ms: ");
+#endif
+
+ // Rotate Points with the shape [2, n]
+ Matrix2xNd p2D_rot = rot * p2D;
+
+#ifdef RFIT_DEBUG
+ printf("Fast fit Tan(theta): %g\n", fast_fit(3));
+ printf("Rotation angle: %g\n", theta);
+ printIt(&rot, "Rotation Matrix:");
+ printIt(&p2D, "Original Hits(s,z):");
+ printIt(&p2D_rot, "Rotated hits(S3D, Z'):");
+ printIt(&rot, "Rotation Matrix:");
+#endif
+
+ // Build the A Matrix
+ Matrix2xNd aMat;
+ aMat << MatrixXd::Ones(1, n), p2D_rot.row(0); // rotated s values
+
+#ifdef RFIT_DEBUG
+ printIt(&aMat, "A Matrix:");
+#endif
+
+ // Build A^T V-1 A, where V-1 is the covariance of only the Y components.
+ MatrixNd vyInvMat;
+ math::cholesky::invert(cov_with_ms, vyInvMat);
+ // MatrixNd vyInvMat = cov_with_ms.inverse();
+ Eigen::Matrix covParamsMat = aMat * vyInvMat * aMat.transpose();
+ // Compute the Covariance Matrix of the fit parameters
+ math::cholesky::invert(covParamsMat, covParamsMat);
+
+ // Now Compute the Parameters in the form [2,1]
+ // The first component is q.
+ // The second component is m.
+ Eigen::Matrix sol = covParamsMat * aMat * vyInvMat * p2D_rot.row(1).transpose();
+
+#ifdef RFIT_DEBUG
+ printIt(&sol, "Rotated solutions:");
+#endif
+
+ // We need now to transfer back the results in the original s-z plane
+ const auto sinTheta = sin(theta);
+ const auto cosTheta = cos(theta);
+ auto common_factor = 1. / (sinTheta - sol(1, 0) * cosTheta);
+ Eigen::Matrix jMat;
+ jMat << 0., common_factor * common_factor, common_factor, sol(0, 0) * cosTheta * common_factor * common_factor;
+
+ double tempM = common_factor * (sol(1, 0) * sinTheta + cosTheta);
+ double tempQ = common_factor * sol(0, 0);
+ auto cov_mq = jMat * covParamsMat * jMat.transpose();
+
+ VectorNd res = p2D_rot.row(1).transpose() - aMat.transpose() * sol;
+ double chi2 = res.transpose() * vyInvMat * res;
+
+ LineFit line;
+ line.par << tempM, tempQ;
+ line.cov << cov_mq;
+ line.chi2 = chi2;
+
+#ifdef RFIT_DEBUG
+ printf("Common_factor: %g\n", common_factor);
+ printIt(&jMat, "Jacobian:");
+ printIt(&sol, "Rotated solutions:");
+ printIt(&covParamsMat, "Cov_params:");
+ printIt(&cov_mq, "Rotated Covariance Matrix:");
+ printIt(&(line.par), "Real Parameters:");
+ printIt(&(line.cov), "Real Covariance Matrix:");
+ printf("Chi2: %g\n", chi2);
+#endif
+
+ return line;
+ }
+
+ /*!
+ \brief Helix fit by three step:
+ -fast pre-fit (see Fast_fit() for further info); \n
+ -circle fit of hits projected in the transverse plane by Riemann-Chernov
+ algorithm (see Circle_fit() for further info); \n
+ -line fit of hits projected on cylinder surface by orthogonal distance
+ regression (see Line_fit for further info). \n
+ Points must be passed ordered (from inner to outer layer).
+ \param hits Matrix3xNd hits coordinates in this form: \n
+ |x0|x1|x2|...|xn| \n
+ |y0|y1|y2|...|yn| \n
+ |z0|z1|z2|...|zn|
+ \param hits_cov Matrix3Nd covariance matrix in this form (()->cov()): \n
+ |(x0,x0)|(x1,x0)|(x2,x0)|.|(y0,x0)|(y1,x0)|(y2,x0)|.|(z0,x0)|(z1,x0)|(z2,x0)| \n
+ |(x0,x1)|(x1,x1)|(x2,x1)|.|(y0,x1)|(y1,x1)|(y2,x1)|.|(z0,x1)|(z1,x1)|(z2,x1)| \n
+ |(x0,x2)|(x1,x2)|(x2,x2)|.|(y0,x2)|(y1,x2)|(y2,x2)|.|(z0,x2)|(z1,x2)|(z2,x2)| \n
+ . . . . . . . . . . . \n
+ |(x0,y0)|(x1,y0)|(x2,y0)|.|(y0,y0)|(y1,y0)|(y2,x0)|.|(z0,y0)|(z1,y0)|(z2,y0)| \n
+ |(x0,y1)|(x1,y1)|(x2,y1)|.|(y0,y1)|(y1,y1)|(y2,x1)|.|(z0,y1)|(z1,y1)|(z2,y1)| \n
+ |(x0,y2)|(x1,y2)|(x2,y2)|.|(y0,y2)|(y1,y2)|(y2,x2)|.|(z0,y2)|(z1,y2)|(z2,y2)| \n
+ . . . . . . . . . . . \n
+ |(x0,z0)|(x1,z0)|(x2,z0)|.|(y0,z0)|(y1,z0)|(y2,z0)|.|(z0,z0)|(z1,z0)|(z2,z0)| \n
+ |(x0,z1)|(x1,z1)|(x2,z1)|.|(y0,z1)|(y1,z1)|(y2,z1)|.|(z0,z1)|(z1,z1)|(z2,z1)| \n
+ |(x0,z2)|(x1,z2)|(x2,z2)|.|(y0,z2)|(y1,z2)|(y2,z2)|.|(z0,z2)|(z1,z2)|(z2,z2)|
+ \param bField magnetic field in the center of the detector in Gev/cm/c
+ unit, in order to perform pt calculation.
+ \param error flag for error computation.
+ \param scattering flag for multiple scattering treatment.
+ (see Circle_fit() documentation for further info).
+ \warning see Circle_fit(), Line_fit() and Fast_fit() warnings.
+ \bug see Circle_fit(), Line_fit() and Fast_fit() bugs.
+*/
+
+ template
+ inline HelixFit helixFit(const Matrix3xNd& hits,
+ const Eigen::Matrix& hits_ge,
+ const double bField,
+ const bool error) {
+ constexpr uint n = N;
+ VectorNd<4> rad = (hits.block(0, 0, 2, n).colwise().norm());
+
+ // Fast_fit gives back (X0, Y0, R, theta) w/o errors, using only 3 points.
+ Vector4d fast_fit;
+ fastFit(hits, fast_fit);
+ riemannFit::Matrix2Nd hits_cov = MatrixXd::Zero(2 * n, 2 * n);
+ riemannFit::loadCovariance2D(hits_ge, hits_cov);
+ CircleFit circle = circleFit(hits.block(0, 0, 2, n), hits_cov, fast_fit, rad, bField, error);
+ LineFit line = lineFit(hits, hits_ge, circle, fast_fit, bField, error);
+
+ par_uvrtopak(circle, bField, error);
+
+ HelixFit helix;
+ helix.par << circle.par, line.par;
+ if (error) {
+ helix.cov = MatrixXd::Zero(5, 5);
+ helix.cov.block(0, 0, 3, 3) = circle.cov;
+ helix.cov.block(3, 3, 2, 2) = line.cov;
+ }
+ helix.qCharge = circle.qCharge;
+ helix.chi2_circle = circle.chi2;
+ helix.chi2_line = line.chi2;
+
+ return helix;
+ }
+
+} // namespace riemannFit
+
+#endif // RecoPixelVertexing_PixelTrackFitting_interface_RiemannFit_h
diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/BuildFile.xml b/RecoPixelVertexing/PixelTrackFitting/plugins/BuildFile.xml
index be113d7a5a3dc..ecfbd99b667fc 100644
--- a/RecoPixelVertexing/PixelTrackFitting/plugins/BuildFile.xml
+++ b/RecoPixelVertexing/PixelTrackFitting/plugins/BuildFile.xml
@@ -1,3 +1,6 @@
+
+
+
diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelNtupletsFitterProducer.cc b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelNtupletsFitterProducer.cc
new file mode 100644
index 0000000000000..f49d2f01f48c6
--- /dev/null
+++ b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelNtupletsFitterProducer.cc
@@ -0,0 +1,44 @@
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/global/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "MagneticField/Engine/interface/MagneticField.h"
+#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
+#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelFitter.h"
+#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelNtupletsFitter.h"
+#include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
+
+class PixelNtupletsFitterProducer : public edm::global::EDProducer<> {
+public:
+ explicit PixelNtupletsFitterProducer(const edm::ParameterSet& iConfig)
+ : useRiemannFit_(iConfig.getParameter("useRiemannFit")), idealMagneticFieldToken_(esConsumes()) {
+ produces();
+ }
+ ~PixelNtupletsFitterProducer() override {}
+
+ static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+ edm::ParameterSetDescription desc;
+ desc.add("useRiemannFit", false)->setComment("true for Riemann, false for BrokenLine");
+ descriptions.add("pixelNtupletsFitterDefault", desc);
+ }
+
+private:
+ bool useRiemannFit_;
+ const edm::ESGetToken idealMagneticFieldToken_;
+ void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
+};
+
+void PixelNtupletsFitterProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
+ auto const& idealField = iSetup.getData(idealMagneticFieldToken_);
+ float bField = 1 / PixelRecoUtilities::fieldInInvGev(iSetup);
+ auto impl = std::make_unique(bField, &idealField, useRiemannFit_);
+ auto prod = std::make_unique(std::move(impl));
+ iEvent.put(std::move(prod));
+}
+
+DEFINE_FWK_MODULE(PixelNtupletsFitterProducer);
diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.cc b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.cc
index bd390f5f65352..91c3a44cc8643 100644
--- a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.cc
+++ b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.cc
@@ -1,23 +1,22 @@
-#include "PixelTrackProducer.h"
+#include
+#include "DataFormats/Common/interface/OrphanHandle.h"
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackExtra.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
+#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
+#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
-#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
-#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
-
-#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
-#include "DataFormats/TrackReco/interface/Track.h"
-#include "DataFormats/TrackReco/interface/TrackFwd.h"
-#include "DataFormats/TrackReco/interface/TrackExtra.h"
-#include "DataFormats/Common/interface/OrphanHandle.h"
-
-#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
-#include
+#include "PixelTrackProducer.h"
+#include "storeTracks.h"
using namespace pixeltrackfitting;
using edm::ParameterSet;
@@ -45,62 +44,9 @@ void PixelTrackProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
TracksWithTTRHs tracks;
theReconstruction.run(tracks, ev, es);
-
edm::ESHandle httopo;
es.get().get(httopo);
// store tracks
- store(ev, tracks, *httopo);
-}
-
-void PixelTrackProducer::store(edm::Event& ev, const TracksWithTTRHs& tracksWithHits, const TrackerTopology& ttopo) {
- auto tracks = std::make_unique();
- auto recHits = std::make_unique();
- auto trackExtras = std::make_unique();
-
- int cc = 0, nTracks = tracksWithHits.size();
-
- for (int i = 0; i < nTracks; i++) {
- reco::Track* track = tracksWithHits.at(i).first;
- const SeedingHitSet& hits = tracksWithHits.at(i).second;
-
- for (unsigned int k = 0; k < hits.size(); k++) {
- TrackingRecHit* hit = hits[k]->hit()->clone();
-
- track->appendHitPattern(*hit, ttopo);
- recHits->push_back(hit);
- }
- tracks->push_back(*track);
- delete track;
- }
-
- LogDebug("TrackProducer") << "put the collection of TrackingRecHit in the event"
- << "\n";
- edm::OrphanHandle ohRH = ev.put(std::move(recHits));
-
- edm::RefProd hitCollProd(ohRH);
- for (int k = 0; k < nTracks; k++) {
- reco::TrackExtra theTrackExtra{};
-
- //fill the TrackExtra with TrackingRecHitRef
- unsigned int nHits = tracks->at(k).numberOfValidHits();
- theTrackExtra.setHits(hitCollProd, cc, nHits);
- cc += nHits;
- AlgebraicVector5 v = AlgebraicVector5(0, 0, 0, 0, 0);
- reco::TrackExtra::TrajParams trajParams(nHits, LocalTrajectoryParameters(v, 1.));
- reco::TrackExtra::Chi2sFive chi2s(nHits, 0);
- theTrackExtra.setTrajParams(std::move(trajParams), std::move(chi2s));
- trackExtras->push_back(theTrackExtra);
- }
-
- LogDebug("TrackProducer") << "put the collection of TrackExtra in the event"
- << "\n";
- edm::OrphanHandle ohTE = ev.put(std::move(trackExtras));
-
- for (int k = 0; k < nTracks; k++) {
- const reco::TrackExtraRef theTrackExtraRef(ohTE, k);
- (tracks->at(k)).setExtra(theTrackExtraRef);
- }
-
- ev.put(std::move(tracks));
+ storeTracks(ev, tracks, *httopo);
}
diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.h b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.h
index d756a9cf963f5..c38fd44c0d7f5 100644
--- a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.h
+++ b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.h
@@ -1,8 +1,7 @@
-#ifndef PixelTrackProducer_H
-#define PixelTrackProducer_H
+#ifndef RecoPixelVertexing_PixelTrackFitting_plugins_PixelTrackProducer_h
+#define RecoPixelVertexing_PixelTrackFitting_plugins_PixelTrackProducer_h
#include "FWCore/Framework/interface/stream/EDProducer.h"
-#include "RecoPixelVertexing/PixelTrackFitting/interface/TracksWithHits.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackReconstruction.h"
namespace edm {
@@ -24,7 +23,7 @@ class PixelTrackProducer : public edm::stream::EDProducer<> {
void produce(edm::Event& ev, const edm::EventSetup& es) override;
private:
- void store(edm::Event& ev, const pixeltrackfitting::TracksWithTTRHs& selectedTracks, const TrackerTopology& ttopo);
PixelTrackReconstruction theReconstruction;
};
-#endif
+
+#endif // RecoPixelVertexing_PixelTrackFitting_plugins_PixelTrackProducer_h
diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc
new file mode 100644
index 0000000000000..94c490e948575
--- /dev/null
+++ b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc
@@ -0,0 +1,205 @@
+#include "DataFormats/BeamSpot/interface/BeamSpot.h"
+#include "DataFormats/Common/interface/OrphanHandle.h"
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackExtra.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
+#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
+#include "DataFormats/GeometrySurface/interface/Plane.h"
+#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/global/EDProducer.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/PluginManager/interface/ModuleDef.h"
+#include "FWCore/Utilities/interface/EDGetToken.h"
+#include "Geometry/Records/interface/TrackerTopologyRcd.h"
+#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
+
+#include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
+#include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
+#include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
+#include "RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h"
+
+#include "CUDADataFormats/Common/interface/HostProduct.h"
+#include "CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h"
+#include "CUDADataFormats/SiPixelCluster/interface/gpuClusteringConstants.h"
+
+#include "storeTracks.h"
+#include "CUDADataFormats/Common/interface/HostProduct.h"
+
+/**
+ * This class creates "leagcy" reco::Track
+ * objects from the output of SoA CA.
+ */
+class PixelTrackProducerFromSoA : public edm::global::EDProducer<> {
+public:
+ using IndToEdm = std::vector;
+
+ explicit PixelTrackProducerFromSoA(const edm::ParameterSet &iConfig);
+ ~PixelTrackProducerFromSoA() override = default;
+
+ static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
+
+ // using HitModuleStart = std::array;
+ using HMSstorage = HostProduct;
+
+private:
+ void produce(edm::StreamID streamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const override;
+
+ // Event Data tokens
+ const edm::EDGetTokenT tBeamSpot_;
+ const edm::EDGetTokenT tokenTrack_;
+ const edm::EDGetTokenT cpuHits_;
+ const edm::EDGetTokenT hmsToken_;
+ // Event Setup tokens
+ const edm::ESGetToken idealMagneticFieldToken_;
+ const edm::ESGetToken ttTopoToken_;
+
+ int32_t const minNumberOfHits_;
+};
+
+PixelTrackProducerFromSoA::PixelTrackProducerFromSoA(const edm::ParameterSet &iConfig)
+ : tBeamSpot_(consumes(iConfig.getParameter("beamSpot"))),
+ tokenTrack_(consumes(iConfig.getParameter("trackSrc"))),
+ cpuHits_(consumes