9 #include <ecl/modules/eclTrackBremFinder/BremFindingMatchCompute.h>
12 #include <framework/utilities/Angle.h>
15 #include <mdst/dataobjects/ECLCluster.h>
18 #include <genfit/MeasuredStateOnPlane.h>
23 bool BremFindingMatchCompute::isMatch()
25 auto fitted_state = m_measuredStateOnPlane;
27 auto clusterPosition = m_eclCluster.getClusterPosition();
29 auto fitted_pos = fitted_state.getPos();
30 auto fitted_mom = fitted_state.getMom();
31 auto fitted_dir = fitted_state.getDir();
33 auto cov = fitted_state.get6DCov();
34 const double err_px = cov[3][3];
35 const double err_py = cov[4][4];
36 const double err_pz = cov[5][5];
37 const double px = fitted_mom.X();
38 const double py = fitted_mom.Y();
39 const double pz = fitted_mom.Z();
41 const double square_perp = px * px + py * py ;
42 const double perp = std::sqrt(square_perp);
43 const double square_sum = square_perp + pz * pz;
45 const double err_phi = std::sqrt(pow((py / square_perp), 2) * err_px + pow((px / square_perp), 2) * err_py +
46 (py / square_perp) * (px / square_perp) * cov[3][4]);
47 const double err_theta = std::sqrt(pow(((px * pz) / (square_sum * perp)), 2) * err_px +
48 pow(((py * pz) / (square_sum * perp)), 2) * err_py +
49 pow((perp / square_sum), 2) * err_pz +
50 ((px * pz) / (square_sum * perp)) * ((py * pz) / (square_sum * perp)) * cov[3][4] +
51 ((px * pz) / (square_sum * perp)) * (perp / square_sum) * cov[3][5] +
52 ((py * pz) / (square_sum * perp)) * (perp / square_sum) * cov[4][5]);
55 const auto hit_theta = fitted_mom.Theta();
56 const auto hit_phi = fitted_mom.Phi();
62 hitPhi =
PhiAngle(hit_phi + TMath::TwoPi(), err_phi);
68 if ((clusterPosition - fitted_pos).Phi() >= 0) {
69 clusterPhi =
PhiAngle((clusterPosition - fitted_pos).Phi(), m_eclCluster.getUncertaintyPhi());
71 clusterPhi =
PhiAngle((clusterPosition - fitted_pos).Phi() + TMath::TwoPi(), m_eclCluster.getUncertaintyPhi());
73 ThetaAngle clusterTheta =
ThetaAngle((clusterPosition - fitted_pos).Theta(), m_eclCluster.getUncertaintyTheta());
75 if (clusterPhi.
containsIn(hitPhi, m_clusterAcceptanceFactor) &&
76 clusterTheta.
containsIn(hitTheta, m_clusterAcceptanceFactor)) {
81 clusterV.SetMagThetaPhi(1.0f, clusterTheta.
getAngle(), clusterPhi.
getAngle());
83 auto distV = hitV - clusterV;
85 m_distanceHitCluster = distV.Mag();
88 double deltaTheta = abs(hitV.Y() - clusterV.Y());
89 m_effAcceptanceFactor = deltaTheta / (err_theta + m_eclCluster.getUncertaintyTheta());
90 double deltaPhi = abs(hitV.Z() - clusterV.Z());
91 double effFactor = deltaPhi / (err_phi + m_eclCluster.getUncertaintyPhi());
92 if (effFactor > m_effAcceptanceFactor) {
93 m_effAcceptanceFactor = effFactor;