Project
Loading...
Searching...
No Matches
testDCAFitterN.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 const std::array<float, 3> genPos{static_cast<float>(vgen[0]), static_cast<float>(vgen[1]), static_cast<float>(vgen[2])};
60 const std::array<float, 3> genMom{static_cast<float>(genPar.Px()), static_cast<float>(genPar.Py()), static_cast<float>(genPar.Pz())};
61 o2::track::TrackPar genParentTrack(genPos, genMom, parentTrack.getCharge(), false);
62 genParentTrack.rotateParam(parentTrack.getAlpha());
63 std::array<float, o2::track::kLabCovMatSize> parentCovGlo{};
64 const bool hasParentCovGlo = parentTrack.getCovXYZPxPyPzGlo(parentCovGlo);
65 const double pullX = hasParentCovGlo && parentCovGlo[0] > 0.f ? df[0] / TMath::Sqrt(parentCovGlo[0]) : 0.;
66 const double pullY = hasParentCovGlo && parentCovGlo[2] > 0.f ? df[1] / TMath::Sqrt(parentCovGlo[2]) : 0.;
67 const double pullZ = hasParentCovGlo && parentCovGlo[5] > 0.f ? df[2] / TMath::Sqrt(parentCovGlo[5]) : 0.;
68 outs << treeName.c_str() << "cand=" << ic << "ncand=" << nCand << "nIter=" << nIter << "chi2=" << chi2
69 << "genPart=" << genPar << "recPart=" << moth
70 << "genX=" << vgen[0] << "genY=" << vgen[1] << "genZ=" << vgen[2]
71 << "dx=" << df[0] << "dy=" << df[1] << "dz=" << df[2] << "dst=" << dst
72 << "pullX=" << pullX << "pullY=" << pullY << "pullZ=" << pullZ
73 << "genParentTrack=" << genParentTrack
74 << "useAbsDCA=" << absDCA << "useWghDCA=" << useWghDCA << "parent=" << parentTrack;
75 for (int i = 0; i < fitter.getNProngs(); i++) {
76 outs << treeName.c_str() << fmt::format("prong{}=", i).c_str() << fitter.getTrack(i, ic);
77 }
78 outs << treeName.c_str() << "\n";
79 }
80 return distMin;
81}
82
83TLorentzVector generate(Vec3D& vtx, std::vector<o2::track::TrackParCov>& vctr, float bz,
84 TGenPhaseSpace& genPHS, double parMass, const std::vector<double>& dtMass, std::vector<int> forceQ)
85{
86 const float errYZ = 1e-2, errSlp = 1e-3, errQPT = 2e-2;
87 std::array<float, 15> covm = {
88 errYZ * errYZ,
89 0., errYZ * errYZ,
90 0, 0., errSlp * errSlp,
91 0., 0., 0., errSlp * errSlp,
92 0., 0., 0., 0., errQPT * errQPT};
93 bool accept = true;
94 TLorentzVector parent, d0, d1, d2;
95 do {
96 accept = true;
97 double y = gRandom->Rndm() - 0.5;
98 double pt = 0.1 + gRandom->Rndm() * 3;
99 double mt = TMath::Sqrt(parMass * parMass + pt * pt);
100 double pz = mt * TMath::SinH(y);
101 double phi = gRandom->Rndm() * TMath::Pi() * 2;
102 double en = mt * TMath::CosH(y);
103 double rdec = 10.; // radius of the decay
104 vtx[0] = rdec * TMath::Cos(phi);
105 vtx[1] = rdec * TMath::Sin(phi);
106 vtx[2] = rdec * pz / pt;
107 parent.SetPxPyPzE(pt * TMath::Cos(phi), pt * TMath::Sin(phi), pz, en);
108 int nd = dtMass.size();
109 genPHS.SetDecay(parent, nd, dtMass.data());
110 genPHS.Generate();
111 vctr.clear();
112 float p[4];
113 for (int i = 0; i < nd; i++) {
114 auto* dt = genPHS.GetDecay(i);
115 if (dt->Pt() < 0.05) {
116 accept = false;
117 break;
118 }
119 dt->GetXYZT(p);
120 float s, c, x;
121 std::array<float, 5> params;
122 o2::math_utils::sincos(dt->Phi(), s, c);
123 o2::math_utils::rotateZInv(vtx[0], vtx[1], x, params[0], s, c);
124
125 params[1] = vtx[2];
126 params[2] = 0.; // since alpha = phi
127 params[3] = 1. / TMath::Tan(dt->Theta());
128 params[4] = (i % 2 ? -1. : 1.) / dt->Pt();
129 covm[14] = errQPT * errQPT * params[4] * params[4];
130 //
131 // randomize
132 float r1, r2;
133 gRandom->Rannor(r1, r2);
134 params[0] += r1 * errYZ;
135 params[1] += r2 * errYZ;
136 gRandom->Rannor(r1, r2);
137 params[2] += r1 * errSlp;
138 params[3] += r2 * errSlp;
139 params[4] *= gRandom->Gaus(1., errQPT);
140 if (forceQ[i] == 0) {
141 params[4] = 0.; // impose straight track
142 }
143 auto& trc = vctr.emplace_back(x, dt->Phi(), params, covm);
144 float rad = forceQ[i] == 0 ? 600. : TMath::Abs(1. / trc.getCurvature(bz));
145 if (!trc.propagateTo(trc.getX() + (gRandom->Rndm() - 0.5) * rad * 0.05, bz) ||
146 !trc.rotate(trc.getAlpha() + (gRandom->Rndm() - 0.5) * 0.2)) {
147 LOGP(error, "Failed to randomize ");
148 trc.print();
149 }
150 }
151 } while (!accept);
152
153 return parent;
154}
155
156static constexpr int NFitStatus{14};
157using FitStatusArray = std::array<std::array<int, 3>, NFitStatus>;
158static constexpr const char* FitStatusNames[NFitStatus] = {
159 "None", "Converged", "MaxIter", "NoCrossing", "RejRadius", "RejTrackX", "RejTrackRoughZ", "RejChi2Max",
160 "FailProp", "FailInvConv", "FailInvWeight", "FailInv2ndDeriv", "FailCorrTracks", "FailCloserAlt"};
161inline void printStat(const FitStatusArray& a)
162{
163 LOGP(info, "FitStatus summary : ....A / ..AWD / ...WD (A=abs.dist;AWD=abs.wghPCA.dist;WD=wgh.dist)");
164 for (int i{0}; i < NFitStatus; ++i) {
165 LOGP(info, "{:2d}={:20s}: {:5d} / {:5d} / {:5d}", i, FitStatusNames[i], a[i][0], a[i][1], a[i][2]);
166 }
167 BOOST_CHECK(a[0][0] == 0); // ensure coverage of all possible states
168 BOOST_CHECK(a[0][1] == 0);
169 BOOST_CHECK(a[0][2] == 0);
170}
171
172BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
173{
174 constexpr int NTest = 10000;
175 o2::utils::TreeStreamRedirector outStream("dcafitterNTest.root");
176
177 TGenPhaseSpace genPHS;
178 constexpr double ele = 0.00051;
179 constexpr double gamma = 2 * ele + 1e-6;
180 constexpr double pion = 0.13957;
181 constexpr double k0 = 0.49761;
182 constexpr double kch = 0.49368;
183 constexpr double dch = 1.86965;
184 std::vector<double> gammadec = {ele, ele};
185 std::vector<double> k0dec = {pion, pion};
186 std::vector<double> dchdec = {pion, kch, pion};
187 std::vector<o2::track::TrackParCov> vctracks;
188 FitStatusArray fitstat;
189 Vec3D vtxGen;
190
191 double bz = 5.0;
192 // 2 prongs vertices
193 {
194 LOG(info) << "\n\nProcessing 2-prong Helix - Helix case";
195 std::vector<int> forceQ{1, 1};
196 std::memset(fitstat.data(), 0, sizeof(fitstat));
197
198 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
199 ft.setBz(bz);
200 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
201 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
202 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
203 ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway
204 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
205 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
206
207 std::string treeName2A = "pr2a", treeName2AW = "pr2aw", treeName2W = "pr2w";
208 TStopwatch swA, swAW, swW;
209 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
210 double meanDA = 0, meanDAW = 0, meanDW = 0;
211 swA.Stop();
212 swAW.Stop();
213 swW.Stop();
214 for (int iev = 0; iev < NTest; iev++) {
215 auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ);
216
217 ft.setUseAbsDCA(true);
218 swA.Start(false);
219 int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
220 swA.Stop();
221 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
222 if (ncA) {
223 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec);
224 meanDA += minD;
225 nfoundA++;
226 }
227 ++fitstat[ft.getFitStatus()][0];
228
229 ft.setUseAbsDCA(true);
230 ft.setWeightedFinalPCA(true);
231 swAW.Start(false);
232 int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
233 swAW.Stop();
234 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
235 if (ncAW) {
236 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec);
237 meanDAW += minD;
238 nfoundAW++;
239 }
240 ++fitstat[ft.getFitStatus()][1];
241
242 ft.setUseAbsDCA(false);
243 ft.setWeightedFinalPCA(false);
244 swW.Start(false);
245 int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
246 swW.Stop();
247 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
248 if (ncW) {
249 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec);
250 meanDW += minD;
251 nfoundW++;
252 }
253 ++fitstat[ft.getFitStatus()][2];
254 }
255 // ft.print();
256 meanDA /= nfoundA ? nfoundA : 1;
257 meanDAW /= nfoundAW ? nfoundAW : 1;
258 meanDW /= nfoundW ? nfoundW : 1;
259 LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix";
260 LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest
261 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
262 LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
263 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
264 LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
265 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
266 printStat(fitstat);
267 BOOST_CHECK(nfoundA > 0.99 * NTest);
268 BOOST_CHECK(nfoundAW > 0.99 * NTest);
269 BOOST_CHECK(nfoundW > 0.99 * NTest);
270 BOOST_CHECK(meanDA < 0.1);
271 BOOST_CHECK(meanDAW < 0.1);
272 BOOST_CHECK(meanDW < 0.1);
273 ft.print();
274 }
275
276 // 2 prongs vertices with collinear tracks (gamma conversion)
277 {
278 LOG(info) << "\n\nProcessing 2-prong Helix - Helix case gamma conversion";
279 std::vector<int> forceQ{1, 1};
280 std::memset(fitstat.data(), 0, sizeof(fitstat));
281
282 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
283 ft.setBz(bz);
284 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
285 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
286 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
287 ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway
288 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
289 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
290 ft.setMaxChi2();
291 ft.setCollinear(true);
292
293 std::string treeName2A = "gpr2a", treeName2AW = "gpr2aw", treeName2W = "gpr2w";
294 TStopwatch swA, swAW, swW;
295 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
296 double meanDA = 0, meanDAW = 0, meanDW = 0;
297 swA.Stop();
298 swAW.Stop();
299 swW.Stop();
300 for (int iev = 0; iev < NTest; iev++) {
301 auto genParent = generate(vtxGen, vctracks, bz, genPHS, gamma, gammadec, forceQ);
302
303 ft.setUseAbsDCA(true);
304 swA.Start(false);
305 int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
306 swA.Stop();
307 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
308 if (ncA) {
309 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, gammadec);
310 meanDA += minD;
311 nfoundA++;
312 }
313 ++fitstat[ft.getFitStatus()][0];
314
315 ft.setUseAbsDCA(true);
316 ft.setWeightedFinalPCA(true);
317 swAW.Start(false);
318 int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
319 swAW.Stop();
320 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
321 if (ncAW) {
322 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, gammadec);
323 meanDAW += minD;
324 nfoundAW++;
325 }
326 ++fitstat[ft.getFitStatus()][1];
327
328 ft.setUseAbsDCA(false);
329 ft.setWeightedFinalPCA(false);
330 swW.Start(false);
331 int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
332 swW.Stop();
333 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
334 if (ncW) {
335 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, gammadec);
336 meanDW += minD;
337 nfoundW++;
338 }
339 ++fitstat[ft.getFitStatus()][2];
340 }
341 // ft.print();
342 meanDA /= nfoundA ? nfoundA : 1;
343 meanDAW /= nfoundAW ? nfoundAW : 1;
344 meanDW /= nfoundW ? nfoundW : 1;
345 LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix from gamma conversion";
346 LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest
347 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
348 LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
349 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
350 LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
351 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
352 printStat(fitstat);
353 BOOST_CHECK(nfoundA > 0.99 * NTest);
354 BOOST_CHECK(nfoundAW > 0.99 * NTest);
355 BOOST_CHECK(nfoundW > 0.99 * NTest);
356 BOOST_CHECK(meanDA < 2.1);
357 BOOST_CHECK(meanDAW < 2.1);
358 BOOST_CHECK(meanDW < 2.1);
359 ft.print();
360 }
361
362 // 2 prongs vertices with one of charges set to 0: Helix : Line
363 {
364 std::vector<int> forceQ{1, 1};
365 LOG(info) << "\n\nProcessing 2-prong Helix - Line case";
366 std::memset(fitstat.data(), 0, sizeof(fitstat));
367
368 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
369 ft.setBz(bz);
370 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
371 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
372 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
373 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
374 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
375
376 std::string treeName2A = "pr2aHL", treeName2AW = "pr2awHL", treeName2W = "pr2wHL";
377 TStopwatch swA, swAW, swW;
378 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
379 double meanDA = 0, meanDAW = 0, meanDW = 0;
380 swA.Stop();
381 swAW.Stop();
382 swW.Stop();
383 for (int iev = 0; iev < NTest; iev++) {
384 forceQ[iev % 2] = 1;
385 forceQ[1 - iev % 2] = 0;
386 auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ);
387
388 ft.setUseAbsDCA(true);
389 swA.Start(false);
390 int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
391 swA.Stop();
392 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
393 if (ncA) {
394 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec);
395 meanDA += minD;
396 nfoundA++;
397 }
398 ++fitstat[ft.getFitStatus()][0];
399
400 ft.setUseAbsDCA(true);
401 ft.setWeightedFinalPCA(true);
402 swAW.Start(false);
403 int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
404 swAW.Stop();
405 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
406 if (ncAW) {
407 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec);
408 meanDAW += minD;
409 nfoundAW++;
410 }
411 ++fitstat[ft.getFitStatus()][1];
412
413 ft.setUseAbsDCA(false);
414 ft.setWeightedFinalPCA(false);
415 swW.Start(false);
416 int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
417 swW.Stop();
418 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
419 if (ncW) {
420 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec);
421 meanDW += minD;
422 nfoundW++;
423 }
424 ++fitstat[ft.getFitStatus()][2];
425 }
426 // ft.print();
427 meanDA /= nfoundA ? nfoundA : 1;
428 meanDAW /= nfoundAW ? nfoundAW : 1;
429 meanDW /= nfoundW ? nfoundW : 1;
430 LOG(info) << "Processed " << NTest << " 2-prong vertices: Helix : Line";
431 LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest
432 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
433 LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
434 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
435 LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
436 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
437 printStat(fitstat);
438 BOOST_CHECK(nfoundA > 0.99 * NTest);
439 BOOST_CHECK(nfoundAW > 0.99 * NTest);
440 BOOST_CHECK(nfoundW > 0.99 * NTest);
441 BOOST_CHECK(meanDA < 0.1);
442 BOOST_CHECK(meanDAW < 0.1);
443 BOOST_CHECK(meanDW < 0.1);
444 ft.print();
445 }
446
447 // 2 prongs vertices with both of charges set to 0: Line : Line
448 {
449 std::vector<int> forceQ{0, 0};
450 LOG(info) << "\n\nProcessing 2-prong Line - Line case";
451 std::memset(fitstat.data(), 0, sizeof(fitstat));
452
453 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
454 ft.setBz(bz);
455 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
456 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
457 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
458 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
459 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
460
461 std::string treeName2A = "pr2aLL", treeName2AW = "pr2awLL", treeName2W = "pr2wLL";
462 TStopwatch swA, swAW, swW;
463 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
464 double meanDA = 0, meanDAW = 0, meanDW = 0;
465 swA.Stop();
466 swAW.Stop();
467 swW.Stop();
468 for (int iev = 0; iev < NTest; iev++) {
469 forceQ[0] = forceQ[1] = 0;
470 auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ);
471
472 ft.setUseAbsDCA(true);
473 swA.Start(false);
474 int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
475 swA.Stop();
476 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
477 if (ncA) {
478 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec);
479 meanDA += minD;
480 nfoundA++;
481 }
482 ++fitstat[ft.getFitStatus()][0];
483
484 ft.setUseAbsDCA(true);
485 ft.setWeightedFinalPCA(true);
486 swAW.Start(false);
487 int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
488 swAW.Stop();
489 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
490 if (ncAW) {
491 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec);
492 meanDAW += minD;
493 nfoundAW++;
494 }
495 ++fitstat[ft.getFitStatus()][1];
496
497 ft.setUseAbsDCA(false);
498 ft.setWeightedFinalPCA(false);
499 swW.Start(false);
500 int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
501 swW.Stop();
502 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
503 if (ncW) {
504 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec);
505 meanDW += minD;
506 nfoundW++;
507 }
508 ++fitstat[ft.getFitStatus()][2];
509 }
510 // ft.print();
511 meanDA /= nfoundA ? nfoundA : 1;
512 meanDAW /= nfoundAW ? nfoundAW : 1;
513 meanDW /= nfoundW ? nfoundW : 1;
514 LOG(info) << "Processed " << NTest << " 2-prong vertices: Line : Line";
515 LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest
516 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
517 LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
518 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
519 LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
520 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
521 printStat(fitstat);
522 BOOST_CHECK(nfoundA > 0.99 * NTest);
523 BOOST_CHECK(nfoundAW > 0.99 * NTest);
524 BOOST_CHECK(nfoundW > 0.99 * NTest);
525 BOOST_CHECK(meanDA < 0.1);
526 BOOST_CHECK(meanDAW < 0.1);
527 BOOST_CHECK(meanDW < 0.1);
528 ft.print();
529 }
530
531 // 3 prongs vertices
532 {
533 LOG(info) << "\n\nProcessing 3-prong vertices";
534 std::vector<int> forceQ{1, 1, 1};
535 std::memset(fitstat.data(), 0, sizeof(fitstat));
536
537 o2::vertexing::DCAFitterN<3> ft; // 3 prong fitter
538 ft.setBz(bz);
539 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
540 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
541 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
542 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
543 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
544
545 std::string treeName3A = "pr3a", treeName3AW = "pr3aw", treeName3W = "pr3w";
546 TStopwatch swA, swAW, swW;
547 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
548 double meanDA = 0, meanDAW = 0, meanDW = 0;
549 swA.Stop();
550 swAW.Stop();
551 swW.Stop();
552 for (int iev = 0; iev < NTest; iev++) {
553 auto genParent = generate(vtxGen, vctracks, bz, genPHS, dch, dchdec, forceQ);
554
555 ft.setUseAbsDCA(true);
556 swA.Start(false);
557 int ncA = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
558 swA.Stop();
559 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
560 if (ncA) {
561 auto minD = checkResults(outStream, treeName3A, ft, vtxGen, genParent, dchdec);
562 meanDA += minD;
563 nfoundA++;
564 }
565 ++fitstat[ft.getFitStatus()][0];
566
567 ft.setUseAbsDCA(true);
568 ft.setWeightedFinalPCA(true);
569 swAW.Start(false);
570 int ncAW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
571 swAW.Stop();
572 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
573 if (ncAW) {
574 auto minD = checkResults(outStream, treeName3AW, ft, vtxGen, genParent, dchdec);
575 meanDAW += minD;
576 nfoundAW++;
577 }
578 ++fitstat[ft.getFitStatus()][1];
579
580 ft.setUseAbsDCA(false);
581 ft.setWeightedFinalPCA(false);
582 swW.Start(false);
583 int ncW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
584 swW.Stop();
585 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
586 if (ncW) {
587 auto minD = checkResults(outStream, treeName3W, ft, vtxGen, genParent, dchdec);
588 meanDW += minD;
589 nfoundW++;
590 }
591 ++fitstat[ft.getFitStatus()][2];
592 }
593 // ft.print();
594 meanDA /= nfoundA ? nfoundA : 1;
595 meanDAW /= nfoundAW ? nfoundAW : 1;
596 meanDW /= nfoundW ? nfoundW : 1;
597 LOG(info) << "Processed " << NTest << " 3-prong vertices";
598 LOG(info) << "3-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest
599 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
600 LOG(info) << "3-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
601 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
602 LOG(info) << "3-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
603 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
604 printStat(fitstat);
605 BOOST_CHECK(nfoundA > 0.99 * NTest);
606 BOOST_CHECK(nfoundAW > 0.99 * NTest);
607 BOOST_CHECK(nfoundW > 0.99 * NTest);
608 BOOST_CHECK(meanDA < 0.1);
609 BOOST_CHECK(meanDAW < 0.1);
610 BOOST_CHECK(meanDW < 0.1);
611 ft.print();
612 }
613 outStream.Close();
614}
615
616} // namespace vertexing
617} // namespace o2
DCAFitterN.h
Defintions for N-prongs secondary vertex fit.
debug
std::ostringstream debug
Definition ExpressionJSONHelpers.cxx:83
i
int32_t i
Definition GPUCommonAlgorithm.h:451
c
uint32_t c
Definition RawData.h:2
o2::vertexing::DCAFitterN::getChi2AtPCACandidate
float getChi2AtPCACandidate(int cand=0) const
Definition DCAFitterN.h:181
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
a
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
o2::math_utils::rotateZInv
std::tuple< float, float > rotateZInv(float xG, float yG, float snAlp, float csAlp)
Definition Utils.h:142
o2::vertexing::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(DCAFitterNProngsBulk)
Definition testDCAFitterNGPU.cxx:557
o2::vertexing::FitStatusArray
std::array< std::array< int, 3 >, NFitStatus > FitStatusArray
Definition testDCAFitterN.cxx:157
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::vertexing::printStat
void printStat(const FitStatusArray &a)
Definition testDCAFitterN.cxx:161
o2
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Definition FlattenRestore.h:23
LOG
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
BOOST_CHECK
BOOST_CHECK(tree)