-
Notifications
You must be signed in to change notification settings - Fork 184
/
TrackAnalysisUtils.cc
81 lines (67 loc) · 1.96 KB
/
TrackAnalysisUtils.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#include "TrackAnalysisUtils.h"
#include "SvtxTrack.h"
#include <cmath>
namespace TrackAnalysisUtils
{
TrackAnalysisUtils::DCAPair get_dca(SvtxTrack* track,
Acts::Vector3& vertex)
{
TrackAnalysisUtils::DCAPair pair;
if (!track)
{
std::cout << "You passed a nullptr, can't calculate DCA" << std::endl;
return pair;
}
Acts::Vector3 pos(track->get_x(),
track->get_y(),
track->get_z());
Acts::Vector3 mom(track->get_px(),
track->get_py(),
track->get_pz());
pos -= vertex;
Acts::ActsSquareMatrix<3> posCov;
for (int i = 0; i < 3; ++i)
{
for (int j = 0; j < 3; ++j)
{
posCov(i, j) = track->get_error(i, j);
}
}
Acts::Vector3 r = mom.cross(Acts::Vector3(0., 0., 1.));
float phi = atan2(r(1), r(0));
phi *= -1;
Acts::RotationMatrix3 rot;
Acts::RotationMatrix3 rot_T;
rot(0, 0) = std::cos(phi);
rot(0, 1) = -std::sin(phi);
rot(0, 2) = 0;
rot(1, 0) = std::sin(phi);
rot(1, 1) = std::cos(phi);
rot(1, 2) = 0;
rot(2, 0) = 0;
rot(2, 1) = 0;
rot(2, 2) = 1;
rot_T = rot.transpose();
Acts::Vector3 pos_R = rot * pos;
Acts::ActsSquareMatrix<3> rotCov = rot * posCov * rot_T;
//! First pair is DCA_xy +/- DCA_xy_err
pair.first.first = pos_R(0);
pair.first.second = sqrt(rotCov(0, 0));
//! Second pair is DCA_z +/- DCA_z_err
pair.second.first = pos_R(2);
pair.second.second = sqrt(rotCov(2, 2));
return pair;
}
std::vector<TrkrDefs::cluskey> get_cluster_keys(SvtxTrack* track)
{
std::vector<TrkrDefs::cluskey> out;
for (const auto& seed : {track->get_silicon_seed(), track->get_tpc_seed()})
{
if (seed)
{
std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
}
}
return out;
}
} // namespace TrackAnalysisUtils