Process the event and put new data products into it.
27 {
28 std::vector<ldmx::TrackDeDxMassEstimate> mass_estimates;
29
30 if (!event.
exists(track_collection_, input_pass_name_)) {
31 ldmx_log(error) << "Track collection " << track_collection_ << "_"
32 << input_pass_name_ << " not in event, exiting...";
33 event.add(out_collection_, mass_estimates);
34 return;
35 }
36 const std::vector<ldmx::Track> tracks{
37 event.getCollection<
ldmx::Track>(track_collection_, input_pass_name_)};
38
39 int track_type;
40 std::string track_coll_str = track_collection_;
41 std::transform(track_coll_str.begin(), track_coll_str.end(),
42 track_coll_str.begin(), ::tolower);
43
44 bool is_truth = track_coll_str.find("truth") != std::string::npos;
45
46 if (is_truth) {
47 if (track_coll_str.find("tagger") != std::string::npos) {
48 track_type = 1;
49 simhit_collection_ = "TaggerSimHits";
50 } else if (track_coll_str.find("recoil") != std::string::npos) {
51 track_type = 2;
52 simhit_collection_ = "RecoilSimHits";
53 } else {
54 track_type = 0;
55 simhit_collection_ = "";
56 }
57 } else {
58
59 track_type = 4;
60 simhit_collection_ = "";
61 }
62
63
64 std::vector<ldmx::SimTrackerHit> simhits;
65 if (is_truth) {
66 if (!event.
exists(simhit_collection_, input_pass_name_)) {
67 ldmx_log(error) << " SimHit collection (" << simhit_collection_ << "_"
68 << input_pass_name_ << ") does not exists, exiting...";
69 event.add(out_collection_, mass_estimates);
70 return;
71 }
73 input_pass_name_);
74 }
75
76
77 for (uint i = 0; i < tracks.size(); i++) {
78 auto track = tracks.at(i);
79
80 auto the_qop = track.getQoP();
81 if (the_qop == 0) {
82 ldmx_log(debug) << "Track " << i << "has zero q/p ";
83 continue;
84 }
85
86 int pdg_id = track.getPdgID();
87 float momentum = 1. / std::abs(the_qop) * 1000;
88 ldmx_log(debug) << "Track " << i << " has momentum " << momentum;
89
91 float sum_dedx_inv2 = 0.;
92 float dedx;
93 int n_hits = 0;
94
95 if (is_truth) {
96
97 for (auto hit : simhits) {
98 if (hit.getTrackID() != track.getTrackID()) continue;
99 if (hit.getEdep() >= 0 && hit.getPathLength() > 0) {
100 dedx = hit.getEdep() / hit.getPathLength() * 10;
101 sum_dedx_inv2 += 1. / (dedx * dedx);
102 n_hits++;
103 }
104 }
105 } else {
106
107 for (auto dedx_meas : track.getDedxMeasurements()) {
108 if (dedx_meas > 0) {
109 dedx = dedx_meas * 10;
110 sum_dedx_inv2 += 1. / (dedx * dedx);
111 n_hits++;
112 }
113 }
114 }
115
116 if (sum_dedx_inv2 == 0) {
117 ldmx_log(debug) << "Track " << i << " has no dEdx measurements";
118 continue;
119 }
120
121
122 float the_ih = 1. / sqrt(1. / n_hits * sum_dedx_inv2);
123
124 float mass2 = (the_ih - fit_res_c_) * momentum * momentum / fit_res_k_;
125
126 float mass = 0.;
127 if (the_ih > fit_res_c_) {
128 mass = momentum * sqrt((the_ih - fit_res_c_) / fit_res_k_);
129 } else {
130 ldmx_log(info) << "Track " << i << " has Ih " << the_ih
131 << " which is less than fit_res_C " << fit_res_c_;
132 mass = -100.;
133 }
134
137 mass_est.
setIh(the_ih);
143 mass_estimates.push_back(mass_est);
144 }
145
146
147 event.add(out_collection_, mass_estimates);
148}
bool exists(const std::string &name, const std::string &passName, bool unique=true) const
Check for the existence of an object or collection with the given name and pass name in the event.
Represents a simulated tracker hit in the simulation.
Represents the estimated mass of a particle using tracker dE/dx information.
void setNhits(int nhits)
Set the number of hits used in the dEdx calculation.
void setTrackType(int track_type)
Set the type of the track.
void setMassSquared(float mass_squared)
Set the estimated mass squared of the particle/track.
void setPdgId(int pdg_id)
Set the PDG ID of the track.
void setIh(float the_ih)
Set the Ih of the particle/track.
void setMomentum(float momentum)
Set the momentum of the particle/track.
void setMass(float mass)
Set the estimated mass of the particle/track.
void setTrackIndex(int track_index)
Set the index of the track.
Implementation of a track object.