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/CUDADataFormats/Vertex/BuildFile.xml b/CUDADataFormats/Vertex/BuildFile.xml
new file mode 100644
index 0000000000000..e3f9a0910bbd8
--- /dev/null
+++ b/CUDADataFormats/Vertex/BuildFile.xml
@@ -0,0 +1,9 @@
+
+
+
+
+
+
+
+
+
diff --git a/CUDADataFormats/Vertex/interface/ZVertexHeterogeneous.h b/CUDADataFormats/Vertex/interface/ZVertexHeterogeneous.h
new file mode 100644
index 0000000000000..aacfddc6fe7e2
--- /dev/null
+++ b/CUDADataFormats/Vertex/interface/ZVertexHeterogeneous.h
@@ -0,0 +1,14 @@
+#ifndef CUDADataFormatsVertexZVertexHeterogeneous_H
+#define CUDADataFormatsVertexZVertexHeterogeneous_H
+
+#include "CUDADataFormats/Vertex/interface/ZVertexSoA.h"
+#include "CUDADataFormats/Common/interface/HeterogeneousSoA.h"
+#include "CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h"
+
+using ZVertexHeterogeneous = HeterogeneousSoA;
+#ifndef __CUDACC__
+#include "CUDADataFormats/Common/interface/Product.h"
+using ZVertexCUDAProduct = cms::cuda::Product;
+#endif
+
+#endif
diff --git a/CUDADataFormats/Vertex/interface/ZVertexSoA.h b/CUDADataFormats/Vertex/interface/ZVertexSoA.h
new file mode 100644
index 0000000000000..5f0699d5831ec
--- /dev/null
+++ b/CUDADataFormats/Vertex/interface/ZVertexSoA.h
@@ -0,0 +1,26 @@
+#ifndef CUDADataFormatsVertexZVertexSoA_H
+#define CUDADataFormatsVertexZVertexSoA_H
+
+#include
+#include "HeterogeneousCore/CUDAUtilities/interface/cudaCompat.h"
+
+// SOA for vertices
+// These vertices are clusterized and fitted only along the beam line (z)
+// to obtain their global coordinate the beam spot position shall be added (eventually correcting for the beam angle as well)
+struct ZVertexSoA {
+ static constexpr uint32_t MAXTRACKS = 32 * 1024;
+ static constexpr uint32_t MAXVTX = 1024;
+
+ int16_t idv[MAXTRACKS]; // vertex index for each associated (original) track (-1 == not associate)
+ float zv[MAXVTX]; // output z-posistion of found vertices
+ float wv[MAXVTX]; // output weight (1/error^2) on the above
+ float chi2[MAXVTX]; // vertices chi2
+ float ptv2[MAXVTX]; // vertices pt^2
+ int32_t ndof[MAXTRACKS]; // vertices number of dof (reused as workspace for the number of nearest neighbours FIXME)
+ uint16_t sortInd[MAXVTX]; // sorted index (by pt2) ascending
+ uint32_t nvFinal; // the number of vertices
+
+ __host__ __device__ void init() { nvFinal = 0; }
+};
+
+#endif // CUDADataFormatsVertexZVertexSoA.H
diff --git a/CUDADataFormats/Vertex/src/classes.h b/CUDADataFormats/Vertex/src/classes.h
new file mode 100644
index 0000000000000..e7fea871f7d39
--- /dev/null
+++ b/CUDADataFormats/Vertex/src/classes.h
@@ -0,0 +1,8 @@
+#ifndef CUDADataFormats__src_classes_h
+#define CUDADataFormats__src_classes_h
+
+#include "CUDADataFormats/Vertex/interface/ZVertexHeterogeneous.h"
+#include "CUDADataFormats/Common/interface/Product.h"
+#include "DataFormats/Common/interface/Wrapper.h"
+
+#endif
diff --git a/CUDADataFormats/Vertex/src/classes_def.xml b/CUDADataFormats/Vertex/src/classes_def.xml
new file mode 100644
index 0000000000000..ea633080af9af
--- /dev/null
+++ b/CUDADataFormats/Vertex/src/classes_def.xml
@@ -0,0 +1,6 @@
+
+
+
+
+
+
diff --git a/Configuration/PyReleaseValidation/python/relval_2017.py b/Configuration/PyReleaseValidation/python/relval_2017.py
index 193ab79bcd384..87130136c154c 100644
--- a/Configuration/PyReleaseValidation/python/relval_2017.py
+++ b/Configuration/PyReleaseValidation/python/relval_2017.py
@@ -5,7 +5,7 @@
# here only define the workflows as a combination of the steps defined above:
workflows = Matrix()
-# each workflow defines a name and a list of steps to be done.
+# each workflow defines a name and a list of steps to be done.
# if no explicit name/label given for the workflow (first arg),
# the name of step1 will be used
@@ -24,16 +24,16 @@
# (HE collapse: TTbar, TTbar PU, TTbar design)
# (ParkingBPH: TTbar)
# (TTbar PU with JME NanoAOD)
-# (Patatrack pixel-only: ZMM - on CPU)
-# (Patatrack pixel-only: TTbar - on CPU)
+# (Patatrack pixel-only: ZMM - on CPU: quadruplets, triplets)
+# (Patatrack pixel-only: TTbar - on CPU: quadruplets, triplets)
# (Patatrack ECAL-only: TTbar - on CPU)
# (Patatrack HCAL-only: TTbar - on CPU)
# 2021 (DD4HEP: TTbar, ZMM)
# (ele guns 10, 35, 1000; pho guns 10, 35; mu guns 1, 10, 100, 1000, QCD 3TeV, QCD Flat)
# (ZMM, TTbar, ZEE, MinBias, TTbar PU, TTbar PU premix, ZEE PU, TTbar design)
# (TTbar trackingOnly, pixelTrackingOnly, trackingMkFit, trackdnn)
-# (Patatrack pixel-only: ZMM - on CPU)
-# (Patatrack pixel-only: TTbar - on CPU)
+# (Patatrack pixel-only: ZMM - on CPU: quadruplets, triplets)
+# (Patatrack pixel-only: TTbar - on CPU: quadruplets, triplets)
# (Patatrack ECAL-only: TTbar - on CPU)
# (Patatrack HCAL-only: TTbar - on CPU)
# (TTbar 0T, TTbar PU 0T)
@@ -51,16 +51,16 @@
10824.6,11024.6,11224.6,
10824.8,
11024.15,
- 10842.501,
- 10824.501,
+ 10842.501,10842.505,
+ 10824.501,10824.505,
10824.511,
10824.521,
11634.911, 11650.911,
11601.0,11602.0,11603.0,11604.0,11605.0,11606.0,11607.0,11608.0,11609.0,11630.0,11643.0,
11650.0,11634.0,11646.0,11640.0,11834.0,11834.99,11846.0,12024.0,
11634.1,11634.5,11634.7,11634.91,
- 11650.501,
- 11634.501,
+ 11650.501,11650.505,
+ 11634.501,11634.505,
11634.511,
11634.521,
11634.24,11834.24,
diff --git a/Configuration/PyReleaseValidation/python/relval_gpu.py b/Configuration/PyReleaseValidation/python/relval_gpu.py
index 4e49467a0e2e8..43353279ea4ad 100644
--- a/Configuration/PyReleaseValidation/python/relval_gpu.py
+++ b/Configuration/PyReleaseValidation/python/relval_gpu.py
@@ -5,7 +5,7 @@
# here only define the workflows as a combination of the steps defined above:
workflows = Matrix()
-# each workflow defines a name and a list of steps to be done.
+# each workflow defines a name and a list of steps to be done.
# if no explicit name/label given for the workflow (first arg),
# the name of step1 will be used
@@ -14,21 +14,29 @@
#just define all of them
#WFs to run in IB:
-# mc 2018 (Patatrack pixel-only: ZMM - on GPU, both CPU and GPU, auto)
-# (Patatrack pixel-only: TTbar - on GPU, both CPU and GPU, auto)
+# mc 2018 (Patatrack pixel-only quadruplets: ZMM - on GPU, both CPU and GPU, auto)
+# (Patatrack pixel-only triplets: ZMM - on GPU, both CPU and GPU, auto)
+# (Patatrack pixel-only quadruplets: TTbar - on GPU, both CPU and GPU, auto)
+# (Patatrack pixel-only triplets: TTbar - on GPU, both CPU and GPU, auto)
# (Patatrack ECAL-only: TTbar - on GPU, both CPU and GPU, auto)
# (Patatrack HCAL-only: TTbar - on GPU, both CPU and GPU, auto)
-# mc 2021 (Patatrack pixel-only: ZMM - on GPU, both CPU and GPU, auto)
-# (Patatrack pixel-only: TTbar - on GPU, both CPU and GPU, auto)
+# mc 2021 (Patatrack pixel-only quadruplets: ZMM - on GPU, both CPU and GPU, auto)
+# (Patatrack pixel-only triplets: ZMM - on GPU, both CPU and GPU, auto)
+# (Patatrack pixel-only quadruplets: TTbar - on GPU, both CPU and GPU, auto)
+# (Patatrack pixel-only triplets: TTbar - on GPU, both CPU and GPU, auto)
# (Patatrack ECAL-only: TTbar - on GPU, both CPU and GPU, auto)
# (Patatrack HCAL-only: TTbar - on GPU, both CPU and GPU, auto)
numWFIB = [
10842.502, # 10842.503,10842.504,
+ 10842.506, # 10842.507,10842.508,
10824.502, # 10824.503,10824.504,
+ 10824.506, # 10824.507,10824.508,
10824.512, # 10824.513,10824.514,
10824.522, # 10824.523,10824.524,
11650.502, # 11650.503,11650.504,
+ 11650.506, # 11650.507,11650.508,
11634.502, # 11634.503,11634.504,
+ 11634.506, # 11634.507,11634.508,
11634.512, # 11634.513,11634.514,
11634.522, # 11634.523,11634.524
]
diff --git a/Configuration/PyReleaseValidation/python/relval_steps.py b/Configuration/PyReleaseValidation/python/relval_steps.py
index 3556eda005ada..0fc667283361b 100644
--- a/Configuration/PyReleaseValidation/python/relval_steps.py
+++ b/Configuration/PyReleaseValidation/python/relval_steps.py
@@ -2186,8 +2186,11 @@ def gen2021HiMix(fragment,howMuch):
'--era' :'Run2_2016'
}
-step3_pixelNtupleFit = {
- '--procModifiers': 'pixelNtupleFit',
+step3_pixel_ntuplet_cpu = {
+ '--customise': 'RecoPixelVertexing/Configuration/customizePixelTracksSoAonCPU.customizePixelTracksSoAonCPU'
+}
+step3_pixel_triplets = {
+ '--customise': 'RecoPixelVertexing/Configuration/customizePixelTracksSoAonCPU.customizePixelTracksForTriplets'
}
step3_gpu = {
'--procModifiers': 'gpu',
@@ -2320,8 +2323,11 @@ def gen2021HiMix(fragment,howMuch):
steps['RECODR2_2018reHLT_Prompt']=merge([{'--conditions':'auto:run2_data'},steps['RECODR2_2018reHLT']])
steps['RECODR2_2018reHLT_ZBPrompt']=merge([{'--conditions':'auto:run2_data','-s':'RAW2DIGI,L1Reco,RECO,EI,PAT,ALCA:SiStripCalZeroBias+SiStripCalMinBias+TkAlMinBias+EcalESAlign,DQM:@rerecoZeroBias+@ExtraHLT+@miniAODDQM'},steps['RECODR2_2018reHLT']])
steps['RECODR2_2018reHLT_Prompt_pixelTrackingOnly']=merge([{'-s': 'RAW2DIGI:RawToDigi_pixelOnly,RECO:reconstruction_pixelTrackingOnly,DQM:@pixelTrackingOnlyDQM'},steps['RECODR2_2018reHLT_Prompt']])
-steps['RECODR2_2018reHLT_Patatrack_PixelOnlyCPU']=merge([step3_pixelNtupleFit, steps['RECODR2_2018reHLT_Prompt_pixelTrackingOnly']])
+steps['RECODR2_2018reHLT_Patatrack_PixelOnlyCPU']=merge([step3_pixel_ntuplet_cpu, steps['RECODR2_2018reHLT_Prompt_pixelTrackingOnly']])
steps['RECODR2_2018reHLT_Patatrack_PixelOnlyGPU']=merge([step3_gpu, steps['RECODR2_2018reHLT_Prompt_pixelTrackingOnly']])
+steps['RECODR2_2018reHLT_Patatrack_PixelOnlyTripletsCPU']=merge([step3_pixel_ntuplet_cpu, step3_pixel_triplets, steps['RECODR2_2018reHLT_Prompt_pixelTrackingOnly']])
+steps['RECODR2_2018reHLT_Patatrack_PixelOnlyTripletsGPU']=merge([step3_gpu, step3_pixel_triplets, steps['RECODR2_2018reHLT_Prompt_pixelTrackingOnly']])
+
steps['RECODR2_2018reHLT_ECALOnlyCPU']=merge([{'-s': 'RAW2DIGI:RawToDigi_ecalOnly,RECO:reconstruction_ecalOnly,DQM:@ecalOnly'},steps['RECODR2_2018reHLT_Prompt']])
steps['RECODR2_2018reHLT_ECALOnlyGPU']=merge([step3_gpu, steps['RECODR2_2018reHLT_ECALOnlyCPU']])
steps['RECODR2_2018reHLT_HCALOnlyCPU']=merge([{'-s': 'RAW2DIGI:RawToDigi_hcalOnly,RECO:reconstruction_hcalOnly,DQM:@hcalOnly+@hcal2Only'},steps['RECODR2_2018reHLT_Prompt']])
diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py
index 9ff0ae9f22e0d..513462593300e 100644
--- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py
+++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py
@@ -453,7 +453,26 @@ def condition_(self, fragment, stepList, key, hasHarvest):
'-s': 'RAW2DIGI:RawToDigi_pixelOnly,RECO:reconstruction_pixelTrackingOnly,VALIDATION:@pixelTrackingOnlyValidation,DQM:@pixelTrackingOnlyDQM',
'--datatier': 'GEN-SIM-RECO,DQMIO',
'--eventcontent': 'RECOSIM,DQM',
- '--procModifiers': 'pixelNtupleFit'
+ '--customise' : 'RecoPixelVertexing/Configuration/customizePixelTracksSoAonCPU.customizePixelTracksSoAonCPU'
+}
+
+upgradeWFs['PatatrackPixelOnlyTripletsCPU'] = UpgradeWorkflowPatatrack_PixelOnlyCPU(
+ steps = [
+ 'Reco',
+ 'HARVEST',
+ 'RecoGlobal',
+ 'HARVESTGlobal',
+ ],
+ PU = [],
+ suffix = 'Patatrack_PixelOnlyTripletsCPU',
+ offset = 0.505,
+)
+
+upgradeWFs['PatatrackPixelOnlyTripletsCPU'].step3 = {
+ '-s': 'RAW2DIGI:RawToDigi_pixelOnly,RECO:reconstruction_pixelTrackingOnly,VALIDATION:@pixelTrackingOnlyValidation,DQM:@pixelTrackingOnlyDQM',
+ '--datatier': 'GEN-SIM-RECO,DQMIO',
+ '--eventcontent': 'RECOSIM,DQM',
+ '--customise' : 'RecoPixelVertexing/Configuration/customizePixelTracksSoAonCPU.customizePixelTracksSoAonCPU,RecoPixelVertexing/Configuration/customizePixelTracksSoAonCPU.customizePixelTracksForTriplets'
}
class UpgradeWorkflowPatatrack_PixelOnlyGPU(UpgradeWorkflowPatatrack):
@@ -487,6 +506,26 @@ def condition_(self, fragment, stepList, key, hasHarvest):
'--procModifiers': 'gpu'
}
+upgradeWFs['PatatrackPixelOnlyTripletsGPU'] = UpgradeWorkflowPatatrack_PixelOnlyGPU(
+ steps = [
+ 'Reco',
+ 'HARVEST',
+ 'RecoGlobal',
+ 'HARVESTGlobal',
+ ],
+ PU = [],
+ suffix = 'Patatrack_PixelOnlyTripletsGPU',
+ offset = 0.506,
+)
+
+upgradeWFs['PatatrackPixelOnlyTripletsGPU'].step3 = {
+ '-s': 'RAW2DIGI:RawToDigi_pixelOnly,RECO:reconstruction_pixelTrackingOnly,VALIDATION:@pixelTrackingOnlyValidation,DQM:@pixelTrackingOnlyDQM',
+ '--datatier': 'GEN-SIM-RECO,DQMIO',
+ '--eventcontent': 'RECOSIM,DQM',
+ '--procModifiers': 'gpu',
+ '--customise': 'RecoPixelVertexing/Configuration/customizePixelTracksSoAonCPU.customizePixelTracksForTriplets'
+}
+
class UpgradeWorkflowPatatrack_ECALOnlyCPU(UpgradeWorkflowPatatrack):
def setup_(self, step, stepName, stepDict, k, properties):
if 'Reco' in step:
diff --git a/DQM/Integration/python/clients/beam_dqm_sourceclient-live_cfg.py b/DQM/Integration/python/clients/beam_dqm_sourceclient-live_cfg.py
index efdde1512fcf7..4846de0887fde 100644
--- a/DQM/Integration/python/clients/beam_dqm_sourceclient-live_cfg.py
+++ b/DQM/Integration/python/clients/beam_dqm_sourceclient-live_cfg.py
@@ -317,8 +317,7 @@
process.pixelTracksTrackingRegions.RegionPSet.originXPos = 0.08
process.pixelTracksTrackingRegions.RegionPSet.originYPos = -0.03
process.pixelTracksTrackingRegions.RegionPSet.originZPos = 0.
-
-process.pixelVertices.TkFilterParameters.minPt = process.pixelTracksTrackingRegions.RegionPSet.ptMin
+process.pixelVertices.PtMin = process.pixelTracksTrackingRegions.RegionPSet.ptMin
process.tracking_FirstStep = cms.Sequence(
process.siPixelDigis
diff --git a/DQM/Integration/python/clients/beampixel_dqm_sourceclient-live_cfg.py b/DQM/Integration/python/clients/beampixel_dqm_sourceclient-live_cfg.py
index f909104a39834..a3eac2069e6ed 100644
--- a/DQM/Integration/python/clients/beampixel_dqm_sourceclient-live_cfg.py
+++ b/DQM/Integration/python/clients/beampixel_dqm_sourceclient-live_cfg.py
@@ -90,12 +90,12 @@
process.siPixelClusterShapeCachePreSplitting = siPixelClusterShapeCache.clone(src = 'siPixelClustersPreSplitting')
process.load("RecoLocalTracker.SiPixelRecHits.PixelCPEGeneric_cfi")
process.load("RecoPixelVertexing.Configuration.RecoPixelVertexing_cff")
-process.pixelVertices.TkFilterParameters.minPt = process.pixelTracksTrackingRegions.RegionPSet.ptMin
process.pixelTracksTrackingRegions.RegionPSet.originRadius = cms.double(0.4)
process.pixelTracksTrackingRegions.RegionPSet.originHalfLength = cms.double(15.)
process.pixelTracksTrackingRegions.RegionPSet.originXPos = cms.double(0.08)
process.pixelTracksTrackingRegions.RegionPSet.originYPos = cms.double(-0.03)
process.pixelTracksTrackingRegions.RegionPSet.originZPos = cms.double(0.)
+process.pixelVertices.PtMin = process.pixelTracksTrackingRegions.RegionPSet.ptMin
#----------------------------
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/TrackingMonitorClient/python/pixelVertexResolutionClient_cfi.py b/DQM/TrackingMonitorClient/python/pixelVertexResolutionClient_cfi.py
new file mode 100644
index 0000000000000..2558e88d26012
--- /dev/null
+++ b/DQM/TrackingMonitorClient/python/pixelVertexResolutionClient_cfi.py
@@ -0,0 +1,7 @@
+import FWCore.ParameterSet.Config as cms
+
+from DQM.TrackingMonitorClient.primaryVertexResolutionClient_cfi import primaryVertexResolutionClient as _primaryVertexResolutionClient
+
+pixelVertexResolutionClient = _primaryVertexResolutionClient.clone(
+ subDirs = ["OfflinePixelPV/Resolution/*"]
+)
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/DQMOffline/Configuration/python/DQMOffline_SecondStep_cff.py b/DQMOffline/Configuration/python/DQMOffline_SecondStep_cff.py
index 368b328632fd8..29bf311c474d4 100644
--- a/DQMOffline/Configuration/python/DQMOffline_SecondStep_cff.py
+++ b/DQMOffline/Configuration/python/DQMOffline_SecondStep_cff.py
@@ -122,6 +122,7 @@
from DQM.CTPPS.ctppsDQM_cff import *
from Validation.RecoTau.DQMSequences_cfi import *
from DQM.TrackingMonitorClient.pixelTrackingEffFromHitPattern_cff import *
+from DQM.TrackingMonitorClient.pixelVertexResolutionClient_cfi import *
DQMHarvestTrackerStrip = cms.Sequence ( SiStripOfflineDQMClient )
@@ -179,7 +180,8 @@
DQMHarvestTrackingZeroBias = cms.Sequence( TrackingOfflineDQMClientZeroBias *
dqmFastTimerServiceClient )
-DQMHarvestPixelTracking = cms.Sequence( pixelTrackingEffFromHitPattern )
+DQMHarvestPixelTracking = cms.Sequence( pixelTrackingEffFromHitPattern *
+ pixelVertexResolutionClient )
DQMHarvestOuterTracker = cms.Sequence(
OuterTrackerClient *
diff --git a/DQMOffline/Configuration/python/DQMOffline_cff.py b/DQMOffline/Configuration/python/DQMOffline_cff.py
index 2001c22352a48..ac28700d4eaf4 100644
--- a/DQMOffline/Configuration/python/DQMOffline_cff.py
+++ b/DQMOffline/Configuration/python/DQMOffline_cff.py
@@ -157,10 +157,12 @@
#DQMOfflineCommon
from DQM.TrackingMonitorSource.pixelTracksMonitoring_cff import *
+from DQMOffline.RecoB.PixelVertexMonitor_cff import *
from DQM.SiOuterTracker.OuterTrackerSourceConfig_cff import *
from Validation.RecoTau.DQMSequences_cfi import *
-DQMOfflinePixelTracking = cms.Sequence( pixelTracksMonitoring )
+DQMOfflinePixelTracking = cms.Sequence( pixelTracksMonitoring *
+ pixelPVMonitor )
DQMOuterTracker = cms.Sequence( DQMOfflineDCS *
OuterTrackerSource *
diff --git a/DQMOffline/RecoB/python/PixelVertexMonitor_cff.py b/DQMOffline/RecoB/python/PixelVertexMonitor_cff.py
new file mode 100644
index 0000000000000..9e293f4478bd6
--- /dev/null
+++ b/DQMOffline/RecoB/python/PixelVertexMonitor_cff.py
@@ -0,0 +1,8 @@
+import FWCore.ParameterSet.Config as cms
+
+from DQMOffline.RecoB.PrimaryVertexMonitor_cff import pvMonitor as _pvMonitor
+pixelPVMonitor = _pvMonitor.clone(
+ TopFolderName = "OfflinePixelPV",
+ vertexLabel = "pixelVertices",
+ ndof = cms.int32( 1 )
+)
diff --git a/FastSimulation/Tracking/python/SeedingMigration.py b/FastSimulation/Tracking/python/SeedingMigration.py
index 751670daa50c8..3a982eba55e36 100644
--- a/FastSimulation/Tracking/python/SeedingMigration.py
+++ b/FastSimulation/Tracking/python/SeedingMigration.py
@@ -13,8 +13,9 @@ def _hitSetProducerToFactoryPSet(producer):
"PixelTripletLargeTipEDProducer": "PixelTripletLargeTipGenerator",
"MultiHitFromChi2EDProducer": "MultiHitGeneratorFromChi2",
"CAHitTripletEDProducer": "CAHitTripletGenerator",
- "CAHitQuadrupletEDProducer": "CAHitQuadrupletGenerator",
- }
+ "CAHitQuadrupletEDProducer": "CAHitQuadrupletGenerator",
+ "CAHitNtupletHeterogeneousEDProducer": "CAHitQuadrupletGenerator",
+ }
ret = cms.PSet()
_copy(producer, ret)
ret.ComponentName = cms.string(_map[producer._TypedParameterizable__type]);
diff --git a/HLTrigger/Configuration/python/customizeHLTforPatatrack.py b/HLTrigger/Configuration/python/customizeHLTforPatatrack.py
new file mode 100644
index 0000000000000..5164188c94997
--- /dev/null
+++ b/HLTrigger/Configuration/python/customizeHLTforPatatrack.py
@@ -0,0 +1,802 @@
+import copy
+import FWCore.ParameterSet.Config as cms
+from HeterogeneousCore.CUDACore.SwitchProducerCUDA import SwitchProducerCUDA
+from HLTrigger.Configuration.common import *
+from Configuration.Eras.Modifier_run3_common_cff import run3_common
+
+
+# force the SwitchProducerCUDA choice to pick a specific backend: True for offloading to a gpu, False for running on cpu
+def forceGpuOffload(status = True):
+ import HeterogeneousCore.CUDACore.SwitchProducerCUDA
+ HeterogeneousCore.CUDACore.SwitchProducerCUDA._cuda_enabled_cached = bool(status)
+
+
+# reset the SwitchProducerCUDA choice to pick a backend depending on the availability of a supported gpu
+def resetGpuOffload():
+ import HeterogeneousCore.CUDACore.SwitchProducerCUDA
+ HeterogeneousCore.CUDACore.SwitchProducerCUDA._cuda_enabled_cached = None
+ HeterogeneousCore.CUDACore.SwitchProducerCUDA._switch_cuda()
+
+
+# customisation for running the Patatrack reconstruction, common parts
+def customiseCommon(process):
+
+ # Services
+
+ process.load("HeterogeneousCore.CUDAServices.CUDAService_cfi")
+
+ process.load("HeterogeneousCore.CUDAServices.NVProfilerService_cfi")
+
+
+ # Paths and EndPaths
+
+ # the hltGetConditions module would force gpu-specific ESProducers to run even if no supported gpu is present
+ if 'hltGetConditions' in process.__dict__:
+ del process.hltGetConditions
+
+ # produce a boolean to track if the events ar being processed on gpu (true) or cpu (false)
+ process.statusOnGPU = SwitchProducerCUDA(
+ cpu = cms.EDProducer("BooleanProducer", value = cms.bool(False)),
+ cuda = cms.EDProducer("BooleanProducer", value = cms.bool(True))
+ )
+
+ process.statusOnGPUFilter = cms.EDFilter("BooleanFilter",
+ src = cms.InputTag("statusOnGPU")
+ )
+
+ if 'Status_OnGPU' in process.__dict__:
+ replace_with(process.Status_OnGPU, cms.Path(process.statusOnGPU + process.statusOnGPUFilter))
+ else:
+ process.Status_OnGPU = cms.Path(process.statusOnGPU + process.statusOnGPUFilter)
+ if 'HLTSchedule' in process.__dict__:
+ process.HLTSchedule.append(process.Status_OnGPU)
+ if process.schedule is not None:
+ process.schedule.append(process.Status_OnGPU)
+
+
+ # make the ScoutingCaloMuonOutput endpath compatible with using Tasks in the Scouting paths
+ if 'hltOutputScoutingCaloMuon' in process.__dict__ and not 'hltPreScoutingCaloMuonOutputSmart' in process.__dict__:
+ process.hltPreScoutingCaloMuonOutputSmart = cms.EDFilter( "TriggerResultsFilter",
+ l1tIgnoreMaskAndPrescale = cms.bool( False ),
+ l1tResults = cms.InputTag( "" ),
+ hltResults = cms.InputTag( 'TriggerResults','','@currentProcess' ),
+ triggerConditions = process.hltOutputScoutingCaloMuon.SelectEvents.SelectEvents,
+ throw = cms.bool( True )
+ )
+ insert_modules_after(process, process.hltPreScoutingCaloMuonOutput, process.hltPreScoutingCaloMuonOutputSmart)
+
+ # make the ScoutingPFOutput endpath compatible with using Tasks in the Scouting paths
+ if 'hltOutputScoutingPF' in process.__dict__ and not 'hltPreScoutingPFOutputSmart' in process.__dict__:
+ process.hltPreScoutingPFOutputSmart = cms.EDFilter( "TriggerResultsFilter",
+ l1tIgnoreMaskAndPrescale = cms.bool( False ),
+ l1tResults = cms.InputTag( "" ),
+ hltResults = cms.InputTag( 'TriggerResults','','@currentProcess' ),
+ triggerConditions = process.hltOutputScoutingPF.SelectEvents.SelectEvents,
+ throw = cms.bool( True )
+ )
+ insert_modules_after(process, process.hltPreScoutingPFOutput, process.hltPreScoutingPFOutputSmart)
+
+
+ # done
+ return process
+
+
+# customisation for running the "Patatrack" pixel local reconstruction
+def customisePixelLocalReconstruction(process):
+
+ if not 'HLTDoLocalPixelSequence' in process.__dict__:
+ return process
+
+
+ # FIXME replace the Sequences with empty ones to avoid exanding them during the (re)definition of Modules and EDAliases
+
+ process.HLTDoLocalPixelSequence = cms.Sequence()
+
+
+ # Event Setup
+
+ process.load("CalibTracker.SiPixelESProducers.siPixelGainCalibrationForHLTGPU_cfi") # this should be used only on GPUs, will crash otherwise
+ process.load("CalibTracker.SiPixelESProducers.siPixelROCsStatusAndMappingWrapperESProducer_cfi") # this should be used only on GPUs, will crash otherwise
+ process.load("RecoLocalTracker.SiPixelRecHits.PixelCPEFastESProducer_cfi")
+
+
+ # Modules and EDAliases
+
+ # referenced in HLTDoLocalPixelTask
+
+ # transfer the beamspot to the gpu
+ from RecoVertex.BeamSpotProducer.offlineBeamSpotToCUDA_cfi import offlineBeamSpotToCUDA as _offlineBeamSpotToCUDA
+ process.hltOnlineBeamSpotToCUDA = _offlineBeamSpotToCUDA.clone(
+ src = "hltOnlineBeamSpot"
+ )
+
+ # reconstruct the pixel digis and clusters on the gpu
+ from RecoLocalTracker.SiPixelClusterizer.siPixelRawToClusterCUDA_cfi import siPixelRawToClusterCUDA as _siPixelRawToClusterCUDA
+ process.hltSiPixelClustersCUDA = _siPixelRawToClusterCUDA.clone()
+ # use the pixel channel calibrations scheme for Run 3
+ run3_common.toModify(process.hltSiPixelClustersCUDA, isRun2 = False)
+
+ # copy the pixel digis errors to the host
+ from EventFilter.SiPixelRawToDigi.siPixelDigiErrorsSoAFromCUDA_cfi import siPixelDigiErrorsSoAFromCUDA as _siPixelDigiErrorsSoAFromCUDA
+ process.hltSiPixelDigiErrorsSoA = _siPixelDigiErrorsSoAFromCUDA.clone(
+ src = "hltSiPixelClustersCUDA"
+ )
+
+ # convert the pixel digis errors to the legacy format
+ from EventFilter.SiPixelRawToDigi.siPixelDigiErrorsFromSoA_cfi import siPixelDigiErrorsFromSoA as _siPixelDigiErrorsFromSoA
+ process.hltSiPixelDigiErrors = _siPixelDigiErrorsFromSoA.clone(
+ digiErrorSoASrc = "hltSiPixelDigiErrorsSoA",
+ UsePhase1 = True
+ )
+
+ # copy the pixel digis (except errors) and clusters to the host
+ from EventFilter.SiPixelRawToDigi.siPixelDigisSoAFromCUDA_cfi import siPixelDigisSoAFromCUDA as _siPixelDigisSoAFromCUDA
+ process.hltSiPixelDigisSoA = _siPixelDigisSoAFromCUDA.clone(
+ src = "hltSiPixelClustersCUDA"
+ )
+
+ # convert the pixel digis (except errors) and clusters to the legacy format
+ from RecoLocalTracker.SiPixelClusterizer.siPixelDigisClustersFromSoA_cfi import siPixelDigisClustersFromSoA as _siPixelDigisClustersFromSoA
+ process.hltSiPixelDigisClusters = _siPixelDigisClustersFromSoA.clone(
+ src = "hltSiPixelDigisSoA"
+ )
+
+ # SwitchProducer wrapping the legacy pixel digis producer or an alias combining the pixel digis information converted from SoA
+ process.hltSiPixelDigis = SwitchProducerCUDA(
+ # legacy producer
+ cpu = process.hltSiPixelDigis,
+ # alias used to access products from multiple conversion modules
+ cuda = cms.EDAlias(
+ hltSiPixelDigisClusters = cms.VPSet(
+ cms.PSet(type = cms.string("PixelDigiedmDetSetVector"))
+ ),
+ hltSiPixelDigiErrors = cms.VPSet(
+ cms.PSet(type = cms.string("DetIdedmEDCollection")),
+ cms.PSet(type = cms.string("SiPixelRawDataErroredmDetSetVector")),
+ cms.PSet(type = cms.string("PixelFEDChanneledmNewDetSetVector"))
+ )
+ )
+ )
+
+ # SwitchProducer wrapping the legacy pixel cluster producer or an alias for the pixel clusters information converted from SoA
+ process.hltSiPixelClusters = SwitchProducerCUDA(
+ # legacy producer
+ cpu = process.hltSiPixelClusters,
+ # alias used to access products from multiple conversion modules
+ cuda = cms.EDAlias(
+ hltSiPixelDigisClusters = cms.VPSet(
+ cms.PSet(type = cms.string("SiPixelClusteredmNewDetSetVector"))
+ )
+ )
+ )
+
+ # reconstruct the pixel rechits on the gpu
+ from RecoLocalTracker.SiPixelRecHits.siPixelRecHitCUDA_cfi import siPixelRecHitCUDA as _siPixelRecHitCUDA
+ process.hltSiPixelRecHitsCUDA = _siPixelRecHitCUDA.clone(
+ src = "hltSiPixelClustersCUDA",
+ beamSpot = "hltOnlineBeamSpotToCUDA"
+ )
+
+ # SwitchProducer wrapping the legacy pixel rechit producer or the transfer of the pixel rechits to the host and the conversion from SoA
+ from RecoLocalTracker.SiPixelRecHits.siPixelRecHitFromCUDA_cfi import siPixelRecHitFromCUDA as _siPixelRecHitFromCUDA
+ process.hltSiPixelRecHits = SwitchProducerCUDA(
+ # legacy producer
+ cpu = process.hltSiPixelRecHits,
+ # converter to legacy format
+ cuda = _siPixelRecHitFromCUDA.clone(
+ pixelRecHitSrc = "hltSiPixelRecHitsCUDA",
+ src = "hltSiPixelClusters"
+ )
+ )
+
+
+ # Tasks and Sequences
+
+ process.HLTDoLocalPixelTask = cms.Task(
+ process.hltOnlineBeamSpotToCUDA, # transfer the beamspot to the gpu
+ process.hltSiPixelClustersCUDA, # reconstruct the pixel digis and clusters on the gpu
+ process.hltSiPixelRecHitsCUDA, # reconstruct the pixel rechits on the gpu
+ process.hltSiPixelDigisSoA, # copy the pixel digis (except errors) and clusters to the host
+ process.hltSiPixelDigisClusters, # convert the pixel digis (except errors) and clusters to the legacy format
+ process.hltSiPixelDigiErrorsSoA, # copy the pixel digis errors to the host
+ process.hltSiPixelDigiErrors, # convert the pixel digis errors to the legacy format
+ process.hltSiPixelDigis, # SwitchProducer wrapping the legacy pixel digis producer or an alias combining the pixel digis information converted from SoA
+ process.hltSiPixelClusters, # SwitchProducer wrapping the legacy pixel cluster producer or an alias for the pixel clusters information converted from SoA
+ process.hltSiPixelClustersCache, # legacy module, used by the legacy pixel quadruplet producer
+ process.hltSiPixelRecHits) # SwitchProducer wrapping the legacy pixel rechit producer or the transfer of the pixel rechits to the host and the conversion from SoA
+
+ process.HLTDoLocalPixelSequence = cms.Sequence(process.HLTDoLocalPixelTask)
+
+
+ # done
+ return process
+
+
+# customisation for running the "Patatrack" pixel track reconstruction
+def customisePixelTrackReconstruction(process):
+
+ if not 'HLTRecoPixelTracksSequence' in process.__dict__:
+ return process
+
+
+ # FIXME replace the Sequences with empty ones to avoid exanding them during the (re)definition of Modules and EDAliases
+
+ process.HLTRecoPixelTracksSequence = cms.Sequence()
+ process.HLTRecopixelvertexingSequence = cms.Sequence()
+
+
+ # Modules and EDAliases
+
+ # referenced in process.HLTRecoPixelTracksTask
+
+ # cpu only: convert the pixel rechits from legacy to SoA format
+ from RecoLocalTracker.SiPixelRecHits.siPixelRecHitSoAFromLegacy_cfi import siPixelRecHitSoAFromLegacy as _siPixelRecHitSoAFromLegacy
+ process.hltSiPixelRecHitSoA = _siPixelRecHitSoAFromLegacy.clone(
+ src = "hltSiPixelClusters",
+ beamSpot = "hltOnlineBeamSpot",
+ convertToLegacy = True
+ )
+
+ # build pixel ntuplets and pixel tracks in SoA format on gpu
+ from RecoPixelVertexing.PixelTriplets.caHitNtupletCUDA_cfi import caHitNtupletCUDA as _caHitNtupletCUDA
+ process.hltPixelTracksCUDA = _caHitNtupletCUDA.clone(
+ idealConditions = False,
+ pixelRecHitSrc = "hltSiPixelRecHitsCUDA",
+ onGPU = True
+ )
+ # use quality cuts tuned for Run 2 ideal conditions for all Run 3 workflows
+ run3_common.toModify(process.hltPixelTracksCUDA, idealConditions = True)
+
+ # SwitchProducer providing the pixel tracks in SoA format on cpu
+ process.hltPixelTracksSoA = SwitchProducerCUDA(
+ # build pixel ntuplets and pixel tracks in SoA format on cpu
+ cpu = _caHitNtupletCUDA.clone(
+ idealConditions = False,
+ pixelRecHitSrc = "hltSiPixelRecHitSoA",
+ onGPU = False
+ ),
+ # transfer the pixel tracks in SoA format to the host
+ cuda = cms.EDProducer("PixelTrackSoAFromCUDA",
+ src = cms.InputTag("hltPixelTracksCUDA")
+ )
+ )
+ # use quality cuts tuned for Run 2 ideal conditions for all Run 3 workflows
+ run3_common.toModify(process.hltPixelTracksSoA.cpu, idealConditions = True)
+
+ # convert the pixel tracks from SoA to legacy format
+ from RecoPixelVertexing.PixelTrackFitting.pixelTrackProducerFromSoA_cfi import pixelTrackProducerFromSoA as _pixelTrackProducerFromSoA
+ process.hltPixelTracks = _pixelTrackProducerFromSoA.clone(
+ beamSpot = "hltOnlineBeamSpot",
+ pixelRecHitLegacySrc = "hltSiPixelRecHits",
+ trackSrc = "hltPixelTracksSoA"
+ )
+
+
+ # referenced in process.HLTRecopixelvertexingTask
+
+ # build pixel vertices in SoA format on gpu
+ from RecoPixelVertexing.PixelVertexFinding.pixelVertexCUDA_cfi import pixelVertexCUDA as _pixelVertexCUDA
+ process.hltPixelVerticesCUDA = _pixelVertexCUDA.clone(
+ pixelTrackSrc = "hltPixelTracksCUDA",
+ onGPU = True
+ )
+
+ # build or transfer pixel vertices in SoA format on cpu
+ process.hltPixelVerticesSoA = SwitchProducerCUDA(
+ # build pixel vertices in SoA format on cpu
+ cpu = _pixelVertexCUDA.clone(
+ pixelTrackSrc = "hltPixelTracksSoA",
+ onGPU = False
+ ),
+ # transfer the pixel vertices in SoA format to cpu
+ cuda = cms.EDProducer("PixelVertexSoAFromCUDA",
+ src = cms.InputTag("hltPixelVerticesCUDA")
+ )
+ )
+
+ # convert the pixel vertices from SoA to legacy format
+ from RecoPixelVertexing.PixelVertexFinding.pixelVertexFromSoA_cfi import pixelVertexFromSoA as _pixelVertexFromSoA
+ process.hltPixelVertices = _pixelVertexFromSoA.clone(
+ src = "hltPixelVerticesSoA",
+ TrackCollection = "hltPixelTracks",
+ beamSpot = "hltOnlineBeamSpot"
+ )
+
+
+ # Tasks and Sequences
+
+ process.HLTRecoPixelTracksTask = cms.Task(
+ process.hltPixelTracksTrackingRegions, # from the original sequence
+ process.hltSiPixelRecHitSoA, # pixel rechits on cpu, converted to SoA
+ process.hltPixelTracksCUDA, # pixel ntuplets on gpu, in SoA format
+ process.hltPixelTracksSoA, # pixel ntuplets on cpu, in SoA format
+ process.hltPixelTracks) # pixel tracks on cpu, in legacy format
+
+
+ process.HLTRecoPixelTracksSequence = cms.Sequence(process.HLTRecoPixelTracksTask)
+
+ process.HLTRecopixelvertexingTask = cms.Task(
+ process.HLTRecoPixelTracksTask,
+ process.hltPixelVerticesCUDA, # pixel vertices on gpu, in SoA format
+ process.hltPixelVerticesSoA, # pixel vertices on cpu, in SoA format
+ process.hltPixelVertices, # pixel vertices on cpu, in legacy format
+ process.hltTrimmedPixelVertices) # from the original sequence
+
+ process.HLTRecopixelvertexingSequence = cms.Sequence(
+ process.hltPixelTracksFitter + # not used here, kept for compatibility with legacy sequences
+ process.hltPixelTracksFilter, # not used here, kept for compatibility with legacy sequences
+ process.HLTRecopixelvertexingTask)
+
+
+ # done
+ return process
+
+
+# customisation for offloading the ECAL local reconstruction via CUDA if a supported gpu is present
+def customiseEcalLocalReconstruction(process):
+
+ if not 'HLTDoFullUnpackingEgammaEcalSequence' in process.__dict__:
+ return process
+
+
+ # FIXME replace the Sequences with empty ones to avoid exanding them during the (re)definition of Modules and EDAliases
+
+ process.HLTDoFullUnpackingEgammaEcalMFSequence = cms.Sequence()
+ process.HLTDoFullUnpackingEgammaEcalWithoutPreshowerSequence = cms.Sequence()
+ process.HLTDoFullUnpackingEgammaEcalSequence = cms.Sequence()
+
+
+ # Event Setup
+
+ process.load("EventFilter.EcalRawToDigi.ecalElectronicsMappingGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalGainRatiosGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalPedestalsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalPulseCovariancesGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalPulseShapesGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalSamplesCorrelationGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalTimeBiasCorrectionsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalTimeCalibConstantsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalMultifitParametersGPUESProducer_cfi")
+
+ process.load("RecoLocalCalo.EcalRecProducers.ecalRechitADCToGeVConstantGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalRechitChannelStatusGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalIntercalibConstantsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalLaserAPDPNRatiosGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalLaserAPDPNRatiosRefGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalLaserAlphasGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalLinearCorrectionsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.EcalRecProducers.ecalRecHitParametersGPUESProducer_cfi")
+
+
+ # Modules and EDAliases
+
+ # ECAL unpacker running on gpu
+ process.hltEcalDigisGPU = cms.EDProducer("EcalRawToDigiGPU",
+ InputLabel = cms.InputTag("rawDataCollector"),
+ FEDs = cms.vint32(
+ 601, 602, 603, 604, 605,
+ 606, 607, 608, 609, 610,
+ 611, 612, 613, 614, 615,
+ 616, 617, 618, 619, 620,
+ 621, 622, 623, 624, 625,
+ 626, 627, 628, 629, 630,
+ 631, 632, 633, 634, 635,
+ 636, 637, 638, 639, 640,
+ 641, 642, 643, 644, 645,
+ 646, 647, 648, 649, 650,
+ 651, 652, 653, 654
+ ),
+ digisLabelEB = cms.string("ebDigis"),
+ digisLabelEE = cms.string("eeDigis"),
+ maxChannelsEB = cms.uint32(61200),
+ maxChannelsEE = cms.uint32(14648),
+ )
+
+ # SwitchProducer wrapping the legacy ECAL unpacker or the ECAL digi converter from SoA format on gpu to legacy format on cpu
+ process.hltEcalDigisLegacy = process.hltEcalDigis.clone()
+
+ process.hltEcalDigis = SwitchProducerCUDA(
+ # legacy producer
+ cpu = cms.EDAlias(
+ hltEcalDigisLegacy = cms.VPSet(
+ cms.PSet(type = cms.string("EBDigiCollection")),
+ cms.PSet(type = cms.string("EEDigiCollection")),
+ cms.PSet(type = cms.string("EBDetIdedmEDCollection")),
+ cms.PSet(type = cms.string("EEDetIdedmEDCollection")),
+ cms.PSet(type = cms.string("EBSrFlagsSorted")),
+ cms.PSet(type = cms.string("EESrFlagsSorted")),
+ cms.PSet(type = cms.string("EcalElectronicsIdedmEDCollection"), fromProductInstance = cms.string("EcalIntegrityBlockSizeErrors")),
+ cms.PSet(type = cms.string("EcalElectronicsIdedmEDCollection"), fromProductInstance = cms.string("EcalIntegrityTTIdErrors"))
+ )
+ ),
+ # convert ECAL digis from SoA format on gpu to legacy format on cpu
+ cuda = cms.EDProducer("EcalCPUDigisProducer",
+ digisInLabelEB = cms.InputTag("hltEcalDigisGPU", "ebDigis"),
+ digisInLabelEE = cms.InputTag("hltEcalDigisGPU", "eeDigis"),
+ digisOutLabelEB = cms.string("ebDigis"),
+ digisOutLabelEE = cms.string("eeDigis"),
+ produceDummyIntegrityCollections = cms.bool(True)
+ )
+ )
+
+ # ECAL multifit running on gpu
+ from RecoLocalCalo.EcalRecProducers.ecalUncalibRecHitProducerGPU_cfi import ecalUncalibRecHitProducerGPU as _ecalUncalibRecHitProducerGPU
+ process.hltEcalUncalibRecHitGPU = _ecalUncalibRecHitProducerGPU.clone(
+ digisLabelEB = ("hltEcalDigisGPU", "ebDigis"),
+ digisLabelEE = ("hltEcalDigisGPU", "eeDigis"),
+ shouldRunTimingComputation = False
+ )
+
+ # copy the ECAL uncalibrated rechits from gpu to cpu in SoA format
+ process.hltEcalUncalibRecHitSoA = cms.EDProducer("EcalCPUUncalibRecHitProducer",
+ containsTimingInformation = cms.bool(False),
+ recHitsInLabelEB = cms.InputTag("hltEcalUncalibRecHitGPU", "EcalUncalibRecHitsEB"),
+ recHitsInLabelEE = cms.InputTag("hltEcalUncalibRecHitGPU", "EcalUncalibRecHitsEE"),
+ recHitsOutLabelEB = cms.string("EcalUncalibRecHitsEB"),
+ recHitsOutLabelEE = cms.string("EcalUncalibRecHitsEE")
+ )
+
+ # SwitchProducer wrapping the legacy ECAL uncalibrated rechits producer or a converter from SoA to legacy format
+ process.hltEcalUncalibRecHit = SwitchProducerCUDA(
+ # legacy producer
+ cpu = process.hltEcalUncalibRecHit,
+ # convert the ECAL uncalibrated rechits from SoA to legacy format
+ cuda = cms.EDProducer("EcalUncalibRecHitConvertGPU2CPUFormat",
+ recHitsLabelGPUEB = cms.InputTag("hltEcalUncalibRecHitSoA", "EcalUncalibRecHitsEB"),
+ recHitsLabelGPUEE = cms.InputTag("hltEcalUncalibRecHitSoA", "EcalUncalibRecHitsEE"),
+ recHitsLabelCPUEB = cms.string("EcalUncalibRecHitsEB"),
+ recHitsLabelCPUEE = cms.string("EcalUncalibRecHitsEE")
+ )
+ )
+
+ # Reconstructing the ECAL calibrated rechits on gpu works, but is extremely slow.
+ # Disable it for the time being, until the performance has been addressed.
+ """
+ process.hltEcalRecHitGPU = cms.EDProducer("EcalRecHitProducerGPU",
+ uncalibrecHitsInLabelEB = cms.InputTag("hltEcalUncalibRecHitGPU","EcalUncalibRecHitsEB"),
+ uncalibrecHitsInLabelEE = cms.InputTag("hltEcalUncalibRecHitGPU","EcalUncalibRecHitsEE"),
+ recHitsLabelEB = cms.string("EcalRecHitsEB"),
+ recHitsLabelEE = cms.string("EcalRecHitsEE"),
+ maxNumberHitsEB = cms.uint32(61200),
+ maxNumberHitsEE = cms.uint32(14648),
+ ChannelStatusToBeExcluded = cms.vstring(
+ "kDAC",
+ "kNoisy",
+ "kNNoisy",
+ "kFixedG6",
+ "kFixedG1",
+ "kFixedG0",
+ "kNonRespondingIsolated",
+ "kDeadVFE",
+ "kDeadFE",
+ "kNoDataNoTP"),
+ killDeadChannels = cms.bool(True),
+ EBLaserMIN = cms.double(0.01),
+ EELaserMIN = cms.double(0.01),
+ EBLaserMAX = cms.double(30.0),
+ EELaserMAX = cms.double(30.0),
+ flagsMapDBReco = cms.PSet(
+ kGood = cms.vstring("kOk","kDAC","kNoLaser","kNoisy"),
+ kNoisy = cms.vstring("kNNoisy","kFixedG6","kFixedG1"),
+ kNeighboursRecovered = cms.vstring("kFixedG0", "kNonRespondingIsolated", "kDeadVFE"),
+ kTowerRecovered = cms.vstring("kDeadFE"),
+ kDead = cms.vstring("kNoDataNoTP")
+ ),
+ recoverEBIsolatedChannels = cms.bool(False),
+ recoverEEIsolatedChannels = cms.bool(False),
+ recoverEBVFE = cms.bool(False),
+ recoverEEVFE = cms.bool(False),
+ recoverEBFE = cms.bool(True),
+ recoverEEFE = cms.bool(True),
+ )
+
+ process.hltEcalRecHitSoA = cms.EDProducer("EcalCPURecHitProducer",
+ recHitsInLabelEB = cms.InputTag("hltEcalRecHitGPU", "EcalRecHitsEB"),
+ recHitsInLabelEE = cms.InputTag("hltEcalRecHitGPU", "EcalRecHitsEE"),
+ recHitsOutLabelEB = cms.string("EcalRecHitsEB"),
+ recHitsOutLabelEE = cms.string("EcalRecHitsEE"),
+ containsTimingInformation = cms.bool(False),
+ )
+
+ # SwitchProducer wrapping the legacy ECAL calibrated rechits producer or a converter from SoA to legacy format
+ process.hltEcalRecHit = SwitchProducerCUDA(
+ # legacy producer
+ cpu = process.hltEcalRecHit,
+ # convert the ECAL calibrated rechits from SoA to legacy format
+ cuda = cms.EDProducer("EcalRecHitConvertGPU2CPUFormat",
+ recHitsLabelGPUEB = cms.InputTag("hltEcalRecHitSoA", "EcalRecHitsEB"),
+ recHitsLabelGPUEE = cms.InputTag("hltEcalRecHitSoA", "EcalRecHitsEE"),
+ recHitsLabelCPUEB = cms.string("EcalRecHitsEB"),
+ recHitsLabelCPUEE = cms.string("EcalRecHitsEE"),
+ )
+ """
+
+
+ # SwitchProducer wrapping the legacy ECAL rechits producer
+ # the gpu unpacker does not produce the TPs used for the recovery, so the SwitchProducer alias does not provide them:
+ # - the cpu uncalibrated rechit producer may mark them for recovery, read the TPs explicitly from the legacy unpacker
+ # - the gpu uncalibrated rechit producer does not flag them for recovery, so the TPs are not necessary
+ process.hltEcalRecHit = SwitchProducerCUDA(
+ cpu = process.hltEcalRecHit.clone(
+ triggerPrimitiveDigiCollection = cms.InputTag('hltEcalDigisLegacy', 'EcalTriggerPrimitives')
+ ),
+ cuda = process.hltEcalRecHit.clone(
+ triggerPrimitiveDigiCollection = cms.InputTag('unused')
+ )
+ )
+
+ # Tasks and Sequences
+
+ process.HLTDoFullUnpackingEgammaEcalWithoutPreshowerTask = cms.Task(
+ process.hltEcalDigisGPU, # unpack ECAL digis on gpu
+ process.hltEcalDigisLegacy, # legacy producer, referenced in the SwitchProducer
+ process.hltEcalDigis, # SwitchProducer
+ process.hltEcalUncalibRecHitGPU, # run ECAL local reconstruction and multifit on gpu
+ process.hltEcalUncalibRecHitSoA, # needed by hltEcalPhiSymFilter - copy to host
+ process.hltEcalUncalibRecHit, # needed by hltEcalPhiSymFilter - convert to legacy format
+ # process.hltEcalRecHitGPU, # make ECAL calibrated rechits on gpu
+ # process.hltEcalRecHitSoA, # copy to host
+ process.hltEcalDetIdToBeRecovered, # legacy producer
+ process.hltEcalRecHit) # legacy producer
+
+ process.HLTDoFullUnpackingEgammaEcalWithoutPreshowerSequence = cms.Sequence(
+ process.HLTDoFullUnpackingEgammaEcalWithoutPreshowerTask)
+
+ process.HLTPreshowerTask = cms.Task(
+ process.hltEcalPreshowerDigis, # unpack ECAL preshower digis on the host
+ process.hltEcalPreshowerRecHit) # build ECAL preshower rechits on the host
+
+ process.HLTPreshowerSequence = cms.Sequence(process.HLTPreshowerTask)
+
+ process.HLTDoFullUnpackingEgammaEcalTask = cms.Task(
+ process.HLTDoFullUnpackingEgammaEcalWithoutPreshowerTask,
+ process.HLTPreshowerTask)
+
+ process.HLTDoFullUnpackingEgammaEcalSequence = cms.Sequence(
+ process.HLTDoFullUnpackingEgammaEcalTask)
+
+ process.HLTDoFullUnpackingEgammaEcalMFSequence = cms.Sequence(
+ process.HLTDoFullUnpackingEgammaEcalTask)
+
+
+ # done
+ return process
+
+# customisation for offloading the HCAL local reconstruction via CUDA if a supported gpu is present
+def customiseHcalLocalReconstruction(process):
+
+ if not 'HLTDoLocalHcalSequence' in process.__dict__:
+ return process
+
+
+ # FIXME replace the Sequences with empty ones to avoid exanding them during the (re)definition of Modules and EDAliases
+
+ process.HLTDoLocalHcalSequence = cms.Sequence()
+ process.HLTStoppedHSCPLocalHcalReco = cms.Sequence()
+
+
+ # Event Setup
+
+ process.load("EventFilter.HcalRawToDigi.hcalElectronicsMappingGPUESProducer_cfi")
+
+ process.load("RecoLocalCalo.HcalRecProducers.hcalGainsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalGainWidthsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalLUTCorrsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalConvertedPedestalsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalConvertedEffectivePedestalsGPUESProducer_cfi")
+ process.hcalConvertedEffectivePedestalsGPUESProducer.label0 = "withTopoEff"
+ process.load("RecoLocalCalo.HcalRecProducers.hcalConvertedPedestalWidthsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalConvertedEffectivePedestalWidthsGPUESProducer_cfi")
+ process.hcalConvertedEffectivePedestalWidthsGPUESProducer.label0 = "withTopoEff"
+ process.hcalConvertedEffectivePedestalWidthsGPUESProducer.label1 = "withTopoEff"
+ process.load("RecoLocalCalo.HcalRecProducers.hcalQIECodersGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalRecoParamsWithPulseShapesGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalRespCorrsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalTimeCorrsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalQIETypesGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalSiPMParametersGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalSiPMCharacteristicsGPUESProducer_cfi")
+ process.load("RecoLocalCalo.HcalRecProducers.hcalMahiPulseOffsetsGPUESProducer_cfi")
+
+
+ # Modules and EDAliases
+
+ # The HCAL unpacker running on the gpu supports only the HB and HE digis.
+ # So, run the legacy unacker on the cpu, then convert the HB and HE digis
+ # to SoA format and copy them to the gpu.
+ process.hltHcalDigisGPU = cms.EDProducer("HcalDigisProducerGPU",
+ hbheDigisLabel = cms.InputTag("hltHcalDigis"),
+ qie11DigiLabel = cms.InputTag("hltHcalDigis"),
+ digisLabelF01HE = cms.string(""),
+ digisLabelF5HB = cms.string(""),
+ digisLabelF3HB = cms.string(""),
+ maxChannelsF01HE = cms.uint32(10000),
+ maxChannelsF5HB = cms.uint32(10000),
+ maxChannelsF3HB = cms.uint32(10000)
+ )
+
+ # run the HCAL local reconstruction (including Method 0 and MAHI) on gpu
+ from RecoLocalCalo.HcalRecProducers.hbheRecHitProducerGPU_cfi import hbheRecHitProducerGPU as _hbheRecHitProducerGPU
+ process.hltHbherecoGPU = _hbheRecHitProducerGPU.clone(
+ digisLabelF01HE = "hltHcalDigisGPU",
+ digisLabelF5HB = "hltHcalDigisGPU",
+ digisLabelF3HB = "hltHcalDigisGPU",
+ recHitsLabelM0HBHE = ""
+ )
+
+ # transfer the HCAL rechits to the cpu, and convert them to the legacy format
+ from RecoLocalCalo.HcalRecProducers.hcalCPURecHitsProducer_cfi import hcalCPURecHitsProducer as _hcalCPURecHitsProducer
+ process.hltHbherecoFromGPU = _hcalCPURecHitsProducer.clone(
+ recHitsM0LabelIn = "hltHbherecoGPU",
+ recHitsM0LabelOut = "",
+ recHitsLegacyLabelOut = ""
+ )
+
+ # SwitchProducer between the legacy producer and the copy from gpu with conversion
+ process.hltHbhereco = SwitchProducerCUDA(
+ # legacy producer
+ cpu = process.hltHbhereco.clone(),
+ # alias to the rechits converted to legacy format
+ cuda = cms.EDAlias(
+ hltHbherecoFromGPU = cms.VPSet(
+ cms.PSet(type = cms.string("HBHERecHitsSorted"))
+ )
+ )
+ )
+
+
+ # Tasks and Sequences
+
+ process.HLTDoLocalHcalTask = cms.Task(
+ process.hltHcalDigis, # legacy producer, unpack HCAL digis on cpu
+ process.hltHcalDigisGPU, # copy to gpu and convert to SoA format
+ process.hltHbherecoGPU, # run the HCAL local reconstruction (including Method 0 and MAHI) on gpu
+ process.hltHbherecoFromGPU, # transfer the HCAL rechits to the cpu, and convert them to the legacy format
+ process.hltHbhereco, # SwitchProducer between the legacy producer and the copy from gpu with conversion
+ process.hltHfprereco, # legacy producer
+ process.hltHfreco, # legacy producer
+ process.hltHoreco) # legacy producer
+
+ process.HLTDoLocalHcalSequence = cms.Sequence(
+ process.HLTDoLocalHcalTask)
+
+ process.HLTStoppedHSCPLocalHcalRecoTask = cms.Task(
+ process.hltHcalDigis, # legacy producer, unpack HCAL digis on cpu
+ process.hltHcalDigisGPU, # copy to gpu and convert to SoA format
+ process.hltHbherecoGPU, # run the HCAL local reconstruction (including Method 0 and MAHI) on gpu
+ process.hltHbherecoFromGPU, # transfer the HCAL rechits to the cpu, and convert them to the legacy format
+ process.hltHbhereco) # SwitchProducer between the legacy producer and the copy from gpu with conversion
+
+ process.HLTStoppedHSCPLocalHcalReco = cms.Sequence(
+ process.HLTStoppedHSCPLocalHcalRecoTask)
+
+
+ # done
+ return process
+
+
+# customisation to enable pixel triplets instead of quadruplets
+def enablePatatrackPixelTriplets(process):
+
+ # configure GPU pixel tracks for triplets
+ process.hltPixelTracksCUDA.minHitsPerNtuplet = 3
+ process.hltPixelTracksCUDA.includeJumpingForwardDoublets = True
+
+ # configure CPU pixel tracks for triplets
+ process.hltPixelTracksSoA.cpu.minHitsPerNtuplet = 3
+ process.hltPixelTracksSoA.cpu.includeJumpingForwardDoublets = True
+
+ # done
+ return process
+
+
+# customisation for running the Patatrack reconstruction, with automatic offload via CUDA when a supported gpu is available
+def customizeHLTforPatatrack(process):
+ process = customiseCommon(process)
+ process = customisePixelLocalReconstruction(process)
+ process = customisePixelTrackReconstruction(process)
+ process = customiseEcalLocalReconstruction(process)
+ process = customiseHcalLocalReconstruction(process)
+ return process
+
+
+# customisation for running the Patatrack triplets reconstruction, with automatic offload via CUDA when a supported gpu is available
+def customizeHLTforPatatrackTriplets(process):
+ process = customiseCommon(process)
+ process = customisePixelLocalReconstruction(process)
+ process = customisePixelTrackReconstruction(process)
+ process = customiseEcalLocalReconstruction(process)
+ process = customiseHcalLocalReconstruction(process)
+ process = enablePatatrackPixelTriplets(process)
+ return process
+
+
+def _addConsumerPath(process):
+ # add to a path all consumers and the tasks that define the producers
+ process.Consumer = cms.Path(
+ process.HLTBeginSequence +
+ process.hltPixelConsumer +
+ process.hltEcalConsumer +
+ process.hltHbheConsumer,
+ process.HLTDoLocalPixelTask,
+ process.HLTRecoPixelTracksTask,
+ process.HLTRecopixelvertexingTask,
+ process.HLTDoFullUnpackingEgammaEcalTask,
+ process.HLTDoLocalHcalTask,
+ )
+
+ if 'HLTSchedule' in process.__dict__:
+ process.HLTSchedule.append(process.Consumer)
+ if process.schedule is not None:
+ process.schedule.append(process.Consumer)
+
+ # done
+ return process
+
+
+def consumeGPUSoAProducts(process):
+ # consume the Pixel tracks and vertices on the GPU in SoA format
+ process.hltPixelConsumer = cms.EDAnalyzer("GenericConsumer",
+ eventProducts = cms.untracked.vstring( 'hltPixelTracksCUDA', 'hltPixelVerticesCUDA' )
+ )
+
+ # consume the ECAL uncalibrated rechits on the GPU in SoA format
+ process.hltEcalConsumer = cms.EDAnalyzer("GenericConsumer",
+ eventProducts = cms.untracked.vstring( 'hltEcalUncalibRecHitGPU' )
+ )
+
+ # consume the HCAL rechits on the GPU in SoA format
+ process.hltHbheConsumer = cms.EDAnalyzer("GenericConsumer",
+ eventProducts = cms.untracked.vstring( 'hltHbherecoGPU' )
+ )
+
+ # add to a path all consumers and the tasks that define the producers
+ process = _addConsumerPath(process)
+
+ # done
+ return process
+
+
+def consumeCPUSoAProducts(process):
+ # consume the Pixel tracks and vertices on the CPU in SoA format
+ process.hltPixelConsumer = cms.EDAnalyzer("GenericConsumer",
+ eventProducts = cms.untracked.vstring( 'hltPixelTracksSoA', 'hltPixelVerticesSoA' )
+ )
+
+ # consume the ECAL uncalibrated rechits on the CPU in SoA format
+ process.hltEcalConsumer = cms.EDAnalyzer("GenericConsumer",
+ eventProducts = cms.untracked.vstring( 'hltEcalUncalibRecHitSoA' )
+ )
+
+ # consume the HCAL rechits on the CPU in legacy format
+ process.hltHbheConsumer = cms.EDAnalyzer("GenericConsumer",
+ eventProducts = cms.untracked.vstring( 'hltHbhereco' )
+ )
+
+ # add to a path all consumers and the tasks that define the producers
+ process = _addConsumerPath(process)
+
+ # done
+ return process
+
+def consumeCPULegacyProducts(process):
+ # consume the Pixel tracks and vertices on the CPU in legacy format
+ process.hltPixelConsumer = cms.EDAnalyzer("GenericConsumer",
+ eventProducts = cms.untracked.vstring( 'hltPixelTracks', 'hltPixelVertices' )
+ )
+
+ # consume the ECAL runcalibrated echits on the CPU in legacy format
+ process.hltEcalConsumer = cms.EDAnalyzer("GenericConsumer",
+ eventProducts = cms.untracked.vstring( 'hltEcalUncalibRecHit' )
+ )
+
+ # consume the HCAL rechits on the CPU in legacy format
+ process.hltHbheConsumer = cms.EDAnalyzer("GenericConsumer",
+ eventProducts = cms.untracked.vstring( 'hltHbhereco' )
+ )
+
+ # add to a path all consumers and the tasks that define the producers
+ process = _addConsumerPath(process)
+
+ # done
+ return process
diff --git a/RecoPixelVertexing/Configuration/python/RecoPixelVertexing_cff.py b/RecoPixelVertexing/Configuration/python/RecoPixelVertexing_cff.py
index 34ee6fadb04de..424ac13a43627 100644
--- a/RecoPixelVertexing/Configuration/python/RecoPixelVertexing_cff.py
+++ b/RecoPixelVertexing/Configuration/python/RecoPixelVertexing_cff.py
@@ -4,7 +4,21 @@
#
# for STARTUP ONLY use try and use Offline 3D PV from pixelTracks, with adaptive vertex
#
-#from RecoPixelVertexing.PixelVertexFinding.PixelVertexes_cff import *
-from RecoVertex.PrimaryVertexProducer.OfflinePixel3DPrimaryVertices_cfi import *
+from RecoPixelVertexing.PixelVertexFinding.PixelVertexes_cff import *
+#from RecoVertex.PrimaryVertexProducer.OfflinePixel3DPrimaryVertices_cfi import *
recopixelvertexingTask = cms.Task(pixelTracksTask,pixelVertices)
recopixelvertexing = cms.Sequence(recopixelvertexingTask)
+
+from Configuration.ProcessModifiers.gpu_cff import gpu
+
+from RecoPixelVertexing.PixelVertexFinding.pixelVertexCUDA_cfi import pixelVertexCUDA
+from RecoPixelVertexing.PixelVertexFinding.pixelVertexSoA_cfi import pixelVertexSoA
+from RecoPixelVertexing.PixelVertexFinding.pixelVertexFromSoA_cfi import pixelVertexFromSoA as _pixelVertexFromSoA
+
+_pixelVertexingCUDATask = cms.Task(pixelTracksTask,pixelVertexCUDA,pixelVertexSoA,pixelVertices)
+
+# pixelVertexSoAonCPU = pixelVertexCUDA.clone()
+# pixelVertexSoAonCPU.onGPU = False;
+
+gpu.toReplaceWith(pixelVertices,_pixelVertexFromSoA)
+gpu.toReplaceWith(recopixelvertexingTask,_pixelVertexingCUDATask)
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..909959f2d81be
--- /dev/null
+++ b/RecoPixelVertexing/Configuration/python/customizePixelTracksSoAonCPU.py
@@ -0,0 +1,61 @@
+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.PixelVertexFinding.pixelVertexCUDA_cfi import pixelVertexCUDA
+ process.pixelVertexSoA = pixelVertexCUDA.clone(
+ onGPU = False,
+ pixelTrackSrc = 'pixelTrackSoA'
+ )
+
+ from RecoPixelVertexing.PixelTrackFitting.pixelTrackProducerFromSoA_cfi import pixelTrackProducerFromSoA
+ process.pixelTracks = pixelTrackProducerFromSoA.clone(
+ pixelRecHitLegacySrc = 'siPixelRecHitsPreSplitting'
+ )
+
+ from RecoPixelVertexing.PixelVertexFinding.pixelVertexFromSoA_cfi import pixelVertexFromSoA
+ process.pixelVertices = pixelVertexFromSoA.clone()
+
+ process.reconstruction_step += process.siPixelRecHitsPreSplitting + process.pixelTrackSoA + process.pixelVertexSoA
+
+ 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.MessageLogger.cerr.FwkReport.reportEvery = 100
+
+ process = customizePixelTracksSoAonCPU(process)
+ process.siPixelRecHitSoAFromLegacy.convertToLegacy = False
+
+ process.TkSoA = cms.Path(process.offlineBeamSpot + process.siPixelDigis + process.siPixelClustersPreSplitting + process.siPixelRecHitSoAFromLegacy + process.pixelTrackSoA + process.pixelVertexSoA)
+ 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/PixelTrackDumpCUDA.cc b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackDumpCUDA.cc
new file mode 100644
index 0000000000000..2f0965be50eb8
--- /dev/null
+++ b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackDumpCUDA.cc
@@ -0,0 +1,86 @@
+#include
+
+#include "CUDADataFormats/Common/interface/Product.h"
+#include "CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h"
+#include "CUDADataFormats/TrackingRecHit/interface/TrackingRecHit2DHeterogeneous.h"
+#include "CUDADataFormats/Vertex/interface/ZVertexHeterogeneous.h"
+#include "DataFormats/Common/interface/Handle.h"
+#include "FWCore/Framework/interface/ConsumesCollector.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/EDAnalyzer.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/PluginManager/interface/ModuleDef.h"
+#include "FWCore/Utilities/interface/EDGetToken.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/Utilities/interface/RunningAverage.h"
+#include "HeterogeneousCore/CUDACore/interface/ScopedContext.h"
+#include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
+
+class PixelTrackDumpCUDA : public edm::global::EDAnalyzer<> {
+public:
+ explicit PixelTrackDumpCUDA(const edm::ParameterSet& iConfig);
+ ~PixelTrackDumpCUDA() override = default;
+
+ static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+ void analyze(edm::StreamID streamID, edm::Event const& iEvent, const edm::EventSetup& iSetup) const override;
+ const bool m_onGPU;
+ edm::EDGetTokenT> tokenGPUTrack_;
+ edm::EDGetTokenT> tokenGPUVertex_;
+ edm::EDGetTokenT tokenSoATrack_;
+ edm::EDGetTokenT tokenSoAVertex_;
+};
+
+PixelTrackDumpCUDA::PixelTrackDumpCUDA(const edm::ParameterSet& iConfig)
+ : m_onGPU(iConfig.getParameter("onGPU")) {
+ if (m_onGPU) {
+ tokenGPUTrack_ =
+ consumes>(iConfig.getParameter("pixelTrackSrc"));
+ tokenGPUVertex_ =
+ consumes>(iConfig.getParameter("pixelVertexSrc"));
+ } else {
+ tokenSoATrack_ = consumes(iConfig.getParameter("pixelTrackSrc"));
+ tokenSoAVertex_ = consumes(iConfig.getParameter("pixelVertexSrc"));
+ }
+}
+
+void PixelTrackDumpCUDA::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+ edm::ParameterSetDescription desc;
+
+ desc.add("onGPU", true);
+ desc.add("pixelTrackSrc", edm::InputTag("caHitNtupletCUDA"));
+ desc.add("pixelVertexSrc", edm::InputTag("pixelVertexCUDA"));
+ descriptions.add("pixelTrackDumpCUDA", desc);
+}
+
+void PixelTrackDumpCUDA::analyze(edm::StreamID streamID,
+ edm::Event const& iEvent,
+ const edm::EventSetup& iSetup) const {
+ if (m_onGPU) {
+ auto const& hTracks = iEvent.get(tokenGPUTrack_);
+ cms::cuda::ScopedContextProduce ctx{hTracks};
+
+ auto const& tracks = ctx.get(hTracks);
+ auto const* tsoa = tracks.get();
+ assert(tsoa);
+
+ auto const& vertices = ctx.get(iEvent.get(tokenGPUVertex_));
+ auto const* vsoa = vertices.get();
+ assert(vsoa);
+
+ } else {
+ auto const* tsoa = iEvent.get(tokenSoATrack_).get();
+ assert(tsoa);
+
+ auto const* vsoa = iEvent.get(tokenSoAVertex_).get();
+ assert(vsoa);
+ }
+}
+
+DEFINE_FWK_MODULE(PixelTrackDumpCUDA);
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