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 printf("Failed to randomize ");
138 trc.print();
139 }
140 }
141 } while (!accept);
142
143 return parent;
144}
145
146BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
147{
148 constexpr int NTest = 10000;
149 o2::utils::TreeStreamRedirector outStream("dcafitterNTest.root");
150
151 TGenPhaseSpace genPHS;
152 constexpr double ele = 0.00051;
153 constexpr double gamma = 2 * ele + 1e-6;
154 constexpr double pion = 0.13957;
155 constexpr double k0 = 0.49761;
156 constexpr double kch = 0.49368;
157 constexpr double dch = 1.86965;
158 std::vector<double> gammadec = {ele, ele};
159 std::vector<double> k0dec = {pion, pion};
160 std::vector<double> dchdec = {pion, kch, pion};
161 std::vector<o2::track::TrackParCov> vctracks;
162 Vec3D vtxGen;
163
164 double bz = 5.0;
165 // 2 prongs vertices
166 {
167 LOG(info) << "\n\nProcessing 2-prong Helix - Helix case";
168 std::vector<int> forceQ{1, 1};
169
170 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
171 ft.setBz(bz);
172 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
173 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
174 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
175 ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway
176 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
177 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
178
179 std::string treeName2A = "pr2a", treeName2AW = "pr2aw", treeName2W = "pr2w";
180 TStopwatch swA, swAW, swW;
181 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
182 double meanDA = 0, meanDAW = 0, meanDW = 0;
183 swA.Stop();
184 swAW.Stop();
185 swW.Stop();
186 for (int iev = 0; iev < NTest; iev++) {
187 auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ);
188
189 ft.setUseAbsDCA(true);
190 swA.Start(false);
191 int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
192 swA.Stop();
193 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
194 if (ncA) {
195 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec);
196 meanDA += minD;
197 nfoundA++;
198 }
199
200 ft.setUseAbsDCA(true);
201 ft.setWeightedFinalPCA(true);
202 swAW.Start(false);
203 int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
204 swAW.Stop();
205 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
206 if (ncAW) {
207 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec);
208 meanDAW += minD;
209 nfoundAW++;
210 }
211
212 ft.setUseAbsDCA(false);
213 ft.setWeightedFinalPCA(false);
214 swW.Start(false);
215 int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
216 swW.Stop();
217 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
218 if (ncW) {
219 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec);
220 meanDW += minD;
221 nfoundW++;
222 }
223 }
224 // ft.print();
225 meanDA /= nfoundA ? nfoundA : 1;
226 meanDAW /= nfoundAW ? nfoundAW : 1;
227 meanDW /= nfoundW ? nfoundW : 1;
228 LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix";
229 LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest
230 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
231 LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
232 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
233 LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
234 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
235 BOOST_CHECK(nfoundA > 0.99 * NTest);
236 BOOST_CHECK(nfoundAW > 0.99 * NTest);
237 BOOST_CHECK(nfoundW > 0.99 * NTest);
238 BOOST_CHECK(meanDA < 0.1);
239 BOOST_CHECK(meanDAW < 0.1);
240 BOOST_CHECK(meanDW < 0.1);
241 ft.print();
242 }
243
244 // 2 prongs vertices with collinear tracks (gamma conversion)
245 {
246 LOG(info) << "\n\nProcessing 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 = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
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 = ft.process(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 = ft.process(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 // ft.print();
304 meanDA /= nfoundA ? nfoundA : 1;
305 meanDAW /= nfoundAW ? nfoundAW : 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 minization: eff= " << float(nfoundA) / NTest
309 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
310 LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
311 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
312 LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
313 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
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 ft.print();
321 }
322
323 // 2 prongs vertices with one of charges set to 0: Helix : Line
324 {
325 std::vector<int> forceQ{1, 1};
326 LOG(info) << "\n\nProcessing 2-prong Helix - Line case";
327 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
328 ft.setBz(bz);
329 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
330 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
331 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
332 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
333 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
334
335 std::string treeName2A = "pr2aHL", treeName2AW = "pr2awHL", treeName2W = "pr2wHL";
336 TStopwatch swA, swAW, swW;
337 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
338 double meanDA = 0, meanDAW = 0, meanDW = 0;
339 swA.Stop();
340 swAW.Stop();
341 swW.Stop();
342 for (int iev = 0; iev < NTest; iev++) {
343 forceQ[iev % 2] = 1;
344 forceQ[1 - iev % 2] = 0;
345 auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ);
346
347 ft.setUseAbsDCA(true);
348 swA.Start(false);
349 int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
350 swA.Stop();
351 LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
352 if (ncA) {
353 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec);
354 meanDA += minD;
355 nfoundA++;
356 }
357
358 ft.setUseAbsDCA(true);
359 ft.setWeightedFinalPCA(true);
360 swAW.Start(false);
361 int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
362 swAW.Stop();
363 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
364 if (ncAW) {
365 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec);
366 meanDAW += minD;
367 nfoundAW++;
368 }
369
370 ft.setUseAbsDCA(false);
371 ft.setWeightedFinalPCA(false);
372 swW.Start(false);
373 int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
374 swW.Stop();
375 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
376 if (ncW) {
377 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec);
378 meanDW += minD;
379 nfoundW++;
380 }
381 }
382 // ft.print();
383 meanDA /= nfoundA ? nfoundA : 1;
384 meanDAW /= nfoundAW ? nfoundAW : 1;
385 meanDW /= nfoundW ? nfoundW : 1;
386 LOG(info) << "Processed " << NTest << " 2-prong vertices: Helix : Line";
387 LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest
388 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
389 LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
390 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
391 LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
392 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
393 BOOST_CHECK(nfoundA > 0.99 * NTest);
394 BOOST_CHECK(nfoundAW > 0.99 * NTest);
395 BOOST_CHECK(nfoundW > 0.99 * NTest);
396 BOOST_CHECK(meanDA < 0.1);
397 BOOST_CHECK(meanDAW < 0.1);
398 BOOST_CHECK(meanDW < 0.1);
399 ft.print();
400 }
401
402 // 2 prongs vertices with both of charges set to 0: Line : Line
403 {
404 std::vector<int> forceQ{0, 0};
405 LOG(info) << "\n\nProcessing 2-prong Line - Line case";
406 o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter
407 ft.setBz(bz);
408 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
409 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
410 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
411 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
412 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
413
414 std::string treeName2A = "pr2aLL", treeName2AW = "pr2awLL", treeName2W = "pr2wLL";
415 TStopwatch swA, swAW, swW;
416 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
417 double meanDA = 0, meanDAW = 0, meanDW = 0;
418 swA.Stop();
419 swAW.Stop();
420 swW.Stop();
421 for (int iev = 0; iev < NTest; iev++) {
422 forceQ[0] = forceQ[1] = 0;
423 auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ);
424
425 ft.setUseAbsDCA(true);
426 swA.Start(false);
427 int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
428 swA.Stop();
429 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
430 if (ncA) {
431 auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec);
432 meanDA += minD;
433 nfoundA++;
434 }
435
436 ft.setUseAbsDCA(true);
437 ft.setWeightedFinalPCA(true);
438 swAW.Start(false);
439 int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
440 swAW.Stop();
441 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
442 if (ncAW) {
443 auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec);
444 meanDAW += minD;
445 nfoundAW++;
446 }
447
448 ft.setUseAbsDCA(false);
449 ft.setWeightedFinalPCA(false);
450 swW.Start(false);
451 int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES
452 swW.Stop();
453 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
454 if (ncW) {
455 auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec);
456 meanDW += minD;
457 nfoundW++;
458 }
459 }
460 // ft.print();
461 meanDA /= nfoundA ? nfoundA : 1;
462 meanDAW /= nfoundAW ? nfoundAW : 1;
463 meanDW /= nfoundW ? nfoundW : 1;
464 LOG(info) << "Processed " << NTest << " 2-prong vertices: Line : Line";
465 LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest
466 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
467 LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
468 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
469 LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
470 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
471 BOOST_CHECK(nfoundA > 0.99 * NTest);
472 BOOST_CHECK(nfoundAW > 0.99 * NTest);
473 BOOST_CHECK(nfoundW > 0.99 * NTest);
474 BOOST_CHECK(meanDA < 0.1);
475 BOOST_CHECK(meanDAW < 0.1);
476 BOOST_CHECK(meanDW < 0.1);
477 ft.print();
478 }
479
480 // 3 prongs vertices
481 {
482 LOG(info) << "\n\nProcessing 3-prong vertices";
483 std::vector<int> forceQ{1, 1, 1};
484 o2::vertexing::DCAFitterN<3> ft; // 3 prong fitter
485 ft.setBz(bz);
486 ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway
487 ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway
488 ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway
489 ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway
490 ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor
491
492 std::string treeName3A = "pr3a", treeName3AW = "pr3aw", treeName3W = "pr3w";
493 TStopwatch swA, swAW, swW;
494 int nfoundA = 0, nfoundAW = 0, nfoundW = 0;
495 double meanDA = 0, meanDAW = 0, meanDW = 0;
496 swA.Stop();
497 swAW.Stop();
498 swW.Stop();
499 for (int iev = 0; iev < NTest; iev++) {
500 auto genParent = generate(vtxGen, vctracks, bz, genPHS, dch, dchdec, forceQ);
501
502 ft.setUseAbsDCA(true);
503 swA.Start(false);
504 int ncA = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
505 swA.Stop();
506 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1);
507 if (ncA) {
508 auto minD = checkResults(outStream, treeName3A, ft, vtxGen, genParent, dchdec);
509 meanDA += minD;
510 nfoundA++;
511 }
512
513 ft.setUseAbsDCA(true);
514 ft.setWeightedFinalPCA(true);
515 swAW.Start(false);
516 int ncAW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
517 swAW.Stop();
518 LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1);
519 if (ncAW) {
520 auto minD = checkResults(outStream, treeName3AW, ft, vtxGen, genParent, dchdec);
521 meanDAW += minD;
522 nfoundAW++;
523 }
524
525 ft.setUseAbsDCA(false);
526 ft.setWeightedFinalPCA(false);
527 swW.Start(false);
528 int ncW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES
529 swW.Stop();
530 LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1);
531 if (ncW) {
532 auto minD = checkResults(outStream, treeName3W, ft, vtxGen, genParent, dchdec);
533 meanDW += minD;
534 nfoundW++;
535 }
536 }
537 // ft.print();
538 meanDA /= nfoundA ? nfoundA : 1;
539 meanDAW /= nfoundAW ? nfoundAW : 1;
540 meanDW /= nfoundW ? nfoundW : 1;
541 LOG(info) << "Processed " << NTest << " 3-prong vertices";
542 LOG(info) << "3-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest
543 << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime() * 1000 << " ms";
544 LOG(info) << "3-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest
545 << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime() * 1000 << " ms";
546 LOG(info) << "3-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest
547 << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime() * 1000 << " ms";
548 BOOST_CHECK(nfoundA > 0.99 * NTest);
549 BOOST_CHECK(nfoundAW > 0.99 * NTest);
550 BOOST_CHECK(nfoundW > 0.99 * NTest);
551 BOOST_CHECK(meanDA < 0.1);
552 BOOST_CHECK(meanDAW < 0.1);
553 BOOST_CHECK(meanDW < 0.1);
554 ft.print();
555 }
556 outStream.Close();
557}
558
559} // namespace vertexing
560} // 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::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::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)