Project
Loading...
Searching...
No Matches
testDCAFitterNGPU.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
12#define BOOST_TEST_MODULE Test DCAFitterN class
13#define BOOST_TEST_MAIN
14#define BOOST_TEST_DYN_LINK
15#include <boost/test/unit_test.hpp>
16
17#include "DCAFitter/DCAFitterN.h"
18#include "CommonUtils/TreeStreamRedirector.h"
19#include <TRandom.h>
20#include <TGenPhaseSpace.h>
21#include <TLorentzVector.h>
22#include <TStopwatch.h>
23#include <Math/SVector.h>
24#include <array>
25
26namespace o2
27{
28namespace vertexing
29{
30
31using Vec3D = ROOT::Math::SVector<double, 3>;
32
33template <class FITTER>
34float checkResults(o2::utils::TreeStreamRedirector& outs, std::string& treeName, FITTER& fitter,
35 Vec3D& vgen, TLorentzVector& genPar, const std::vector<double>& dtMass)
36{
37 int nCand = fitter.getNCandidates();
38 std::array<float, 3> p;
39 float distMin = 1e9;
40 bool absDCA = fitter.getUseAbsDCA();
41 bool useWghDCA = fitter.getWeightedFinalPCA();
42 for (int ic = 0; ic < nCand; ic++) {
43 const auto& vtx = fitter.getPCACandidate(ic);
44 auto df = vgen;
45 df -= vtx;
46
47 TLorentzVector moth, prong;
48 for (int i = 0; i < fitter.getNProngs(); i++) {
49 const auto& trc = fitter.getTrack(i, ic);
50 trc.getPxPyPzGlo(p);
51 prong.SetVectM({p[0], p[1], p[2]}, dtMass[i]);
52 moth += prong;
53 }
54 auto nIter = fitter.getNIterations(ic);
55 auto chi2 = fitter.getChi2AtPCACandidate(ic);
56 double dst = TMath::Sqrt(df[0] * df[0] + df[1] * df[1] + df[2] * df[2]);
57 distMin = dst < distMin ? dst : distMin;
58 auto parentTrack = fitter.createParentTrackParCov(ic);
59 outs << treeName.c_str() << "cand=" << ic << "ncand=" << nCand << "nIter=" << nIter << "chi2=" << chi2
60 << "genPart=" << genPar << "recPart=" << moth
61 << "genX=" << vgen[0] << "genY=" << vgen[1] << "genZ=" << vgen[2]
62 << "dx=" << df[0] << "dy=" << df[1] << "dz=" << df[2] << "dst=" << dst
63 << "useAbsDCA=" << absDCA << "useWghDCA=" << useWghDCA << "parent=" << parentTrack;
64 for (int i = 0; i < fitter.getNProngs(); i++) {
65 outs << treeName.c_str() << fmt::format("prong{}=", i).c_str() << fitter.getTrack(i, ic);
66 }
67 outs << treeName.c_str() << "\n";
68 }
69 return distMin;
70}
71
72TLorentzVector generate(Vec3D& vtx, std::vector<o2::track::TrackParCov>& vctr, float bz,
73 TGenPhaseSpace& genPHS, double parMass, const std::vector<double>& dtMass, std::vector<int> forceQ)
74{
75 const float errYZ = 1e-2, errSlp = 1e-3, errQPT = 2e-2;
76 std::array<float, 15> covm = {
77 errYZ * errYZ,
78 0., errYZ * errYZ,
79 0, 0., errSlp * errSlp,
80 0., 0., 0., errSlp * errSlp,
81 0., 0., 0., 0., errQPT * errQPT};
82 bool accept = true;
83 TLorentzVector parent, d0, d1, d2;
84 do {
85 accept = true;
86 double y = gRandom->Rndm() - 0.5;
87 double pt = 0.1 + gRandom->Rndm() * 3;
88 double mt = TMath::Sqrt(parMass * parMass + pt * pt);
89 double pz = mt * TMath::SinH(y);
90 double phi = gRandom->Rndm() * TMath::Pi() * 2;
91 double en = mt * TMath::CosH(y);
92 double rdec = 10.; // radius of the decay
93 vtx[0] = rdec * TMath::Cos(phi);
94 vtx[1] = rdec * TMath::Sin(phi);
95 vtx[2] = rdec * pz / pt;
96 parent.SetPxPyPzE(pt * TMath::Cos(phi), pt * TMath::Sin(phi), pz, en);
97 int nd = dtMass.size();
98 genPHS.SetDecay(parent, nd, dtMass.data());
99 genPHS.Generate();
100 vctr.clear();
101 float p[4];
102 for (int i = 0; i < nd; i++) {
103 auto* dt = genPHS.GetDecay(i);
104 if (dt->Pt() < 0.05) {
105 accept = false;
106 break;
107 }
108 dt->GetXYZT(p);
109 float s, c, x;
110 std::array<float, 5> params;
111 o2::math_utils::sincos(dt->Phi(), s, c);
112 o2::math_utils::rotateZInv(vtx[0], vtx[1], x, params[0], s, c);
113
114 params[1] = vtx[2];
115 params[2] = 0.; // since alpha = phi
116 params[3] = 1. / TMath::Tan(dt->Theta());
117 params[4] = (i % 2 ? -1. : 1.) / dt->Pt();
118 covm[14] = errQPT * errQPT * params[4] * params[4];
119 //
120 // randomize
121 float r1, r2;
122 gRandom->Rannor(r1, r2);
123 params[0] += r1 * errYZ;
124 params[1] += r2 * errYZ;
125 gRandom->Rannor(r1, r2);
126 params[2] += r1 * errSlp;
127 params[3] += r2 * errSlp;
128 params[4] *= gRandom->Gaus(1., errQPT);
129 if (forceQ[i] == 0) {
130 params[4] = 0.; // impose straight track
131 }
132 auto& trc = vctr.emplace_back(x, dt->Phi(), params, covm);
133 float rad = forceQ[i] == 0 ? 600. : TMath::Abs(1. / trc.getCurvature(bz));
134 if (!trc.propagateTo(trc.getX() + (gRandom->Rndm() - 0.5) * rad * 0.05, bz) ||
135 !trc.rotate(trc.getAlpha() + (gRandom->Rndm() - 0.5) * 0.2)) {
136 printf("Failed to randomize ");
137 trc.print();
138 }
139 }
140 } while (!accept);
141
142 return parent;
143}
144
145#ifdef DO_SINGLE_THREAD_TEST
146BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
147{
148 // gRandom->Delete();
149 // gRandom = new TRandom(42);
150 o2::utils::TreeStreamRedirector outStream("dcafitterNTest.root");
151
152 TGenPhaseSpace genPHS;
153 constexpr double ele = 0.00051;
154 constexpr double gamma = 2 * ele + 1e-6;
155 constexpr double pion = 0.13957;
156 constexpr double k0 = 0.49761;
157 constexpr double kch = 0.49368;
158 constexpr double dch = 1.86965;
159 std::vector<double> gammadec = {ele, ele};
160 std::vector<double> k0dec = {pion, pion};
161 std::vector<double> dchdec = {pion, kch, pion};
162 std::vector<o2::track::TrackParCov> vctracks;
163 Vec3D vtxGen;
164
165 double bz = 5.0;
166 // 2 prongs vertices
167 {
168 LOG(info) << "Processing 2-prong Helix - Helix case";
169 std::vector<int> forceQ{1, 1};
170
171 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
172 ft.setBz(bz);
173 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
174 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
175 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
176 ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway
177 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
178 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
179
180 std::string treeName2A = "pr2a", treeName2AW = "pr2aw", treeName2W = "pr2w";
181 TStopwatch swA, swAW, swW;
182 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
183 double meanDA = 0, meanDAW = 0, meanDW = 0;
184 swA.Stop();
185 swAW.Stop();
186 swW.Stop();
187 for (int iev = 0; iev < NTest; iev++) {
188 auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ);
189
190 ft.setUseAbsDCA(true);
191 swA.Start(false);
192 int ncA = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
193 swA.Stop();
194 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
195 if (ncA) {
196 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec);
197 meanDA += minD;
198 nfoundA++;
199 }
200
201 ft.setUseAbsDCA(true);
202 ft.setWeightedFinalPCA(true);
203 swAW.Start(false);
204 int ncAW = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
205 swAW.Stop();
206 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
207 if (ncAW) {
208 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec);
209 meanDAW += minD;
210 nfoundAW++;
211 }
212
213 ft.setUseAbsDCA(false);
214 ft.setWeightedFinalPCA(false);
215 swW.Start(false);
216 int ncW = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
217 swW.Stop();
218 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
219 if (ncW) {
220 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec);
221 meanDW += minD;
222 nfoundW++;
223 }
224 }
225
226 meanDA /= nfoundA ? nfoundA : 1;
227 meanDAW /= nfoundAW ? nfoundA : 1;
228 meanDW /= nfoundW ? nfoundW : 1;
229 LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix";
230 LOG(info) << "2-prongs with abs.dist minimization: eff = " << float(nfoundA) / NTest
231 << " mean.dist to truth: " << meanDA << " Total time: " << swA.CpuTime() * 1000 << " ms";
232 LOG(info) << "2-prongs with abs.dist but wghPCA: eff = " << float(nfoundAW) / NTest
233 << " mean.dist to truth: " << meanDAW << " Total time: " << swAW.CpuTime() * 1000 << " ms";
234 LOG(info) << "2-prongs with wgh.dist minimization: eff = " << float(nfoundW) / NTest
235 << " mean.dist to truth: " << meanDW << " Total time: " << swW.CpuTime() * 1000 << " ms";
236 BOOST_CHECK(nfoundA > 0.99 * NTest);
237 BOOST_CHECK(nfoundAW > 0.99 * NTest);
238 BOOST_CHECK(nfoundW > 0.99 * NTest);
239 BOOST_CHECK(meanDA < 0.1);
240 BOOST_CHECK(meanDAW < 0.1);
241 BOOST_CHECK(meanDW < 0.1);
242 }
243
244 // 2 prongs vertices with collinear tracks (gamma conversion)
245 {
246 LOG(info) << "Processing 2-prong Helix - Helix case gamma conversion";
247 std::vector<int> forceQ{1, 1};
248
249 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
250 ft.setBz(bz);
251 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
252 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
253 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
254 ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway
255 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
256 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
257
258 std::string treeName2A = "gpr2a", treeName2AW = "gpr2aw", treeName2W = "gpr2w";
259 TStopwatch swA, swAW, swW;
260 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
261 double meanDA = 0, meanDAW = 0, meanDW = 0;
262 swA.Stop();
263 swAW.Stop();
264 swW.Stop();
265 for (int iev = 0; iev < NTest; iev++) {
266 auto genParent = generate(vtxGen, vctracks, bz, genPHS, gamma, gammadec, forceQ);
267
268 ft.setUseAbsDCA(true);
269 swA.Start(false);
270 int ncA = device::process(1, 1, ft, vctracks[0], vctracks[1]);
271 swA.Stop();
272 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
273 if (ncA) {
274 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, gammadec);
275 meanDA += minD;
276 nfoundA++;
277 }
278
279 ft.setUseAbsDCA(true);
280 ft.setWeightedFinalPCA(true);
281 swAW.Start(false);
282 int ncAW = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
283 swAW.Stop();
284 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
285 if (ncAW) {
286 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, gammadec);
287 meanDAW += minD;
288 nfoundAW++;
289 }
290
291 ft.setUseAbsDCA(false);
292 ft.setWeightedFinalPCA(false);
293 swW.Start(false);
294 int ncW = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
295 swW.Stop();
296 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
297 if (ncW) {
298 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, gammadec);
299 meanDW += minD;
300 nfoundW++;
301 }
302 }
303
304 meanDA /= nfoundA ? nfoundA : 1;
305 meanDAW /= nfoundA ? nfoundA : 1;
306 meanDW /= nfoundW ? nfoundW : 1;
307 LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix from gamma conversion";
308 LOG(info) << "2-prongs with abs.dist minimization: eff = " << float(nfoundA) / NTest
309 << " mean.dist to truth: " << meanDA << " Total time: " << swA.CpuTime();
310 LOG(info) << "2-prongs with abs.dist but wghPCA: eff = " << float(nfoundAW) / NTest
311 << " mean.dist to truth: " << meanDAW << " Total time: " << swAW.CpuTime();
312 LOG(info) << "2-prongs with wgh.dist minimization: eff = " << float(nfoundW) / NTest
313 << " mean.dist to truth: " << meanDW << " Total time: " << swW.CpuTime();
314 BOOST_CHECK(nfoundA > 0.99 * NTest);
315 BOOST_CHECK(nfoundAW > 0.99 * NTest);
316 BOOST_CHECK(nfoundW > 0.99 * NTest);
317 BOOST_CHECK(meanDA < 2.1);
318 BOOST_CHECK(meanDAW < 2.1);
319 BOOST_CHECK(meanDW < 2.1);
320 }
321
322 // 2 prongs vertices with one of charges set to 0: Helix : Line
323 {
324 std::vector<int> forceQ{1, 1};
325 LOG(info) << "Processing 2-prong Helix - Line case";
326 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
327 ft.setBz(bz);
328 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
329 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
330 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
331 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
332 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
333
334 std::string treeName2A = "pr2aHL", treeName2AW = "pr2awHL", treeName2W = "pr2wHL";
335 TStopwatch swA, swAW, swW;
336 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
337 double meanDA = 0, meanDAW = 0, meanDW = 0;
338 swA.Stop();
339 swAW.Stop();
340 swW.Stop();
341 for (int iev = 0; iev < NTest; iev++) {
342 forceQ[iev % 2] = 1;
343 forceQ[1 - iev % 2] = 0;
344 auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ);
345
346 ft.setUseAbsDCA(true);
347 swA.Start(false);
348 int ncA = device::process(1, 1, ft, vctracks[0], vctracks[1]);
349 swA.Stop();
350 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
351 if (ncA) {
352 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec);
353 meanDA += minD;
354 nfoundA++;
355 }
356
357 ft.setUseAbsDCA(true);
358 ft.setWeightedFinalPCA(true);
359 swAW.Start(false);
360 int ncAW = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
361 swAW.Stop();
362 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
363 if (ncAW) {
364 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec);
365 meanDAW += minD;
366 nfoundAW++;
367 }
368
369 ft.setUseAbsDCA(false);
370 ft.setWeightedFinalPCA(false);
371 swW.Start(false);
372 int ncW = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
373 swW.Stop();
374 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
375 if (ncW) {
376 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec);
377 meanDW += minD;
378 nfoundW++;
379 }
380 }
381
382 meanDA /= nfoundA ? nfoundA : 1;
383 meanDAW /= nfoundAW ? nfoundAW : 1;
384 meanDW /= nfoundW ? nfoundW : 1;
385 LOG(info) << "Processed " << NTest << " 2-prong vertices: Helix : Line";
386 LOG(info) << "2-prongs with abs.dist minimization: eff = " << float(nfoundA) / NTest
387 << " mean.dist to truth: " << meanDA << " Total time: " << swA.CpuTime();
388 LOG(info) << "2-prongs with abs.dist but wghPCA: eff = " << float(nfoundAW) / NTest
389 << " mean.dist to truth: " << meanDAW << " Total time: " << swAW.CpuTime();
390 LOG(info) << "2-prongs with wgh.dist minimization: eff = " << float(nfoundW) / NTest
391 << " mean.dist to truth: " << meanDW << " Total time: " << swW.CpuTime();
392 BOOST_CHECK(nfoundA > 0.99 * NTest);
393 BOOST_CHECK(nfoundAW > 0.99 * NTest);
394 BOOST_CHECK(nfoundW > 0.99 * NTest);
395 BOOST_CHECK(meanDA < 0.1);
396 BOOST_CHECK(meanDAW < 0.1);
397 BOOST_CHECK(meanDW < 0.1);
398 }
399
400 // 2 prongs vertices with both of charges set to 0: Line : Line
401 {
402 std::vector<int> forceQ{0, 0};
403 LOG(info) << "Processing 2-prong Line - Line case";
404 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
405 ft.setBz(bz);
406 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
407 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
408 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
409 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
410 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
411
412 std::string treeName2A = "pr2aLL", treeName2AW = "pr2awLL", treeName2W = "pr2wLL";
413 TStopwatch swA, swAW, swW;
414 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
415 double meanDA = 0, meanDAW = 0, meanDW = 0;
416 swA.Stop();
417 swAW.Stop();
418 swW.Stop();
419 for (int iev = 0; iev < NTest; iev++) {
420 forceQ[0] = forceQ[1] = 0;
421 auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ);
422
423 ft.setUseAbsDCA(true);
424 swA.Start(false);
425 int ncA = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
426 swA.Stop();
427 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
428 if (ncA) {
429 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec);
430 meanDA += minD;
431 nfoundA++;
432 }
433
434 ft.setUseAbsDCA(true);
435 ft.setWeightedFinalPCA(true);
436 swAW.Start(false);
437 int ncAW = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
438 swAW.Stop();
439 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
440 if (ncAW) {
441 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec);
442 meanDAW += minD;
443 nfoundAW++;
444 }
445
446 ft.setUseAbsDCA(false);
447 ft.setWeightedFinalPCA(false);
448 swW.Start(false);
449 int ncW = device::process(1, 1, ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
450 swW.Stop();
451 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
452 if (ncW) {
453 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec);
454 meanDW += minD;
455 nfoundW++;
456 }
457 }
458
459 meanDA /= nfoundA ? nfoundA : 1;
460 meanDAW /= nfoundAW ? nfoundAW : 1;
461 meanDW /= nfoundW ? nfoundW : 1;
462 LOG(info) << "Processed " << NTest << " 2-prong vertices: Line : Line";
463 LOG(info) << "2-prongs with abs.dist minimization: eff = " << float(nfoundA) / NTest
464 << " mean.dist to truth: " << meanDA << " Total time: " << swA.CpuTime();
465 LOG(info) << "2-prongs with abs.dist but wghPCA: eff = " << float(nfoundAW) / NTest
466 << " mean.dist to truth: " << meanDAW << " Total time: " << swAW.CpuTime();
467 LOG(info) << "2-prongs with wgh.dist minimization: eff = " << float(nfoundW) / NTest
468 << " mean.dist to truth: " << meanDW << " Total time: " << swW.CpuTime();
469 BOOST_CHECK(nfoundA > 0.99 * NTest);
470 BOOST_CHECK(nfoundAW > 0.99 * NTest);
471 BOOST_CHECK(nfoundW > 0.99 * NTest);
472 BOOST_CHECK(meanDA < 0.1);
473 BOOST_CHECK(meanDAW < 0.1);
474 BOOST_CHECK(meanDW < 0.1);
475 }
476
477 // 3 prongs vertices
478 {
479 std::vector<int> forceQ{1, 1, 1};
480
481 o2::vertexing::DCAFitterN<3> ft; // 3 prong fitter
482 ft.setBz(bz);
483 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
484 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
485 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
486 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
487 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
488
489 std::string treeName3A = "pr3a", treeName3AW = "pr3aw", treeName3W = "pr3w";
490 TStopwatch swA, swAW, swW;
491 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
492 double meanDA = 0, meanDAW = 0, meanDW = 0;
493 swA.Stop();
494 swAW.Stop();
495 swW.Stop();
496 for (int iev = 0; iev < NTest; iev++) {
497 auto genParent = generate(vtxGen, vctracks, bz, genPHS, dch, dchdec, forceQ);
498
499 ft.setUseAbsDCA(true);
500 swA.Start(false);
501 int ncA = device::process(1, 1, ft, vctracks[0], vctracks[1], vctracks[2]);
502 swA.Stop();
503 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
504 if (ncA) {
505 auto minD = checkResults(outStream, treeName3A, ft, vtxGen, genParent, dchdec);
506 meanDA += minD;
507 nfoundA++;
508 }
509
510 ft.setUseAbsDCA(true);
511 ft.setWeightedFinalPCA(true);
512 swAW.Start(false);
513 int ncAW = device::process(1, 1, ft, vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
514 swAW.Stop();
515 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
516 if (ncAW) {
517 auto minD = checkResults(outStream, treeName3AW, ft, vtxGen, genParent, dchdec);
518 meanDAW += minD;
519 nfoundAW++;
520 }
521
522 ft.setUseAbsDCA(false);
523 ft.setWeightedFinalPCA(false);
524 swW.Start(false);
525 int ncW = device::process(1, 1, ft, vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
526 swW.Stop();
527 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
528 if (ncW) {
529 auto minD = checkResults(outStream, treeName3W, ft, vtxGen, genParent, dchdec);
530 meanDW += minD;
531 nfoundW++;
532 }
533 }
534
535 meanDA /= nfoundA ? nfoundA : 1;
536 meanDAW /= nfoundAW ? nfoundAW : 1;
537 meanDW /= nfoundW ? nfoundW : 1;
538 LOG(info) << "Processed " << NTest << " 3-prong vertices";
539 LOG(info) << "3-prongs with abs.dist minimization: eff = " << float(nfoundA) / NTest
540 << " mean.dist to truth: " << meanDA << " Total time: " << swA.CpuTime();
541 LOG(info) << "3-prongs with abs.dist but wghPCA: eff = " << float(nfoundAW) / NTest
542 << " mean.dist to truth: " << meanDAW << " Total time: " << swAW.CpuTime();
543 LOG(info) << "3-prongs with wgh.dist minimization: eff = " << float(nfoundW) / NTest
544 << " mean.dist to truth: " << meanDW << " Total time: " << swW.CpuTime();
545 BOOST_CHECK(nfoundA > 0.99 * NTest);
546 BOOST_CHECK(nfoundAW > 0.99 * NTest);
547 BOOST_CHECK(nfoundW > 0.99 * NTest);
548 BOOST_CHECK(meanDA < 0.1);
549 BOOST_CHECK(meanDAW < 0.1);
550 BOOST_CHECK(meanDW < 0.1);
551 }
552
553 outStream.Close();
554}
555#endif
556
557BOOST_AUTO_TEST_CASE(DCAFitterNProngsBulk)
558{
559 const char* nThreadsEnvVarName = "DCAFITTERGPU_TEST_NTHREADS";
560 const char* nBlocksEnvVarName = "DCAFITTERGPU_TEST_NBLOCKS";
561 const char* nBatchesEnvVarName = "DCAFITTERGPU_TEST_NBATCHES";
562 const char* nTestsEnvVarName = "DCAFITTERGPU_TEST_NTESTS";
563 int nBlocks = std::getenv(nBlocksEnvVarName) == nullptr ? 30 : std::stoi(std::getenv(nBlocksEnvVarName));
564 int nThreads = std::getenv(nThreadsEnvVarName) == nullptr ? 256 : std::stoi(std::getenv(nThreadsEnvVarName));
565 int nBatches = std::getenv(nBatchesEnvVarName) == nullptr ? 8 : std::stoi(std::getenv(nBatchesEnvVarName));
566 int NTest = std::getenv(nTestsEnvVarName) == nullptr ? 100001 : std::stoi(std::getenv(nTestsEnvVarName));
567
568 o2::utils::TreeStreamRedirector outStreamB("dcafitterNTestBulk.root");
569
570 TGenPhaseSpace genPHS;
571 constexpr double ele = 0.00051;
572 constexpr double gamma = 2 * ele + 1e-6;
573 constexpr double pion = 0.13957;
574 constexpr double k0 = 0.49761;
575 constexpr double kch = 0.49368;
576 constexpr double dch = 1.86965;
577 std::vector<double> gammadec = {ele, ele};
578 std::vector<double> k0dec = {pion, pion};
579 std::vector<double> dchdec = {pion, kch, pion};
580 std::vector<std::vector<o2::track::TrackParCov>> vctracks(3, std::vector<o2::track::TrackParCov>(NTest));
581 std::vector<Vec3D> vtxGen(NTest);
582
583 double bz = 5.0;
584 { // 2 prongs vertices bulk processing
585 LOG(info) << "\n\nBulk-processing 2-prong Helix - Helix case";
586 std::vector<int> forceQ{1, 1};
587
588 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
589 ft.setBz(bz);
590 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
591 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
592 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
593 ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway
594 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
595 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
596
597 std::vector<o2::vertexing::DCAFitterN<2>> fitters_host(NTest);
598 std::vector<TLorentzVector> genParents(NTest);
599
600 std::string treeName2Abulk = "pr2aBulk", treeName2AWbulk = "pr2awBulk", treeName2Wbulk = "pr2wBulk";
601 TStopwatch swAb, swAWb, swWb;
602 int nfoundAb = 0, nfoundAWb = 0, nfoundWb = 0;
603 double meanDAb = 0, meanDAWb = 0, meanDWb = 0;
604 swAb.Stop();
605 swAWb.Stop();
606 swWb.Stop();
607
608 ft.setUseAbsDCA(true);
609 std::fill(fitters_host.begin(), fitters_host.end(), ft);
610 for (int iev = 0; iev < NTest; iev++) {
611 std::vector<o2::track::TrackParCov> vc(2);
612 genParents[iev] = generate(vtxGen[iev], vc, bz, genPHS, k0, k0dec, forceQ);
613 vctracks[0][iev] = vc[0];
614 vctracks[1][iev] = vc[1];
615 }
616
617 swAb.Start(false);
618 std::vector<int> ncAb(NTest, 0);
619 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
620 swAb.Stop();
621
622 for (int iev = 0; iev < NTest; iev++) {
623 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAb[iev] << " Chi2: " << (ncAb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
624 if (ncAb[iev]) {
625 auto minDb = checkResults(outStreamB, treeName2Abulk, fitters_host[iev], vtxGen[iev], genParents[iev], k0dec);
626 meanDAb += minDb;
627 nfoundAb++;
628 }
629 }
630
631 ft.setUseAbsDCA(true);
632 ft.setWeightedFinalPCA(true);
633 std::fill(fitters_host.begin(), fitters_host.end(), ft);
634 swAWb.Start(false);
635 std::vector<int> ncAWb(NTest, 0);
636 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAWb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
637 swAWb.Stop();
638
639 for (int iev = 0; iev < NTest; iev++) {
640 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAWb[iev] << " Chi2: " << (ncAWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
641 if (ncAWb[iev]) {
642 auto minDb = checkResults(outStreamB, treeName2AWbulk, fitters_host[iev], vtxGen[iev], genParents[iev], k0dec);
643 meanDAWb += minDb;
644 nfoundAWb++;
645 }
646 }
647
648 ft.setUseAbsDCA(false);
649 ft.setWeightedFinalPCA(false);
650 std::fill(fitters_host.begin(), fitters_host.end(), ft);
651 swWb.Start(false);
652 std::vector<int> ncWb(NTest, 0);
653 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncWb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
654 swWb.Stop();
655
656 for (int iev = 0; iev < NTest; iev++) {
657 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncWb[iev] << " Chi2: " << (ncWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
658 if (ncWb[iev]) {
659 auto minDb = checkResults(outStreamB, treeName2Wbulk, fitters_host[iev], vtxGen[iev], genParents[iev], k0dec);
660 meanDWb += minDb;
661 nfoundWb++;
662 }
663 }
664 //
665 meanDAb /= nfoundAb ? nfoundAb : 1;
666 meanDAWb /= nfoundAWb ? nfoundAWb : 1;
667 meanDWb /= nfoundWb ? nfoundWb : 1;
668 LOGP(info, "Bulk-processed {} 2-prong vertices Helix : Helix", NTest);
669 LOG(info) << "2-prongs with abs.dist minimization: eff = " << float(nfoundAb) / NTest
670 << " mean.dist to truth: " << meanDAb << " Total time: " << swAb.CpuTime() * 1000 << " ms";
671 LOG(info) << "2-prongs with abs.dist but wghPCA: eff = " << float(nfoundAWb) / NTest
672 << " mean.dist to truth: " << meanDAWb << " Total time: " << swAWb.CpuTime() * 1000 << " ms";
673 LOG(info) << "2-prongs with wgh.dist minimization: eff = " << float(nfoundWb) / NTest
674 << " mean.dist to truth: " << meanDWb << " Total time: " << swWb.CpuTime() * 1000 << " ms";
675 BOOST_CHECK(nfoundAb > 0.99 * NTest);
676 BOOST_CHECK(nfoundAWb > 0.99 * NTest);
677 BOOST_CHECK(nfoundWb > 0.99 * NTest);
678 BOOST_CHECK(meanDAb < 0.1);
679 BOOST_CHECK(meanDAWb < 0.1);
680 BOOST_CHECK(meanDWb < 0.1);
681 }
682
683 { // 2 prongs vertices bulk processing for gamma conversion
684 LOG(info) << "\n\nBulk-processing 2-prong Helix - Helix case gamma conversion";
685 std::vector<int> forceQ{1, 1};
686
687 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
688 ft.setBz(bz);
689 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
690 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
691 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
692 ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway
693 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
694 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
695
696 std::vector<o2::vertexing::DCAFitterN<2>> fitters_host(NTest);
697 std::vector<TLorentzVector> genParents(NTest);
698
699 std::string treeName2Abulk = "gpr2aBulk", treeName2AWbulk = "gpr2awBulk", treeName2Wbulk = "gpr2wBulk";
700 TStopwatch swAb, swAWb, swWb;
701 int nfoundAb = 0, nfoundAWb = 0, nfoundWb = 0;
702 double meanDAb = 0, meanDAWb = 0, meanDWb = 0;
703 swAb.Stop();
704 swAWb.Stop();
705 swWb.Stop();
706
707 ft.setUseAbsDCA(true);
708 std::fill(fitters_host.begin(), fitters_host.end(), ft);
709 for (int iev = 0; iev < NTest; iev++) {
710 std::vector<o2::track::TrackParCov> vc(2);
711 genParents[iev] = generate(vtxGen[iev], vc, bz, genPHS, gamma, gammadec, forceQ);
712 vctracks[0][iev] = vc[0];
713 vctracks[1][iev] = vc[1];
714 }
715
716 swAb.Start(false);
717 std::vector<int> ncAb(NTest, 0);
718 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
719 swAb.Stop();
720
721 for (int iev = 0; iev < NTest; iev++) {
722 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAb[iev] << " Chi2: " << (ncAb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
723 if (ncAb[iev]) {
724 auto minDb = checkResults(outStreamB, treeName2Abulk, fitters_host[iev], vtxGen[iev], genParents[iev], gammadec);
725 meanDAb += minDb;
726 nfoundAb++;
727 }
728 }
729 //
730 ft.setUseAbsDCA(true);
731 ft.setWeightedFinalPCA(true);
732 std::fill(fitters_host.begin(), fitters_host.end(), ft);
733 swAWb.Start(false);
734 std::vector<int> ncAWb(NTest, 0);
735 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAWb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
736 swAWb.Stop();
737
738 for (int iev = 0; iev < NTest; iev++) {
739 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAWb[iev] << " Chi2: " << (ncAWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
740 if (ncAWb[iev]) {
741 auto minDb = checkResults(outStreamB, treeName2AWbulk, fitters_host[iev], vtxGen[iev], genParents[iev], gammadec);
742 meanDAWb += minDb;
743 nfoundAWb++;
744 }
745 }
746
747 ft.setUseAbsDCA(false);
748 ft.setWeightedFinalPCA(false);
749 std::fill(fitters_host.begin(), fitters_host.end(), ft);
750 swWb.Start(false);
751 std::vector<int> ncWb(NTest, 0);
752 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncWb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
753 swWb.Stop();
754
755 for (int iev = 0; iev < NTest; iev++) {
756 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncWb[iev] << " Chi2: " << (ncWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
757 if (ncWb[iev]) {
758 auto minDb = checkResults(outStreamB, treeName2Wbulk, fitters_host[iev], vtxGen[iev], genParents[iev], gammadec);
759 meanDWb += minDb;
760 nfoundWb++;
761 }
762 }
763 //
764
765 meanDAb /= nfoundAb ? nfoundAb : 1;
766 meanDAWb /= nfoundAWb ? nfoundAWb : 1;
767 meanDWb /= nfoundWb ? nfoundWb : 1;
768 LOGP(info, "Bulk-processed {} 2-prong vertices Helix : Helix from gamma conversion", NTest);
769 LOG(info) << "2-prongs with abs.dist minimization: eff = " << float(nfoundAb) / NTest
770 << " mean.dist to truth: " << meanDAb << " Total time: " << swAb.CpuTime() * 1000 << " ms";
771 LOG(info) << "2-prongs with abs.dist but wghPCA: eff = " << float(nfoundAWb) / NTest
772 << " mean.dist to truth: " << meanDAWb << " Total time: " << swAWb.CpuTime() * 1000 << " ms";
773 LOG(info) << "2-prongs with wgh.dist minimization: eff = " << float(nfoundWb) / NTest
774 << " mean.dist to truth: " << meanDWb << " Total time: " << swWb.CpuTime() * 1000 << " ms";
775 BOOST_CHECK(nfoundAb > 0.99 * NTest);
776 BOOST_CHECK(nfoundAWb > 0.99 * NTest);
777 BOOST_CHECK(nfoundWb > 0.99 * NTest);
778 BOOST_CHECK(meanDAb < 2.1);
779 BOOST_CHECK(meanDAWb < 2.1);
780 BOOST_CHECK(meanDWb < 2.1);
781 }
782
783 // 2 prongs vertices bulk processing with one of charges set to 0: Helix : Line
784 {
785 std::vector<int> forceQ{1, 1};
786 LOG(info) << "\n\nBulk-processing 2-prong Helix - Line case";
787 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
788 ft.setBz(bz);
789 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
790 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
791 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
792 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
793 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
794
795 std::vector<o2::vertexing::DCAFitterN<2>> fitters_host(NTest);
796 std::vector<TLorentzVector> genParents(NTest);
797
798 std::string treeName2Abulk = "pr2aHLb", treeName2AWbulk = "pr2awHLb", treeName2Wbulk = "pr2wHLb";
799 TStopwatch swAb, swAWb, swWb;
800 int nfoundAb = 0, nfoundAWb = 0, nfoundWb = 0;
801 double meanDAb = 0, meanDAWb = 0, meanDWb = 0;
802 swAb.Stop();
803 swAWb.Stop();
804 swWb.Stop();
805
806 for (int iev = 0; iev < NTest; iev++) {
807 forceQ[iev % 2] = 1;
808 forceQ[1 - iev % 2] = 0;
809 std::vector<o2::track::TrackParCov> vc(2);
810 genParents[iev] = generate(vtxGen[iev], vc, bz, genPHS, k0, k0dec, forceQ);
811 vctracks[0][iev] = vc[0];
812 vctracks[1][iev] = vc[1];
813 }
814 ft.setUseAbsDCA(true);
815 std::fill(fitters_host.begin(), fitters_host.end(), ft);
816
817 swAb.Start(false);
818 std::vector<int> ncAb(NTest, 0);
819 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
820 swAb.Stop();
821
822 for (int iev = 0; iev < NTest; iev++) {
823 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAb[iev] << " Chi2: " << (ncAb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
824 if (ncAb[iev]) {
825 auto minDb = checkResults(outStreamB, treeName2Abulk, fitters_host[iev], vtxGen[iev], genParents[iev], k0dec);
826 meanDAb += minDb;
827 nfoundAb++;
828 }
829 }
830
831 ft.setUseAbsDCA(true);
832 ft.setWeightedFinalPCA(true);
833 std::fill(fitters_host.begin(), fitters_host.end(), ft);
834 swAWb.Start(false);
835 std::vector<int> ncAWb(NTest, 0);
836 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAWb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
837 swAWb.Stop();
838
839 for (int iev = 0; iev < NTest; iev++) {
840 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAWb[iev] << " Chi2: " << (ncAWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
841 if (ncAWb[iev]) {
842 auto minDb = checkResults(outStreamB, treeName2AWbulk, fitters_host[iev], vtxGen[iev], genParents[iev], k0dec);
843 meanDAWb += minDb;
844 nfoundAWb++;
845 }
846 }
847
848 ft.setUseAbsDCA(false);
849 ft.setWeightedFinalPCA(false);
850 std::fill(fitters_host.begin(), fitters_host.end(), ft);
851 swWb.Start(false);
852 std::vector<int> ncWb(NTest, 0);
853 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncWb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
854 swWb.Stop();
855
856 for (int iev = 0; iev < NTest; iev++) {
857 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncWb[iev] << " Chi2: " << (ncWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
858 if (ncWb[iev]) {
859 auto minDb = checkResults(outStreamB, treeName2Wbulk, fitters_host[iev], vtxGen[iev], genParents[iev], k0dec);
860 meanDWb += minDb;
861 nfoundWb++;
862 }
863 }
864
865 //
866 meanDAb /= nfoundAb ? nfoundAb : 1;
867 meanDAWb /= nfoundAWb ? nfoundAWb : 1;
868 meanDWb /= nfoundWb ? nfoundWb : 1;
869 LOG(info) << "Bulk-processed " << NTest << " 2-prong vertices: Helix : Line";
870 LOG(info) << "2-prongs with abs.dist minimization: eff = " << float(nfoundAb) / NTest
871 << " mean.dist to truth: " << meanDAb << " Total time: " << swAb.CpuTime() * 1000 << " ms";
872 LOG(info) << "2-prongs with abs.dist but wghPCA: eff = " << float(nfoundAWb) / NTest
873 << " mean.dist to truth: " << meanDAWb << " Total time: " << swAWb.CpuTime() * 1000 << " ms";
874 LOG(info) << "2-prongs with wgh.dist minimization: eff = " << float(nfoundWb) / NTest
875 << " mean.dist to truth: " << meanDWb << " Total time: " << swWb.CpuTime() * 1000 << " ms";
876 BOOST_CHECK(nfoundAb > 0.99 * NTest);
877 BOOST_CHECK(nfoundAWb > 0.99 * NTest);
878 BOOST_CHECK(nfoundWb > 0.99 * NTest);
879 BOOST_CHECK(meanDAb < 0.1);
880 BOOST_CHECK(meanDAWb < 0.1);
881 BOOST_CHECK(meanDWb < 0.1);
882 }
883
884 // 2 prongs vertices with both of charges set to 0: Line : Line
885 {
886 std::vector<int> forceQ{0, 0};
887 LOG(info) << "\n\nBulk-processing 2-prong Line - Line case";
888 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
889 ft.setBz(bz);
890 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
891 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
892 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
893 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
894 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
895
896 std::vector<o2::vertexing::DCAFitterN<2>> fitters_host(NTest);
897 std::vector<TLorentzVector> genParents(NTest);
898
899 std::string treeName2Abulk = "pr2aLL", treeName2AWbulk = "pr2awLL", treeName2Wbulk = "pr2wLL";
900 TStopwatch swAb, swAWb, swWb;
901 int nfoundAb = 0, nfoundAWb = 0, nfoundWb = 0;
902 double meanDAb = 0, meanDAWb = 0, meanDWb = 0;
903 swAb.Stop();
904 swAWb.Stop();
905 swWb.Stop();
906 for (int iev = 0; iev < NTest; iev++) {
907 forceQ[0] = forceQ[1] = 0;
908 std::vector<o2::track::TrackParCov> vc(2);
909 genParents[iev] = generate(vtxGen[iev], vc, bz, genPHS, k0, k0dec, forceQ);
910 vctracks[0][iev] = vc[0];
911 vctracks[1][iev] = vc[1];
912 }
913
914 ft.setUseAbsDCA(true);
915 std::fill(fitters_host.begin(), fitters_host.end(), ft);
916
917 swAb.Start(false);
918 std::vector<int> ncAb(NTest, 0);
919 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
920 swAb.Stop();
921
922 for (int iev = 0; iev < NTest; iev++) {
923 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAb[iev] << " Chi2: " << (ncAb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
924 if (ncAb[iev]) {
925 auto minDb = checkResults(outStreamB, treeName2Abulk, fitters_host[iev], vtxGen[iev], genParents[iev], k0dec);
926 meanDAb += minDb;
927 nfoundAb++;
928 }
929 }
930
931 ft.setUseAbsDCA(true);
932 ft.setWeightedFinalPCA(true);
933 std::fill(fitters_host.begin(), fitters_host.end(), ft);
934 swAWb.Start(false);
935 std::vector<int> ncAWb(NTest, 0);
936 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAWb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
937 swAWb.Stop();
938 for (int iev = 0; iev < NTest; iev++) {
939 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAWb[iev] << " Chi2: " << (ncAWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
940 if (ncAWb[iev]) {
941 auto minDb = checkResults(outStreamB, treeName2AWbulk, fitters_host[iev], vtxGen[iev], genParents[iev], k0dec);
942 meanDAWb += minDb;
943 nfoundAWb++;
944 }
945 }
946
947 ft.setUseAbsDCA(false);
948 ft.setWeightedFinalPCA(false);
949 std::fill(fitters_host.begin(), fitters_host.end(), ft);
950
951 swWb.Start(false);
952 std::vector<int> ncWb(NTest, 0);
953 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncWb, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
954 swWb.Stop();
955
956 for (int iev = 0; iev < NTest; iev++) {
957 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncWb[iev] << " Chi2: " << (ncWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
958 if (ncWb[iev]) {
959 auto minDb = checkResults(outStreamB, treeName2Wbulk, fitters_host[iev], vtxGen[iev], genParents[iev], k0dec);
960 meanDWb += minDb;
961 nfoundWb++;
962 }
963 }
964 // ft.print();
965 meanDAb /= nfoundAb ? nfoundAb : 1;
966 meanDAWb /= nfoundAWb ? nfoundAWb : 1;
967 meanDWb /= nfoundWb ? nfoundWb : 1;
968 LOG(info) << "Bulk-processed " << NTest << " 2-prong vertices: Line : Line";
969 LOG(info) << "2-prongs with abs.dist minimization: eff = " << float(nfoundAb) / NTest
970 << " mean.dist to truth: " << meanDAb << " Total time: " << swAb.CpuTime() * 1000 << " ms";
971 LOG(info) << "2-prongs with abs.dist but wghPCA: eff = " << float(nfoundAWb) / NTest
972 << " mean.dist to truth: " << meanDAWb << " Total time: " << swAWb.CpuTime() * 1000 << " ms";
973 LOG(info) << "2-prongs with wgh.dist minimization: eff = " << float(nfoundWb) / NTest
974 << " mean.dist to truth: " << meanDWb << " Total time: " << swWb.CpuTime() * 1000 << " ms";
975 BOOST_CHECK(nfoundAb > 0.99 * NTest);
976 BOOST_CHECK(nfoundAWb > 0.99 * NTest);
977 BOOST_CHECK(nfoundWb > 0.99 * NTest);
978 BOOST_CHECK(meanDAb < 0.1);
979 BOOST_CHECK(meanDAWb < 0.1);
980 BOOST_CHECK(meanDWb < 0.1);
981 }
982
983 // Bulk-process 3 prongs vertices
984 {
985 LOG(info) << "\n\nBulk-processing 3-prongs";
986 std::vector<int> forceQ{1, 1, 1};
987
988 o2::vertexing::DCAFitterN<3> ft; // 3 prong fitter
989 ft.setBz(bz);
990 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
991 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
992 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
993 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
994 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
995
996 std::vector<o2::vertexing::DCAFitterN<3>> fitters_host(NTest);
997 std::vector<TLorentzVector> genParents(NTest);
998
999 std::string treeName3Abulk = "pr3a", treeName3AWbulk = "pr3aw", treeName3Wbulk = "pr3w";
1000 TStopwatch swAb, swAWb, swWb;
1001 int nfoundAb = 0, nfoundAWb = 0, nfoundWb = 0;
1002 double meanDAb = 0, meanDAWb = 0, meanDWb = 0;
1003 swAb.Stop();
1004 swAWb.Stop();
1005 swWb.Stop();
1006 for (int iev = 0; iev < NTest; iev++) {
1007 std::vector<o2::track::TrackParCov> vc(3);
1008 genParents[iev] = generate(vtxGen[iev], vc, bz, genPHS, dch, dchdec, forceQ);
1009
1010 vctracks[0][iev] = vc[0];
1011 vctracks[1][iev] = vc[1];
1012 vctracks[2][iev] = vc[2];
1013 }
1014
1015 ft.setUseAbsDCA(true);
1016 std::fill(fitters_host.begin(), fitters_host.end(), ft);
1017 swAb.Start(false);
1018 std::vector<int> ncAb(NTest, 0);
1019 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAb, vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
1020 swAb.Stop();
1021 for (int iev = 0; iev < NTest; iev++) {
1022 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAb[iev] << " Chi2: " << (ncAb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
1023 if (ncAb[iev]) {
1024 auto minDb = checkResults(outStreamB, treeName3Abulk, fitters_host[iev], vtxGen[iev], genParents[iev], dchdec);
1025 meanDAb += minDb;
1026 nfoundAb++;
1027 }
1028 }
1029
1030 ft.setUseAbsDCA(true);
1031 ft.setWeightedFinalPCA(true);
1032 std::fill(fitters_host.begin(), fitters_host.end(), ft);
1033
1034 swAWb.Start(false);
1035 std::vector<int> ncAWb(NTest, 0);
1036 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncAWb, vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
1037 swAWb.Stop();
1038 for (int iev = 0; iev < NTest; iev++) {
1039 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAWb[iev] << " Chi2: " << (ncAWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
1040 if (ncAWb[iev]) {
1041 auto minDb = checkResults(outStreamB, treeName3AWbulk, fitters_host[iev], vtxGen[iev], genParents[iev], dchdec);
1042 meanDAWb += minDb;
1043 nfoundAWb++;
1044 }
1045 }
1046
1047 ft.setUseAbsDCA(false);
1048 ft.setWeightedFinalPCA(false);
1049 std::fill(fitters_host.begin(), fitters_host.end(), ft);
1050
1051 swWb.Start(false);
1052 std::vector<int> ncWb(NTest, 0);
1053 device::processBulk(nBlocks, nThreads, nBatches, fitters_host, ncWb, vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
1054 swWb.Stop();
1055 for (int iev = 0; iev < NTest; iev++) {
1056 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncWb[iev] << " Chi2: " << (ncWb[iev] ? fitters_host[iev].getChi2AtPCACandidate(0) : -1);
1057 if (ncWb[iev]) {
1058 auto minDb = checkResults(outStreamB, treeName3Wbulk, fitters_host[iev], vtxGen[iev], genParents[iev], dchdec);
1059 meanDWb += minDb;
1060 nfoundWb++;
1061 }
1062 }
1063
1064 // ft.print();
1065 meanDAb /= nfoundAb ? nfoundAb : 1;
1066 meanDAWb /= nfoundAWb ? nfoundAWb : 1;
1067 meanDWb /= nfoundWb ? nfoundWb : 1;
1068 LOG(info) << "Bulk-processed " << NTest << " 3-prong vertices";
1069 LOG(info) << "3-prongs with abs.dist minimization: eff = " << float(nfoundAb) / NTest
1070 << " mean.dist to truth: " << meanDAb << " Total time: " << swAb.CpuTime() * 1000 << " ms";
1071 LOG(info) << "3-prongs with abs.dist but wghPCA: eff = " << float(nfoundAWb) / NTest
1072 << " mean.dist to truth: " << meanDAWb << " Total time: " << swAWb.CpuTime() * 1000 << " ms";
1073 LOG(info) << "3-prongs with wgh.dist minimization: eff = " << float(nfoundWb) / NTest
1074 << " mean.dist to truth: " << meanDWb << " Total time: " << swWb.CpuTime() * 1000 << " ms";
1075 BOOST_CHECK(nfoundAb > 0.99 * NTest);
1076 BOOST_CHECK(nfoundAWb > 0.99 * NTest);
1077 BOOST_CHECK(nfoundWb > 0.99 * NTest);
1078 BOOST_CHECK(meanDAb < 0.1);
1079 BOOST_CHECK(meanDAWb < 0.1);
1080 BOOST_CHECK(meanDWb < 0.1);
1081 }
1082 outStreamB.Close();
1083}
1084
1085} // namespace vertexing
1086} // namespace o2
DCAFitterN.h
Defintions for N-prongs secondary vertex fit.
i
int32_t i
Definition GPUCommonAlgorithm.h:443
c
uint32_t c
Definition RawData.h:2
debug
std::ostringstream debug
Definition VariantJSONHelpers.h:307
o2::vertexing::DCAFitterN::getChi2AtPCACandidate
float getChi2AtPCACandidate(int cand=0) const
Definition DCAFitterN.h:148
x
GLint GLenum GLint x
Definition glcorearb.h:403
params
GLenum const GLfloat * params
Definition glcorearb.h:272
dst
GLenum GLenum dst
Definition glcorearb.h:1767
o2::its::math_utils::float
const float
Definition MathUtils.h:38
o2::math_utils::rotateZInv
std::tuple< float, float > rotateZInv(float xG, float yG, float snAlp, float csAlp)
Definition Utils.h:142
o2::vertexing::device::process
int process(const int nBlocks, const int nThreads, Fitter &, Tr &... args)
o2::vertexing::device::processBulk
void processBulk(const int nBlocks, const int nThreads, const int nBatches, std::vector< Fitter > &fitters, std::vector< int > &results, std::vector< Tr > &... args)
o2::vertexing::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(DCAFitterNProngsBulk)
Definition testDCAFitterNGPU.cxx:557
o2::vertexing::generate
TLorentzVector generate(Vec3D &vtx, std::vector< o2::track::TrackParCov > &vctr, float bz, TGenPhaseSpace &genPHS, double parMass, const std::vector< double > &dtMass, std::vector< int > forceQ)
Definition testDCAFitterNGPU.cxx:72
o2::vertexing::Vec3D
ROOT::Math::SVector< double, 3 > Vec3D
Definition testDCAFitterNGPU.cxx:31
o2::vertexing::checkResults
float checkResults(o2::utils::TreeStreamRedirector &outs, std::string &treeName, FITTER &fitter, Vec3D &vgen, TLorentzVector &genPar, const std::vector< double > &dtMass)
Definition testDCAFitterNGPU.cxx:34
o2
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Definition BitstreamReader.h:24
LOG
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
BOOST_CHECK
BOOST_CHECK(tree)