add photon conversion to the DCAFitter test

Author: shahor02

Common/DCAFitter/test/testDCAFitterN.cxx

@@ -61,7 +61,11 @@ float checkResults(o2::utils::TreeStreamRedirector& outs, std::string& treeName,
          << "genPart=" << genPar << "recPart=" << moth
          << "genX=" << vgen[0] << "genY=" << vgen[1] << "genZ=" << vgen[2]
          << "dx=" << df[0] << "dy=" << df[1] << "dz=" << df[2] << "dst=" << dst
-         << "useAbsDCA=" << absDCA << "useWghDCA=" << useWghDCA << "parent=" << parentTrack << "\n";
+         << "useAbsDCA=" << absDCA << "useWghDCA=" << useWghDCA << "parent=" << parentTrack;
+    for (int i = 0; i < fitter.getNProngs(); i++) {
+      outs << treeName.c_str() << fmt::format("prong{}=", i).c_str() << fitter.getTrack(i, ic);
+    }
+    outs << treeName.c_str() << "\n";
   }
   return distMin;
 }
@@ -145,10 +149,13 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
   o2::utils::TreeStreamRedirector outStream("dcafitterNTest.root");
 
   TGenPhaseSpace genPHS;
+  constexpr double ele = 0.00051;
+  constexpr double gamma = 2 * ele + 1e-6;
   constexpr double pion = 0.13957;
   constexpr double k0 = 0.49761;
   constexpr double kch = 0.49368;
   constexpr double dch = 1.86965;
+  std::vector<double> gammadec = {ele, ele};
   std::vector<double> k0dec = {pion, pion};
   std::vector<double> dchdec = {pion, kch, pion};
   std::vector<o2::track::TrackParCov> vctracks;
@@ -233,6 +240,84 @@ BOOST_AUTO_TEST_CASE(DCAFitterNProngs)
     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<int> 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<int> forceQ{1, 1};