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