Process the event and put new data products into it.
26 {
27 std::vector<ldmx::TrackDeDxMassEstimate> mass_estimates;
28
29 if (!event.
exists(track_collection_, input_pass_name_)) {
30 ldmx_log(error) << "Track collection " << track_collection_ << "_"
31 << input_pass_name_ << " not in event, exiting...";
32 event.add("TrackDeDxMassEstimate", mass_estimates);
33 return;
34 }
35 const std::vector<ldmx::Track> tracks{
36 event.getCollection<
ldmx::Track>(track_collection_, input_pass_name_)};
37
38 int track_type;
39 std::string track_coll_str = track_collection_;
40 std::transform(track_coll_str.begin(), track_coll_str.end(),
41 track_coll_str.begin(), ::tolower);
42
43 bool is_truth = track_coll_str.find("truth") != std::string::npos;
44
45 if (is_truth) {
46 if (track_coll_str.find("tagger") != std::string::npos) {
47 track_type = 1;
48 simhit_collection_ = "TaggerSimHits";
49 } else if (track_coll_str.find("recoil") != std::string::npos) {
50 track_type = 2;
51 simhit_collection_ = "RecoilSimHits";
52 } else {
53 track_type = 0;
54 simhit_collection_ = "";
55 }
56 } else {
57
58 track_type = 4;
59 simhit_collection_ = "";
60 }
61
62
63 std::vector<ldmx::SimTrackerHit> simhits;
64 if (is_truth) {
65 if (!event.
exists(simhit_collection_, input_pass_name_)) {
66 ldmx_log(error) << " SimHit collection (" << simhit_collection_ << "_"
67 << input_pass_name_ << ") does not exists, exiting...";
68 event.add("TrackDeDxMassEstimate", mass_estimates);
69 return;
70 }
72 input_pass_name_);
73 }
74
75
76 for (uint i = 0; i < tracks.size(); i++) {
77 auto track = tracks.at(i);
78
79 auto the_qop = track.getQoP();
80 if (the_qop == 0) {
81 ldmx_log(debug) << "Track " << i << "has zero q/p ";
82 continue;
83 }
84
85 int pdg_id = track.getPdgID();
86 float momentum = 1. / std::abs(the_qop) * 1000;
87 ldmx_log(debug) << "Track " << i << " has momentum " << momentum;
88
90 float sum_dedx_inv2 = 0.;
91 float dedx;
92 int n_hits = 0;
93
94 if (is_truth) {
95
96 for (auto hit : simhits) {
97 if (hit.getTrackID() != track.getTrackID()) continue;
98 if (hit.getEdep() >= 0 && hit.getPathLength() > 0) {
99 dedx = hit.getEdep() / hit.getPathLength() * 10;
100 sum_dedx_inv2 += 1. / (dedx * dedx);
101 n_hits++;
102 }
103 }
104 } else {
105
106 for (auto dedx_meas : track.getDedxMeasurements()) {
107 if (dedx_meas > 0) {
108 dedx = dedx_meas * 10;
109 sum_dedx_inv2 += 1. / (dedx * dedx);
110 n_hits++;
111 }
112 }
113 }
114
115 if (sum_dedx_inv2 == 0) {
116 ldmx_log(debug) << "Track " << i << " has no dEdx measurements";
117 continue;
118 }
119
120
121 float the_ih = 1. / sqrt(1. / n_hits * sum_dedx_inv2);
122
123 float mass = 0.;
124 if (the_ih > fit_res_c_) {
125 mass = momentum * sqrt((the_ih - fit_res_c_) / fit_res_k_);
126 } else {
127 ldmx_log(info) << "Track " << i << " has Ih " << the_ih
128 << " which is less than fit_res_C " << fit_res_c_;
129 mass = -100.;
130 }
131
134 mass_est.
setIh(the_ih);
139 mass_estimates.push_back(mass_est);
140 }
141
142
143 event.add("TrackDeDxMassEstimate", mass_estimates);
144}
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 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.