CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
BrokenLineFitOnGPU.h
Go to the documentation of this file.
1 //
2 // Author: Felice Pantaleo, CERN
3 //
4 
5 // #define BROKENLINE_DEBUG
6 
7 #include <cstdint>
8 
9 #include <cuda_runtime.h>
10 
16 
17 #include "HelixFitOnGPU.h"
18 
24 
25 // #define BL_DUMP_HITS
26 
27 template <int N>
28 __global__ void kernel_BLFastFit(Tuples const *__restrict__ foundNtuplets,
30  HitsOnGPU const *__restrict__ hhp,
31  tindex_type *__restrict__ ptkids,
32  double *__restrict__ phits,
33  float *__restrict__ phits_ge,
34  double *__restrict__ pfast_fit,
35  uint32_t nHitsL,
36  uint32_t nHitsH,
37  int32_t offset) {
38  constexpr uint32_t hitsInFit = N;
39 
40  assert(hitsInFit <= nHitsL);
41  assert(nHitsL <= nHitsH);
42  assert(hhp);
43  assert(phits);
44  assert(pfast_fit);
45  assert(foundNtuplets);
46  assert(tupleMultiplicity);
47 
48  // look in bin for this hit multiplicity
50  int totTK = tupleMultiplicity->end(nHitsH) - tupleMultiplicity->begin(nHitsL);
51  assert(totTK <= int(tupleMultiplicity->size()));
52  assert(totTK >= 0);
53 
54 #ifdef BROKENLINE_DEBUG
55  if (0 == local_start) {
56  printf("%d total Ntuple\n", tupleMultiplicity->size());
57  printf("%d Ntuple of size %d/%d for %d hits to fit\n", totTK, nHitsL, nHitsH, hitsInFit);
58  }
59 #endif
60 
61  for (int local_idx = local_start, nt = riemannFit::maxNumberOfConcurrentFits; local_idx < nt;
62  local_idx += gridDim.x * blockDim.x) {
63  int tuple_idx = local_idx + offset;
64  if (tuple_idx >= totTK) {
65  ptkids[local_idx] = invalidTkId;
66  break;
67  }
68  // get it from the ntuple container (one to one to helix)
69  auto tkid = *(tupleMultiplicity->begin(nHitsL) + tuple_idx);
70  assert(tkid < foundNtuplets->nOnes());
71 
72  ptkids[local_idx] = tkid;
73 
74  auto nHits = foundNtuplets->size(tkid);
75 
76  assert(nHits >= nHitsL);
77  assert(nHits <= nHitsH);
78 
79  riemannFit::Map3xNd<N> hits(phits + local_idx);
80  riemannFit::Map4d fast_fit(pfast_fit + local_idx);
81  riemannFit::Map6xNf<N> hits_ge(phits_ge + local_idx);
82 
83 #ifdef BL_DUMP_HITS
84  __shared__ int done;
85  done = 0;
86  __syncthreads();
87  bool dump = (foundNtuplets->size(tkid) == 5 && 0 == atomicAdd(&done, 1));
88 #endif
89 
90  // Prepare data structure
91  auto const *hitId = foundNtuplets->begin(tkid);
92 
93  // #define YERR_FROM_DC
94 #ifdef YERR_FROM_DC
95  // try to compute more precise error in y
96  auto dx = hhp->xGlobal(hitId[hitsInFit - 1]) - hhp->xGlobal(hitId[0]);
97  auto dy = hhp->yGlobal(hitId[hitsInFit - 1]) - hhp->yGlobal(hitId[0]);
98  auto dz = hhp->zGlobal(hitId[hitsInFit - 1]) - hhp->zGlobal(hitId[0]);
99  float ux, uy, uz;
100 #endif
101 
102  float incr = std::max(1.f, float(nHits) / float(hitsInFit));
103  float n = 0;
104  for (uint32_t i = 0; i < hitsInFit; ++i) {
105  int j = int(n + 0.5f); // round
106  if (hitsInFit - 1 == i)
107  j = nHits - 1; // force last hit to ensure max lever arm.
108  assert(j < int(nHits));
109  n += incr;
110  auto hit = hitId[j];
111  float ge[6];
112 
113 #ifdef YERR_FROM_DC
114  auto const &dp = hhp->cpeParams().detParams(hhp->detectorIndex(hit));
115  auto status = hhp->status(hit);
116  int qbin = CPEFastParametrisation::kGenErrorQBins - 1 - status.qBin;
117  assert(qbin >= 0 && qbin < 5);
118  bool nok = (status.isBigY | status.isOneY);
119  // compute cotanbeta and use it to recompute error
120  dp.frame.rotation().multiply(dx, dy, dz, ux, uy, uz);
121  auto cb = std::abs(uy / uz);
122  int bin =
123  int(cb * (float(phase1PixelTopology::pixelThickess) / float(phase1PixelTopology::pixelPitchY)) * 8.f) - 4;
124  int low_value = 0;
125  int high_value = CPEFastParametrisation::kNumErrorBins - 1;
126  // return estimated bin value truncated to [0, 15]
127  bin = std::clamp(bin, low_value, high_value);
128  float yerr = dp.sigmay[bin] * 1.e-4f; // toCM
129  yerr *= dp.yfact[qbin]; // inflate
130  yerr *= yerr;
131  yerr += dp.apeYY;
132  yerr = nok ? hhp->yerrLocal(hit) : yerr;
133  dp.frame.toGlobal(hhp->xerrLocal(hit), 0, yerr, ge);
134 #else
135  hhp->cpeParams()
136  .detParams(hhp->detectorIndex(hit))
137  .frame.toGlobal(hhp->xerrLocal(hit), 0, hhp->yerrLocal(hit), ge);
138 #endif
139 
140 #ifdef BL_DUMP_HITS
141  bool dump = foundNtuplets->size(tkid) == 5;
142  if (dump) {
143  printf("Track id %d %d Hit %d on %d\nGlobal: hits.col(%d) << %f,%f,%f\n",
144  local_idx,
145  tkid,
146  hit,
147  hhp->detectorIndex(hit),
148  i,
149  hhp->xGlobal(hit),
150  hhp->yGlobal(hit),
151  hhp->zGlobal(hit));
152  printf("Error: hits_ge.col(%d) << %e,%e,%e,%e,%e,%e\n", i, ge[0], ge[1], ge[2], ge[3], ge[4], ge[5]);
153  }
154 #endif
155 
156  hits.col(i) << hhp->xGlobal(hit), hhp->yGlobal(hit), hhp->zGlobal(hit);
157  hits_ge.col(i) << ge[0], ge[1], ge[2], ge[3], ge[4], ge[5];
158  }
159  brokenline::fastFit(hits, fast_fit);
160 
161  // no NaN here....
162  assert(fast_fit(0) == fast_fit(0));
163  assert(fast_fit(1) == fast_fit(1));
164  assert(fast_fit(2) == fast_fit(2));
165  assert(fast_fit(3) == fast_fit(3));
166  }
167 }
168 
169 template <int N>
170 __global__ void kernel_BLFit(caConstants::TupleMultiplicity const *__restrict__ tupleMultiplicity,
171  double bField,
173  tindex_type const *__restrict__ ptkids,
174  double *__restrict__ phits,
175  float *__restrict__ phits_ge,
176  double *__restrict__ pfast_fit) {
177  assert(results);
178  assert(pfast_fit);
179 
180  // same as above...
181 
182  // look in bin for this hit multiplicity
183  auto local_start = blockIdx.x * blockDim.x + threadIdx.x;
184  for (int local_idx = local_start, nt = riemannFit::maxNumberOfConcurrentFits; local_idx < nt;
185  local_idx += gridDim.x * blockDim.x) {
186  if (invalidTkId == ptkids[local_idx])
187  break;
188 
189  auto tkid = ptkids[local_idx];
190 
192 
193  riemannFit::Map3xNd<N> hits(phits + local_idx);
194  riemannFit::Map4d fast_fit(pfast_fit + local_idx);
195  riemannFit::Map6xNf<N> hits_ge(phits_ge + local_idx);
196 
198 
201 
202  brokenline::prepareBrokenLineData(hits, fast_fit, bField, data);
203  brokenline::lineFit(hits_ge, fast_fit, bField, data, line);
204  brokenline::circleFit(hits, hits_ge, fast_fit, bField, data, circle);
205 
206  results->stateAtBS.copyFromCircle(circle.par, circle.cov, line.par, line.cov, 1.f / float(bField), tkid);
207  results->pt(tkid) = float(bField) / float(std::abs(circle.par(2)));
208  results->eta(tkid) = asinhf(line.par(0));
209  results->chi2(tkid) = (circle.chi2 + line.chi2) / (2 * N - 5);
210 
211 #ifdef BROKENLINE_DEBUG
212  if (!(circle.chi2 >= 0) || !(line.chi2 >= 0))
213  printf("kernelBLFit failed! %f/%f\n", circle.chi2, line.chi2);
214  printf("kernelBLFit size %d for %d hits circle.par(0,1,2): %d %f,%f,%f\n",
215  N,
216  nHits,
217  tkid,
218  circle.par(0),
219  circle.par(1),
220  circle.par(2));
221  printf("kernelBLHits line.par(0,1): %d %f,%f\n", tkid, line.par(0), line.par(1));
222  printf("kernelBLHits chi2 cov %f/%f %e,%e,%e,%e,%e\n",
223  circle.chi2,
224  line.chi2,
225  circle.cov(0, 0),
226  circle.cov(1, 1),
227  circle.cov(2, 2),
228  line.cov(0, 0),
229  line.cov(1, 1));
230 #endif
231  }
232 }
const dim3 threadIdx
Definition: cudaCompat.h:29
__host__ __device__ void lineFit(const M6xN &hits_ge, const V4 &fast_fit, const double bField, const PreparedBrokenLineData< n > &data, riemannFit::LineFit &line_results)
Performs the Broken Line fit in the straight track case (that is, the fit parameters are only the int...
Definition: BrokenLine.h:463
caConstants::tindex_type tindex_type
uint16_t tindex_type
Definition: CAConstants.h:69
constexpr uint16_t pixelPitchY
caConstants::TupleMultiplicity const *__restrict__ HitsOnGPU const *__restrict__ tindex_type *__restrict__ ptkids
const dim3 gridDim
Definition: cudaCompat.h:33
caConstants::TupleMultiplicity const *__restrict__ HitsOnGPU const *__restrict__ tindex_type *__restrict__ double *__restrict__ float *__restrict__ phits_ge
constexpr uint32_t maxTuples
Definition: CAConstants.h:46
#define __global__
Definition: cudaCompat.h:19
dictionary results
const dim3 blockDim
Definition: cudaCompat.h:30
auto local_start
list status
Definition: mps_update.py:107
constexpr int kGenErrorQBins
Vector3d par
parameter: (X0,Y0,R)
Definition: FitResult.h:27
assert(be >=bs)
caConstants::TupleMultiplicity const *__restrict__ HitsOnGPU const *__restrict__ tindex_type *__restrict__ double *__restrict__ float *__restrict__ double *__restrict__ uint32_t nHitsL
eigenSoA::ScalarSoA< float, S > pt
__host__ __device__ void prepareBrokenLineData(const M3xN &hits, const V4 &fast_fit, const double bField, PreparedBrokenLineData< n > &results)
Computes the data needed for the Broken Line fit procedure that are mainly common for the circle and ...
Definition: BrokenLine.h:150
eigenSoA::ScalarSoA< float, S > chi2
auto const & foundNtuplets
Eigen::Map< Matrix3xNd< N >, 0, Eigen::Stride< 3 *stride, stride > > Map3xNd
Definition: HelixFitOnGPU.h:23
caConstants::TupleMultiplicity const *__restrict__ tupleMultiplicity
printf("params %d %f %f %f\n", minT, eps, errmax, chi2max)
TrajectoryStateSoAT< S > stateAtBS
__host__ __device__ void fastFit(const M3xN &hits, V4 &result)
A very fast helix fit.
Definition: BrokenLine.h:258
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
TrackSoA::HitContainer HitContainer
TrackSoAHeterogeneousT< maxNumber()> TrackSoA
const dim3 blockIdx
Definition: cudaCompat.h:32
Eigen::Map< Matrix6xNf< N >, 0, Eigen::Stride< 6 *stride, stride > > Map6xNf
Definition: HelixFitOnGPU.h:28
int nt
Definition: AMPTWrapper.h:42
caConstants::TupleMultiplicity const *__restrict__ HitsOnGPU const *__restrict__ tindex_type *__restrict__ double *__restrict__ float *__restrict__ double *__restrict__ pfast_fit
#define N
Definition: blowfish.cc:9
int totTK
constexpr auto invalidTkId
caConstants::TupleMultiplicity const *__restrict__ HitsOnGPU const *__restrict__ hhp
Eigen::Map< Vector4d, 0, Eigen::InnerStride< stride > > Map4d
Definition: HelixFitOnGPU.h:30
caConstants::TupleMultiplicity const CAHitNtupletGeneratorKernelsGPU::HitToTuple const cms::cuda::AtomicPairCounter GPUCACell const *__restrict__ uint32_t const *__restrict__ gpuPixelDoublets::CellNeighborsVector const gpuPixelDoublets::CellTracksVector const GPUCACell::OuterHitOfCell const int32_t nHits
void __syncthreads()
Definition: cudaCompat.h:108
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:79
caConstants::TupleMultiplicity const *__restrict__ HitsOnGPU const *__restrict__ tindex_type *__restrict__ double *__restrict__ float *__restrict__ double *__restrict__ uint32_t uint32_t nHitsH
Vector2d par
(cotan(theta),Zip)
Definition: FitResult.h:39
__host__ __device__ void circleFit(const M3xN &hits, const M6xN &hits_ge, const V4 &fast_fit, const double bField, PreparedBrokenLineData< n > &data, karimaki_circle_fit &circle_results)
Performs the Broken Line fit in the curved track case (that is, the fit parameters are the intercepti...
Definition: BrokenLine.h:314
eigenSoA::ScalarSoA< float, S > eta
data needed for the Broken Line fit procedure.
Definition: BrokenLine.h:24
T1 atomicAdd(T1 *a, T2 b)
Definition: cudaCompat.h:61
constexpr int kNumErrorBins
tuple dump
OutputFilePath = cms.string(&#39;/tmp/zhokin/&#39;), OutputFileExt = cms.string(&#39;&#39;),.
caConstants::TupleMultiplicity const *__restrict__ HitsOnGPU const *__restrict__ tindex_type *__restrict__ double *__restrict__ phits
constexpr uint32_t maxNumberOfConcurrentFits
Definition: HelixFitOnGPU.h:12