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