From 6fecdd1bb41cbfa88bffa9dac20d8fe601766d25 Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Thu, 12 Sep 2024 19:01:23 +0200 Subject: [PATCH] checkpoint --- Common/DCAFitter/GPU/cuda/DCAFitterN.cu | 42 ++ .../GPU/cuda/test/testDCAFitterNGPU.cxx | 508 ++++++------- .../DCAFitter/include/DCAFitter/DCAFitterN.h | 6 + Common/DCAFitter/test/testDCAFitterN.cxx | 694 +++++++++--------- 4 files changed, 647 insertions(+), 603 deletions(-) diff --git a/Common/DCAFitter/GPU/cuda/DCAFitterN.cu b/Common/DCAFitter/GPU/cuda/DCAFitterN.cu index 68e6aaf988ac8..31fc34614e2ec 100644 --- a/Common/DCAFitter/GPU/cuda/DCAFitterN.cu +++ b/Common/DCAFitter/GPU/cuda/DCAFitterN.cu @@ -103,4 +103,46 @@ int process(o2::vertexing::DCAFitterN<2>& fitter, return result; } +template +int process(o2::vertexing::DCAFitterN<2>&, + const int nBlocks = 1, + const int nThreads = 1, + Tr&... args) +{ + DCAFitterN* ft_device; + std::array tracks_device; + // o2::track::TrackParCov* t1_device; + // o2::track::TrackParCov* t2_device; + int result, *result_device; + + gpuCheckError(cudaMalloc(reinterpret_cast(&ft_device), sizeof(o2::vertexing::DCAFitterN))); + // gpuCheckError(cudaMalloc(reinterpret_cast(&t1_device), sizeof(o2::track::TrackParCov))); + // gpuCheckError(cudaMalloc(reinterpret_cast(&t2_device), sizeof(o2::track::TrackParCov))); + for (int iT{0}; iT < N; ++iT) { + gpuCheckError(cudaMalloc(reinterpret_cast(&(tracks_device[iT])), sizeof(o2::track::TrackParCov))); + } + gpuCheckError(cudaMalloc(reinterpret_cast(&result_device), sizeof(int))); + + gpuCheckError(cudaMemcpy(ft_device, &fitter, sizeof(o2::vertexing::DCAFitterN<2>), cudaMemcpyHostToDevice)); + gpuCheckError(cudaMemcpy(t1_device, &track1, sizeof(o2::track::TrackParCov), cudaMemcpyHostToDevice)); + gpuCheckError(cudaMemcpy(t2_device, &track2, sizeof(o2::track::TrackParCov), cudaMemcpyHostToDevice)); + + kernel::processKernel<<>>(ft_device, t1_device, t2_device, result_device); + + gpuCheckError(cudaPeekAtLastError()); + gpuCheckError(cudaDeviceSynchronize()); + + gpuCheckError(cudaMemcpy(&result, result_device, sizeof(int), cudaMemcpyDeviceToHost)); + gpuCheckError(cudaMemcpy(&fitter, ft_device, sizeof(o2::vertexing::DCAFitterN<2>), cudaMemcpyDeviceToHost)); + gpuCheckError(cudaMemcpy(&track1, t1_device, sizeof(o2::track::TrackParCov), cudaMemcpyDeviceToHost)); + gpuCheckError(cudaMemcpy(&track2, t2_device, sizeof(o2::track::TrackParCov), cudaMemcpyDeviceToHost)); + gpuCheckError(cudaFree(ft_device)); + gpuCheckError(cudaFree(t1_device)); + gpuCheckError(cudaFree(t2_device)); + + gpuCheckError(cudaFree(result_device)); + + return result; +} + } // namespace o2::vertexing::device \ No newline at end of file diff --git a/Common/DCAFitter/GPU/cuda/test/testDCAFitterNGPU.cxx b/Common/DCAFitter/GPU/cuda/test/testDCAFitterNGPU.cxx index 2a629af9c43e8..57172079ccf09 100644 --- a/Common/DCAFitter/GPU/cuda/test/testDCAFitterNGPU.cxx +++ b/Common/DCAFitter/GPU/cuda/test/testDCAFitterNGPU.cxx @@ -41,7 +41,6 @@ float checkResults(o2::utils::TreeStreamRedirector& outs, std::string& treeName, bool useWghDCA = fitter.getWeightedFinalPCA(); for (int ic = 0; ic < nCand; ic++) { const auto& vtx = fitter.getPCACandidate(ic); - LOGP(info, "x: {} y: {} z: {}", vtx[0], vtx[1], vtx[2]); auto df = vgen; df -= vtx; @@ -147,7 +146,7 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) { gRandom->Delete(); gRandom = new TRandom(42); - constexpr int NTest = 1; + constexpr int NTest = 10000; o2::utils::TreeStreamRedirector outStream("dcafitterNTest.root"); TGenPhaseSpace genPHS; @@ -188,16 +187,16 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) for (int iev = 0; iev < NTest; iev++) { auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); - // ft.setUseAbsDCA(true); - // swA.Start(false); - // int ncA = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swA.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncA) { - // auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); - // meanDA += minD; - // nfoundA++; - // } + ft.setUseAbsDCA(true); + swA.Start(false); + int ncA = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swA.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); + if (ncA) { + auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); + meanDA += minD; + nfoundA++; + } ft.setUseAbsDCA(true); ft.setWeightedFinalPCA(true); @@ -211,270 +210,271 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) nfoundAW++; } - // ft.setUseAbsDCA(false); - // ft.setWeightedFinalPCA(false); - // swW.Start(false); - // int ncW = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swW.Stop(); - // LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncW) { - // auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); - // meanDW += minD; - // nfoundW++; - // } + ft.setUseAbsDCA(false); + ft.setWeightedFinalPCA(false); + swW.Start(false); + int ncW = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swW.Stop(); + LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncW) { + auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); + meanDW += minD; + nfoundW++; + } } // ft.print(); - // device::print(ft, 1, 1); - // meanDA /= nfoundA ? nfoundA : 1; + device::print(ft, 1, 1); + meanDA /= nfoundA ? nfoundA : 1; meanDAW /= nfoundA ? nfoundA : 1; - // meanDW /= nfoundW ? nfoundW : 1; - LOG(debug) << "Processed " << NTest << " 2-prong vertices Helix : Helix"; - // LOG(debug) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - // << " mean.dist to truth: " << meanDA << " GPU time: " << swA.CpuTime(); - LOG(debug) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest + meanDW /= nfoundW ? nfoundW : 1; + LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix"; + LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest + << " mean.dist to truth: " << meanDA << " GPU time: " << swA.CpuTime(); + LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest << " mean.dist to truth: " << meanDAW << " GPU time: " << swAW.CpuTime(); - // LOG(debug) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - // << " mean.dist to truth: " << meanDW << " GPU time: " << swW.CpuTime(); - // BOOST_CHECK(nfoundA > 0.99 * NTest); + LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest + << " mean.dist to truth: " << meanDW << " GPU time: " << swW.CpuTime(); + BOOST_CHECK(nfoundA > 0.99 * NTest); BOOST_CHECK(nfoundAW > 0.99 * NTest); - // BOOST_CHECK(nfoundW > 0.99 * NTest); - // BOOST_CHECK(meanDA < 0.1); - // BOOST_CHECK(meanDAW < 0.1); - // BOOST_CHECK(meanDW < 0.1); + BOOST_CHECK(nfoundW > 0.99 * NTest); + BOOST_CHECK(meanDA < 0.1); + BOOST_CHECK(meanDAW < 0.1); + BOOST_CHECK(meanDW < 0.1); } - // // 2 prongs vertices with collinear tracks (gamma conversion) - // { - // LOG(debug) << "Processing 2-prong Helix - Helix case gamma conversion"; - // std::vector forceQ{1, 1}; + // 2 prongs vertices with collinear tracks (gamma conversion) + { + LOG(debug) << "Processing 2-prong Helix - Helix case gamma conversion"; + std::vector forceQ{1, 1}; - // o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter - // ft.setBz(bz); - // ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - // ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - // ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - // ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway - // ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - // ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor + o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter + ft.setBz(bz); + ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway + ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway + ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway + ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway + ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway + ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - // std::string treeName2A = "gpr2a", treeName2AW = "gpr2aw", treeName2W = "gpr2w"; - // TStopwatch swA, swAW, swW; - // int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - // double meanDA = 0, meanDAW = 0, meanDW = 0; - // swA.Stop(); - // swAW.Stop(); - // swW.Stop(); - // for (int iev = 0; iev < NTest; iev++) { - // auto genParent = generate(vtxGen, vctracks, bz, genPHS, gamma, gammadec, forceQ); + std::string treeName2A = "gpr2a", treeName2AW = "gpr2aw", treeName2W = "gpr2w"; + TStopwatch swA, swAW, swW; + int nfoundA = 0, nfoundAW = 0, nfoundW = 0; + double meanDA = 0, meanDAW = 0, meanDW = 0; + swA.Stop(); + swAW.Stop(); + swW.Stop(); + for (int iev = 0; iev < NTest; iev++) { + auto genParent = generate(vtxGen, vctracks, bz, genPHS, gamma, gammadec, forceQ); - // ft.setUseAbsDCA(true); - // swA.Start(false); - // int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swA.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncA) { - // auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, gammadec); - // meanDA += minD; - // nfoundA++; - // } + ft.setUseAbsDCA(true); + swA.Start(false); + int ncA = device::process(ft, vctracks[0], vctracks[1]); + swA.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); + if (ncA) { + auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, gammadec); + meanDA += minD; + nfoundA++; + } - // ft.setUseAbsDCA(true); - // ft.setWeightedFinalPCA(true); - // swAW.Start(false); - // int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swAW.Stop(); - // LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncAW) { - // auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, gammadec); - // meanDAW += minD; - // nfoundAW++; - // } + ft.setUseAbsDCA(true); + ft.setWeightedFinalPCA(true); + swAW.Start(false); + int ncAW = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swAW.Stop(); + LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncAW) { + auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, gammadec); + meanDAW += minD; + nfoundAW++; + } - // ft.setUseAbsDCA(false); - // ft.setWeightedFinalPCA(false); - // swW.Start(false); - // int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swW.Stop(); - // LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncW) { - // auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, gammadec); - // meanDW += minD; - // nfoundW++; - // } - // } - // ft.print(); - // meanDA /= nfoundA ? nfoundA : 1; - // meanDAW /= nfoundA ? nfoundA : 1; - // meanDW /= nfoundW ? nfoundW : 1; - // LOG(debug) << "Processed " << NTest << " 2-prong vertices Helix : Helix from gamma conversion"; - // LOG(debug) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - // << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - // LOG(debug) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - // << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - // LOG(debug) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - // << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - // BOOST_CHECK(nfoundA > 0.99 * NTest); - // BOOST_CHECK(nfoundAW > 0.99 * NTest); - // BOOST_CHECK(nfoundW > 0.99 * NTest); - // BOOST_CHECK(meanDA < 2.1); - // BOOST_CHECK(meanDAW < 2.1); - // BOOST_CHECK(meanDW < 2.1); - // } + ft.setUseAbsDCA(false); + ft.setWeightedFinalPCA(false); + swW.Start(false); + int ncW = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swW.Stop(); + LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncW) { + auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, gammadec); + meanDW += minD; + nfoundW++; + } + } - // // 2 prongs vertices with one of charges set to 0: Helix : Line - // { - // std::vector forceQ{1, 1}; - // LOG(debug) << "Processing 2-prong Helix - Line case"; - // o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter - // ft.setBz(bz); - // ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - // ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - // ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - // ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - // ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor + device::print(ft, 1, 1); + meanDA /= nfoundA ? nfoundA : 1; + meanDAW /= nfoundA ? nfoundA : 1; + meanDW /= nfoundW ? nfoundW : 1; + LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix from gamma conversion"; + LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest + << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); + LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest + << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); + LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest + << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); + BOOST_CHECK(nfoundA > 0.99 * NTest); + BOOST_CHECK(nfoundAW > 0.99 * NTest); + BOOST_CHECK(nfoundW > 0.99 * NTest); + BOOST_CHECK(meanDA < 2.1); + BOOST_CHECK(meanDAW < 2.1); + BOOST_CHECK(meanDW < 2.1); + } - // std::string treeName2A = "pr2aHL", treeName2AW = "pr2awHL", treeName2W = "pr2wHL"; - // TStopwatch swA, swAW, swW; - // int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - // double meanDA = 0, meanDAW = 0, meanDW = 0; - // swA.Stop(); - // swAW.Stop(); - // swW.Stop(); - // for (int iev = 0; iev < NTest; iev++) { - // forceQ[iev % 2] = 1; - // forceQ[1 - iev % 2] = 0; - // auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); + // 2 prongs vertices with one of charges set to 0: Helix : Line + { + std::vector forceQ{1, 1}; + LOG(debug) << "Processing 2-prong Helix - Line case"; + o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter + ft.setBz(bz); + ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway + ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway + ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway + ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway + ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - // ft.setUseAbsDCA(true); - // swA.Start(false); - // int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swA.Stop(); - // LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncA) { - // auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); - // meanDA += minD; - // nfoundA++; - // } + std::string treeName2A = "pr2aHL", treeName2AW = "pr2awHL", treeName2W = "pr2wHL"; + TStopwatch swA, swAW, swW; + int nfoundA = 0, nfoundAW = 0, nfoundW = 0; + double meanDA = 0, meanDAW = 0, meanDW = 0; + swA.Stop(); + swAW.Stop(); + swW.Stop(); + for (int iev = 0; iev < NTest; iev++) { + forceQ[iev % 2] = 1; + forceQ[1 - iev % 2] = 0; + auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); - // ft.setUseAbsDCA(true); - // ft.setWeightedFinalPCA(true); - // swAW.Start(false); - // int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swAW.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncAW) { - // auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); - // meanDAW += minD; - // nfoundAW++; - // } + ft.setUseAbsDCA(true); + swA.Start(false); + int ncA = device::process(ft, vctracks[0], vctracks[1]); + swA.Stop(); + LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); + if (ncA) { + auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); + meanDA += minD; + nfoundA++; + } - // ft.setUseAbsDCA(false); - // ft.setWeightedFinalPCA(false); - // swW.Start(false); - // int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swW.Stop(); - // LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncW) { - // auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); - // meanDW += minD; - // nfoundW++; - // } - // } - // ft.print(); - // meanDA /= nfoundA ? nfoundA : 1; - // meanDAW /= nfoundAW ? nfoundAW : 1; - // meanDW /= nfoundW ? nfoundW : 1; - // LOG(debug) << "Processed " << NTest << " 2-prong vertices: Helix : Line"; - // LOG(debug) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - // << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - // LOG(debug) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - // << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - // LOG(debug) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - // << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - // BOOST_CHECK(nfoundA > 0.99 * NTest); - // BOOST_CHECK(nfoundAW > 0.99 * NTest); - // BOOST_CHECK(nfoundW > 0.99 * NTest); - // BOOST_CHECK(meanDA < 0.1); - // BOOST_CHECK(meanDAW < 0.1); - // BOOST_CHECK(meanDW < 0.1); - // } + ft.setUseAbsDCA(true); + ft.setWeightedFinalPCA(true); + swAW.Start(false); + int ncAW = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swAW.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncAW) { + auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); + meanDAW += minD; + nfoundAW++; + } - // // 2 prongs vertices with both of charges set to 0: Line : Line - // { - // std::vector forceQ{0, 0}; - // LOG(debug) << "Processing 2-prong Line - Line case"; - // o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter - // ft.setBz(bz); - // ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - // ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - // ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - // ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - // ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor + ft.setUseAbsDCA(false); + ft.setWeightedFinalPCA(false); + swW.Start(false); + int ncW = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swW.Stop(); + LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncW) { + auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); + meanDW += minD; + nfoundW++; + } + } + device::print(ft, 1, 1); + meanDA /= nfoundA ? nfoundA : 1; + meanDAW /= nfoundAW ? nfoundAW : 1; + meanDW /= nfoundW ? nfoundW : 1; + LOG(info) << "Processed " << NTest << " 2-prong vertices: Helix : Line"; + LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest + << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); + LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest + << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); + LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest + << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); + BOOST_CHECK(nfoundA > 0.99 * NTest); + BOOST_CHECK(nfoundAW > 0.99 * NTest); + BOOST_CHECK(nfoundW > 0.99 * NTest); + BOOST_CHECK(meanDA < 0.1); + BOOST_CHECK(meanDAW < 0.1); + BOOST_CHECK(meanDW < 0.1); + } - // std::string treeName2A = "pr2aLL", treeName2AW = "pr2awLL", treeName2W = "pr2wLL"; - // TStopwatch swA, swAW, swW; - // int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - // double meanDA = 0, meanDAW = 0, meanDW = 0; - // swA.Stop(); - // swAW.Stop(); - // swW.Stop(); - // for (int iev = 0; iev < NTest; iev++) { - // forceQ[0] = forceQ[1] = 0; - // auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); + // 2 prongs vertices with both of charges set to 0: Line : Line + { + std::vector forceQ{0, 0}; + LOG(debug) << "Processing 2-prong Line - Line case"; + o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter + ft.setBz(bz); + ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway + ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway + ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway + ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway + ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - // ft.setUseAbsDCA(true); - // swA.Start(false); - // int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swA.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncA) { - // auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); - // meanDA += minD; - // nfoundA++; - // } + std::string treeName2A = "pr2aLL", treeName2AW = "pr2awLL", treeName2W = "pr2wLL"; + TStopwatch swA, swAW, swW; + int nfoundA = 0, nfoundAW = 0, nfoundW = 0; + double meanDA = 0, meanDAW = 0, meanDW = 0; + swA.Stop(); + swAW.Stop(); + swW.Stop(); + for (int iev = 0; iev < NTest; iev++) { + forceQ[0] = forceQ[1] = 0; + auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); - // ft.setUseAbsDCA(true); - // ft.setWeightedFinalPCA(true); - // swAW.Start(false); - // int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swAW.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncAW) { - // auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); - // meanDAW += minD; - // nfoundAW++; - // } + ft.setUseAbsDCA(true); + swA.Start(false); + int ncA = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swA.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); + if (ncA) { + auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); + meanDA += minD; + nfoundA++; + } - // ft.setUseAbsDCA(false); - // ft.setWeightedFinalPCA(false); - // swW.Start(false); - // int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swW.Stop(); - // LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncW) { - // auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); - // meanDW += minD; - // nfoundW++; - // } - // } - // ft.print(); - // meanDA /= nfoundA ? nfoundA : 1; - // meanDAW /= nfoundAW ? nfoundAW : 1; - // meanDW /= nfoundW ? nfoundW : 1; - // LOG(debug) << "Processed " << NTest << " 2-prong vertices: Line : Line"; - // LOG(debug) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - // << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - // LOG(debug) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - // << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - // LOG(debug) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - // << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - // BOOST_CHECK(nfoundA > 0.99 * NTest); - // BOOST_CHECK(nfoundAW > 0.99 * NTest); - // BOOST_CHECK(nfoundW > 0.99 * NTest); - // BOOST_CHECK(meanDA < 0.1); - // BOOST_CHECK(meanDAW < 0.1); - // BOOST_CHECK(meanDW < 0.1); - // } + ft.setUseAbsDCA(true); + ft.setWeightedFinalPCA(true); + swAW.Start(false); + int ncAW = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swAW.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncAW) { + auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); + meanDAW += minD; + nfoundAW++; + } + + ft.setUseAbsDCA(false); + ft.setWeightedFinalPCA(false); + swW.Start(false); + int ncW = device::process(ft, vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swW.Stop(); + LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncW) { + auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); + meanDW += minD; + nfoundW++; + } + } + device::print(ft, 1, 1); + meanDA /= nfoundA ? nfoundA : 1; + meanDAW /= nfoundAW ? nfoundAW : 1; + meanDW /= nfoundW ? nfoundW : 1; + LOG(info) << "Processed " << NTest << " 2-prong vertices: Line : Line"; + LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest + << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); + LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest + << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); + LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest + << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); + BOOST_CHECK(nfoundA > 0.99 * NTest); + BOOST_CHECK(nfoundAW > 0.99 * NTest); + BOOST_CHECK(nfoundW > 0.99 * NTest); + BOOST_CHECK(meanDA < 0.1); + BOOST_CHECK(meanDAW < 0.1); + BOOST_CHECK(meanDW < 0.1); + } // // 3 prongs vertices // { diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h index e84d825715610..732bf058bd013 100644 --- a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h +++ b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h @@ -1154,6 +1154,12 @@ int process(o2::vertexing::DCAFitterN<2>&, o2::track::TrackParCov&, const int nBlocks = 1, const int nThreads = 1); + +template +int process(o2::vertexing::DCAFitterN<2>&, + const int nBlocks = 1, + const int nThreads = 1, + Args&... args); } // namespace device } // namespace vertexing } // namespace o2 diff --git a/Common/DCAFitter/test/testDCAFitterN.cxx b/Common/DCAFitter/test/testDCAFitterN.cxx index cec2344583911..c9e502bd4727f 100644 --- a/Common/DCAFitter/test/testDCAFitterN.cxx +++ b/Common/DCAFitter/test/testDCAFitterN.cxx @@ -41,7 +41,6 @@ float checkResults(o2::utils::TreeStreamRedirector& outs, std::string& treeName, bool useWghDCA = fitter.getWeightedFinalPCA(); for (int ic = 0; ic < nCand; ic++) { const auto& vtx = fitter.getPCACandidate(ic); - LOGP(info, "x: {} y: {} z: {}", vtx[0], vtx[1], vtx[2]); auto df = vgen; df -= vtx; @@ -146,9 +145,7 @@ TLorentzVector generate(Vec3D& vtx, std::vector& vctr, f BOOST_AUTO_TEST_CASE(DCAFitterNProngs) { - gRandom->Delete(); - gRandom = new TRandom(42); - constexpr int NTest = 1; + constexpr int NTest = 10; o2::utils::TreeStreamRedirector outStream("dcafitterNTest.root"); TGenPhaseSpace genPHS; @@ -189,17 +186,16 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) for (int iev = 0; iev < NTest; iev++) { auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); - - // ft.setUseAbsDCA(true); - // swA.Start(false); - // int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swA.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncA) { - // auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); - // meanDA += minD; - // nfoundA++; - // } + ft.setUseAbsDCA(true); + swA.Start(false); + int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swA.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); + if (ncA) { + auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); + meanDA += minD; + nfoundA++; + } ft.setUseAbsDCA(true); ft.setWeightedFinalPCA(true); @@ -213,345 +209,345 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs) nfoundAW++; } - // ft.setUseAbsDCA(false); - // ft.setWeightedFinalPCA(false); - // swW.Start(false); - // int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swW.Stop(); - // LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncW) { - // auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); - // meanDW += minD; - // nfoundW++; - // } + ft.setUseAbsDCA(false); + ft.setWeightedFinalPCA(false); + swW.Start(false); + int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swW.Stop(); + LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncW) { + auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); + meanDW += minD; + nfoundW++; + } } - // ft.print(); - // meanDA /= nfoundA ? nfoundA : 1; - // meanDAW /= nfoundA ? nfoundA : 1; - // meanDW /= nfoundW ? nfoundW : 1; + ft.print(); + meanDA /= nfoundA ? nfoundA : 1; + meanDAW /= nfoundA ? nfoundA : 1; + meanDW /= nfoundW ? nfoundW : 1; LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix"; - // LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - // << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - // LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - // << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - // LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - // << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - // BOOST_CHECK(nfoundA > 0.99 * NTest); - // BOOST_CHECK(nfoundAW > 0.99 * NTest); - // BOOST_CHECK(nfoundW > 0.99 * NTest); - // BOOST_CHECK(meanDA < 0.1); - // BOOST_CHECK(meanDAW < 0.1); - // BOOST_CHECK(meanDW < 0.1); + LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest + << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); + LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest + << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); + LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest + << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); + BOOST_CHECK(nfoundA > 0.99 * NTest); + BOOST_CHECK(nfoundAW > 0.99 * NTest); + BOOST_CHECK(nfoundW > 0.99 * NTest); + BOOST_CHECK(meanDA < 0.1); + BOOST_CHECK(meanDAW < 0.1); + BOOST_CHECK(meanDW < 0.1); } // 2 prongs vertices with collinear tracks (gamma conversion) - // { - // LOG(info) << "Processing 2-prong Helix - Helix case gamma conversion"; - // std::vector forceQ{1, 1}; - - // o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter - // ft.setBz(bz); - // ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - // ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - // ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - // ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway - // ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - // ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - - // std::string treeName2A = "gpr2a", treeName2AW = "gpr2aw", treeName2W = "gpr2w"; - // TStopwatch swA, swAW, swW; - // int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - // double meanDA = 0, meanDAW = 0, meanDW = 0; - // swA.Stop(); - // swAW.Stop(); - // swW.Stop(); - // for (int iev = 0; iev < NTest; iev++) { - // auto genParent = generate(vtxGen, vctracks, bz, genPHS, gamma, gammadec, forceQ); - - // ft.setUseAbsDCA(true); - // swA.Start(false); - // int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swA.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncA) { - // auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, gammadec); - // meanDA += minD; - // nfoundA++; - // } - - // ft.setUseAbsDCA(true); - // ft.setWeightedFinalPCA(true); - // swAW.Start(false); - // int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swAW.Stop(); - // LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncAW) { - // auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, gammadec); - // meanDAW += minD; - // nfoundAW++; - // } - - // ft.setUseAbsDCA(false); - // ft.setWeightedFinalPCA(false); - // swW.Start(false); - // int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swW.Stop(); - // LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncW) { - // auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, gammadec); - // meanDW += minD; - // nfoundW++; - // } - // } - // ft.print(); - // meanDA /= nfoundA ? nfoundA : 1; - // meanDAW /= nfoundA ? nfoundA : 1; - // meanDW /= nfoundW ? nfoundW : 1; - // LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix from gamma conversion"; - // LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - // << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - // LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - // << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - // LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - // << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - // BOOST_CHECK(nfoundA > 0.99 * NTest); - // BOOST_CHECK(nfoundAW > 0.99 * NTest); - // BOOST_CHECK(nfoundW > 0.99 * NTest); - // BOOST_CHECK(meanDA < 2.1); - // BOOST_CHECK(meanDAW < 2.1); - // BOOST_CHECK(meanDW < 2.1); - // } - - // // 2 prongs vertices with one of charges set to 0: Helix : Line - // { - // std::vector forceQ{1, 1}; - // LOG(info) << "Processing 2-prong Helix - Line case"; - // o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter - // ft.setBz(bz); - // ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - // ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - // ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - // ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - // ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - - // std::string treeName2A = "pr2aHL", treeName2AW = "pr2awHL", treeName2W = "pr2wHL"; - // TStopwatch swA, swAW, swW; - // int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - // double meanDA = 0, meanDAW = 0, meanDW = 0; - // swA.Stop(); - // swAW.Stop(); - // swW.Stop(); - // for (int iev = 0; iev < NTest; iev++) { - // forceQ[iev % 2] = 1; - // forceQ[1 - iev % 2] = 0; - // auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); - - // ft.setUseAbsDCA(true); - // swA.Start(false); - // int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swA.Stop(); - // LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncA) { - // auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); - // meanDA += minD; - // nfoundA++; - // } - - // ft.setUseAbsDCA(true); - // ft.setWeightedFinalPCA(true); - // swAW.Start(false); - // int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swAW.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncAW) { - // auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); - // meanDAW += minD; - // nfoundAW++; - // } - - // ft.setUseAbsDCA(false); - // ft.setWeightedFinalPCA(false); - // swW.Start(false); - // int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swW.Stop(); - // LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncW) { - // auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); - // meanDW += minD; - // nfoundW++; - // } - // } - // ft.print(); - // meanDA /= nfoundA ? nfoundA : 1; - // meanDAW /= nfoundAW ? nfoundAW : 1; - // meanDW /= nfoundW ? nfoundW : 1; - // LOG(info) << "Processed " << NTest << " 2-prong vertices: Helix : Line"; - // LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - // << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - // LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - // << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - // LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - // << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - // BOOST_CHECK(nfoundA > 0.99 * NTest); - // BOOST_CHECK(nfoundAW > 0.99 * NTest); - // BOOST_CHECK(nfoundW > 0.99 * NTest); - // BOOST_CHECK(meanDA < 0.1); - // BOOST_CHECK(meanDAW < 0.1); - // BOOST_CHECK(meanDW < 0.1); - // } - - // // 2 prongs vertices with both of charges set to 0: Line : Line - // { - // std::vector forceQ{0, 0}; - // LOG(info) << "Processing 2-prong Line - Line case"; - // o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter - // ft.setBz(bz); - // ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - // ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - // ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - // ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - // ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - - // std::string treeName2A = "pr2aLL", treeName2AW = "pr2awLL", treeName2W = "pr2wLL"; - // TStopwatch swA, swAW, swW; - // int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - // double meanDA = 0, meanDAW = 0, meanDW = 0; - // swA.Stop(); - // swAW.Stop(); - // swW.Stop(); - // for (int iev = 0; iev < NTest; iev++) { - // forceQ[0] = forceQ[1] = 0; - // auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); - - // ft.setUseAbsDCA(true); - // swA.Start(false); - // int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swA.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncA) { - // auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); - // meanDA += minD; - // nfoundA++; - // } - - // ft.setUseAbsDCA(true); - // ft.setWeightedFinalPCA(true); - // swAW.Start(false); - // int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swAW.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncAW) { - // auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); - // meanDAW += minD; - // nfoundAW++; - // } - - // ft.setUseAbsDCA(false); - // ft.setWeightedFinalPCA(false); - // swW.Start(false); - // int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - // swW.Stop(); - // LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncW) { - // auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); - // meanDW += minD; - // nfoundW++; - // } - // } - // ft.print(); - // meanDA /= nfoundA ? nfoundA : 1; - // meanDAW /= nfoundAW ? nfoundAW : 1; - // meanDW /= nfoundW ? nfoundW : 1; - // LOG(info) << "Processed " << NTest << " 2-prong vertices: Line : Line"; - // LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - // << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - // LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - // << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - // LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - // << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - // BOOST_CHECK(nfoundA > 0.99 * NTest); - // BOOST_CHECK(nfoundAW > 0.99 * NTest); - // BOOST_CHECK(nfoundW > 0.99 * NTest); - // BOOST_CHECK(meanDA < 0.1); - // BOOST_CHECK(meanDAW < 0.1); - // BOOST_CHECK(meanDW < 0.1); - // } - - // // 3 prongs vertices - // { - // std::vector forceQ{1, 1, 1}; - - // o2::vertexing::DCAFitterN<3> ft; // 3 prong fitter - // ft.setBz(bz); - // ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - // ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - // ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - // ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - // ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - - // std::string treeName3A = "pr3a", treeName3AW = "pr3aw", treeName3W = "pr3w"; - // TStopwatch swA, swAW, swW; - // int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - // double meanDA = 0, meanDAW = 0, meanDW = 0; - // swA.Stop(); - // swAW.Stop(); - // swW.Stop(); - // for (int iev = 0; iev < NTest; iev++) { - // auto genParent = generate(vtxGen, vctracks, bz, genPHS, dch, dchdec, forceQ); - - // ft.setUseAbsDCA(true); - // swA.Start(false); - // int ncA = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES - // swA.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncA) { - // auto minD = checkResults(outStream, treeName3A, ft, vtxGen, genParent, dchdec); - // meanDA += minD; - // nfoundA++; - // } - - // ft.setUseAbsDCA(true); - // ft.setWeightedFinalPCA(true); - // swAW.Start(false); - // int ncAW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES - // swAW.Stop(); - // LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncAW) { - // auto minD = checkResults(outStream, treeName3AW, ft, vtxGen, genParent, dchdec); - // meanDAW += minD; - // nfoundAW++; - // } - - // ft.setUseAbsDCA(false); - // ft.setWeightedFinalPCA(false); - // swW.Start(false); - // int ncW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES - // swW.Stop(); - // LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - // if (ncW) { - // auto minD = checkResults(outStream, treeName3W, ft, vtxGen, genParent, dchdec); - // meanDW += minD; - // nfoundW++; - // } - // } - // ft.print(); - // meanDA /= nfoundA ? nfoundA : 1; - // meanDAW /= nfoundAW ? nfoundAW : 1; - // meanDW /= nfoundW ? nfoundW : 1; - // LOG(info) << "Processed " << NTest << " 3-prong vertices"; - // LOG(info) << "3-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - // << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - // LOG(info) << "3-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - // << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - // LOG(info) << "3-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - // << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - // BOOST_CHECK(nfoundA > 0.99 * NTest); - // BOOST_CHECK(nfoundAW > 0.99 * NTest); - // BOOST_CHECK(nfoundW > 0.99 * NTest); - // BOOST_CHECK(meanDA < 0.1); - // BOOST_CHECK(meanDAW < 0.1); - // BOOST_CHECK(meanDW < 0.1); - // } + { + LOG(info) << "Processing 2-prong Helix - Helix case gamma conversion"; + std::vector forceQ{1, 1}; + + o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter + ft.setBz(bz); + ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway + ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway + ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway + ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway + ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway + ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor + + std::string treeName2A = "gpr2a", treeName2AW = "gpr2aw", treeName2W = "gpr2w"; + TStopwatch swA, swAW, swW; + int nfoundA = 0, nfoundAW = 0, nfoundW = 0; + double meanDA = 0, meanDAW = 0, meanDW = 0; + swA.Stop(); + swAW.Stop(); + swW.Stop(); + for (int iev = 0; iev < NTest; iev++) { + auto genParent = generate(vtxGen, vctracks, bz, genPHS, gamma, gammadec, forceQ); + + ft.setUseAbsDCA(true); + swA.Start(false); + int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swA.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); + if (ncA) { + auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, gammadec); + meanDA += minD; + nfoundA++; + } + + ft.setUseAbsDCA(true); + ft.setWeightedFinalPCA(true); + swAW.Start(false); + int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swAW.Stop(); + LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncAW) { + auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, gammadec); + meanDAW += minD; + nfoundAW++; + } + + ft.setUseAbsDCA(false); + ft.setWeightedFinalPCA(false); + swW.Start(false); + int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swW.Stop(); + LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncW) { + auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, gammadec); + meanDW += minD; + nfoundW++; + } + } + ft.print(); + meanDA /= nfoundA ? nfoundA : 1; + meanDAW /= nfoundA ? nfoundA : 1; + meanDW /= nfoundW ? nfoundW : 1; + LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix from gamma conversion"; + LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest + << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); + LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest + << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); + LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest + << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); + BOOST_CHECK(nfoundA > 0.99 * NTest); + BOOST_CHECK(nfoundAW > 0.99 * NTest); + BOOST_CHECK(nfoundW > 0.99 * NTest); + BOOST_CHECK(meanDA < 2.1); + BOOST_CHECK(meanDAW < 2.1); + BOOST_CHECK(meanDW < 2.1); + } + + // 2 prongs vertices with one of charges set to 0: Helix : Line + { + std::vector forceQ{1, 1}; + LOG(info) << "Processing 2-prong Helix - Line case"; + o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter + ft.setBz(bz); + ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway + ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway + ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway + ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway + ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor + + std::string treeName2A = "pr2aHL", treeName2AW = "pr2awHL", treeName2W = "pr2wHL"; + TStopwatch swA, swAW, swW; + int nfoundA = 0, nfoundAW = 0, nfoundW = 0; + double meanDA = 0, meanDAW = 0, meanDW = 0; + swA.Stop(); + swAW.Stop(); + swW.Stop(); + for (int iev = 0; iev < NTest; iev++) { + forceQ[iev % 2] = 1; + forceQ[1 - iev % 2] = 0; + auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); + + ft.setUseAbsDCA(true); + swA.Start(false); + int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swA.Stop(); + LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); + if (ncA) { + auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); + meanDA += minD; + nfoundA++; + } + + ft.setUseAbsDCA(true); + ft.setWeightedFinalPCA(true); + swAW.Start(false); + int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swAW.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncAW) { + auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); + meanDAW += minD; + nfoundAW++; + } + + ft.setUseAbsDCA(false); + ft.setWeightedFinalPCA(false); + swW.Start(false); + int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swW.Stop(); + LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncW) { + auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); + meanDW += minD; + nfoundW++; + } + } + ft.print(); + meanDA /= nfoundA ? nfoundA : 1; + meanDAW /= nfoundAW ? nfoundAW : 1; + meanDW /= nfoundW ? nfoundW : 1; + LOG(info) << "Processed " << NTest << " 2-prong vertices: Helix : Line"; + LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest + << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); + LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest + << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); + LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest + << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); + BOOST_CHECK(nfoundA > 0.99 * NTest); + BOOST_CHECK(nfoundAW > 0.99 * NTest); + BOOST_CHECK(nfoundW > 0.99 * NTest); + BOOST_CHECK(meanDA < 0.1); + BOOST_CHECK(meanDAW < 0.1); + BOOST_CHECK(meanDW < 0.1); + } + + // 2 prongs vertices with both of charges set to 0: Line : Line + { + std::vector forceQ{0, 0}; + LOG(info) << "Processing 2-prong Line - Line case"; + o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter + ft.setBz(bz); + ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway + ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway + ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway + ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway + ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor + + std::string treeName2A = "pr2aLL", treeName2AW = "pr2awLL", treeName2W = "pr2wLL"; + TStopwatch swA, swAW, swW; + int nfoundA = 0, nfoundAW = 0, nfoundW = 0; + double meanDA = 0, meanDAW = 0, meanDW = 0; + swA.Stop(); + swAW.Stop(); + swW.Stop(); + for (int iev = 0; iev < NTest; iev++) { + forceQ[0] = forceQ[1] = 0; + auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); + + ft.setUseAbsDCA(true); + swA.Start(false); + int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swA.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); + if (ncA) { + auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); + meanDA += minD; + nfoundA++; + } + + ft.setUseAbsDCA(true); + ft.setWeightedFinalPCA(true); + swAW.Start(false); + int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swAW.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncAW) { + auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); + meanDAW += minD; + nfoundAW++; + } + + ft.setUseAbsDCA(false); + ft.setWeightedFinalPCA(false); + swW.Start(false); + int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES + swW.Stop(); + LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncW) { + auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); + meanDW += minD; + nfoundW++; + } + } + ft.print(); + meanDA /= nfoundA ? nfoundA : 1; + meanDAW /= nfoundAW ? nfoundAW : 1; + meanDW /= nfoundW ? nfoundW : 1; + LOG(info) << "Processed " << NTest << " 2-prong vertices: Line : Line"; + LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest + << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); + LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest + << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); + LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest + << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); + BOOST_CHECK(nfoundA > 0.99 * NTest); + BOOST_CHECK(nfoundAW > 0.99 * NTest); + BOOST_CHECK(nfoundW > 0.99 * NTest); + BOOST_CHECK(meanDA < 0.1); + BOOST_CHECK(meanDAW < 0.1); + BOOST_CHECK(meanDW < 0.1); + } + + // 3 prongs vertices + { + std::vector forceQ{1, 1, 1}; + + o2::vertexing::DCAFitterN<3> ft; // 3 prong fitter + ft.setBz(bz); + ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway + ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway + ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway + ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway + ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor + + std::string treeName3A = "pr3a", treeName3AW = "pr3aw", treeName3W = "pr3w"; + TStopwatch swA, swAW, swW; + int nfoundA = 0, nfoundAW = 0, nfoundW = 0; + double meanDA = 0, meanDAW = 0, meanDW = 0; + swA.Stop(); + swAW.Stop(); + swW.Stop(); + for (int iev = 0; iev < NTest; iev++) { + auto genParent = generate(vtxGen, vctracks, bz, genPHS, dch, dchdec, forceQ); + + ft.setUseAbsDCA(true); + swA.Start(false); + int ncA = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES + swA.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); + if (ncA) { + auto minD = checkResults(outStream, treeName3A, ft, vtxGen, genParent, dchdec); + meanDA += minD; + nfoundA++; + } + + ft.setUseAbsDCA(true); + ft.setWeightedFinalPCA(true); + swAW.Start(false); + int ncAW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES + swAW.Stop(); + LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncAW) { + auto minD = checkResults(outStream, treeName3AW, ft, vtxGen, genParent, dchdec); + meanDAW += minD; + nfoundAW++; + } + + ft.setUseAbsDCA(false); + ft.setWeightedFinalPCA(false); + swW.Start(false); + int ncW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES + swW.Stop(); + LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); + if (ncW) { + auto minD = checkResults(outStream, treeName3W, ft, vtxGen, genParent, dchdec); + meanDW += minD; + nfoundW++; + } + } + ft.print(); + meanDA /= nfoundA ? nfoundA : 1; + meanDAW /= nfoundAW ? nfoundAW : 1; + meanDW /= nfoundW ? nfoundW : 1; + LOG(info) << "Processed " << NTest << " 3-prong vertices"; + LOG(info) << "3-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest + << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); + LOG(info) << "3-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest + << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); + LOG(info) << "3-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest + << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); + BOOST_CHECK(nfoundA > 0.99 * NTest); + BOOST_CHECK(nfoundAW > 0.99 * NTest); + BOOST_CHECK(nfoundW > 0.99 * NTest); + BOOST_CHECK(meanDA < 0.1); + BOOST_CHECK(meanDAW < 0.1); + BOOST_CHECK(meanDW < 0.1); + } outStream.Close(); }