Process the event and put new data products into it.
81 {
83
84 clearProcessor();
85
87 "SimParticles", sim_particles_pass_name_)};
88
89
90
91
92
93 std::vector<double> gamma_p(3);
94 std::vector<double> gamma_x0(3);
95 if (recoil_from_tracking_) {
96 const auto &recoil_tracks{
97 event.getCollection<
ldmx::Track>(track_collection_, track_pass_name_)};
98
99 for (auto &track : recoil_tracks) {
100
101 if (track.q() == 1 && track.getNhits() == 5) {
102 gamma_x0 = track.getPosition();
103 gamma_p[0] = -1. * track.getMomentum()[0];
104 gamma_p[1] = -1. * track.getMomentum()[1];
105 gamma_p[2] = beam_energy_mev_ - track.getMomentum()[2];
106 }
107 }
108 } else {
109 if (event.
exists(sp_collection_, sp_pass_name_)) {
111 sp_collection_, sp_pass_name_);
112 for (auto const &it : particle_map) {
113 for (auto const &sphit : target_sp_hits) {
114 if (sphit.getPosition()[2] > 0) {
115 if (it.first == sphit.getTrackID()) {
116 if (it.second.getPdgID() == 11 &&
117 inList(it.second.getParents(), 0)) {
118
119
120
121 std::vector<float> x0_float = sphit.getPosition();
122 std::vector<double> x0_double(x0_float.begin(), x0_float.end());
123 gamma_x0 = x0_double;
124 gamma_p[0] = -1. * sphit.getMomentum()[0];
125 gamma_p[1] = -1. * sphit.getMomentum()[1];
126 gamma_p[2] = beam_energy_mev_ - sphit.getMomentum()[2];
127 }
128 }
129 }
130 }
131 }
132 }
133 }
134
135
136 const auto &hcal_rec_hits =
137 event.getCollection<
ldmx::HcalHit>(rec_coll_name_, rec_pass_name_);
138
139 double z_mean = 0.;
140 std::vector<int> layers_hit;
142 if (hit.getEnergy() > 0.) {
144 if (det_id.getSection() != 0) {
145 continue;
146 }
147 n_readout_hits_ += 1;
148 double hit_x = hit.getXPos();
149 double hit_y = hit.getYPos();
150 double hit_z = hit.getZPos();
151 double hit_r = sqrt(hit_x * hit_x + hit_y * hit_y);
152
153 summed_det_ += hit.getEnergy();
154
155 x_mean_ += hit_x * hit.getEnergy();
156 y_mean_ += hit_y * hit.getEnergy();
157 z_mean += hit_z * hit.getEnergy();
158 r_mean_ += hit_r * hit.getEnergy();
159
160
161 if (!(std::find(layers_hit.begin(), layers_hit.end(),
162 det_id.getLayerID()) != layers_hit.end())) {
163 layers_hit.push_back(det_id.getLayerID());
164 }
165
166 double proj_x =
167 gamma_x0[0] + (hit_z - gamma_x0[2]) * gamma_p[0] / gamma_p[2];
168 double proj_y =
169 gamma_x0[1] + (hit_z - gamma_x0[2]) * gamma_p[1] / gamma_p[2];
170
171 r_mean_from_photon_track_ +=
172 hit.getEnergy() * sqrt((hit_x - proj_x) * (hit_x - proj_x) +
173 (hit_y - proj_y) * (hit_y - proj_y));
174
175
176 double closest_point = 9999.;
178 if (hit2.getEnergy() > 0.) {
180 if (det_i_d2.getLayerID() == det_id.getLayerID()) {
181
182
183
184 if (hit.isOrientationX()) {
185 if (abs(hit2.getYPos() - hit_y) > 0) {
186 if (abs(hit2.getYPos() - hit_y) < closest_point) {
187 closest_point = abs(hit2.getYPos() - hit_y);
188 }
189 }
190 }
191 if (hit.isOrientationY()) {
192 if (abs(hit2.getXPos() - hit_x) > 0) {
193 if (abs(hit2.getXPos() - hit_x) < closest_point) {
194 closest_point = abs(hit2.getXPos() - hit_x);
195 }
196 }
197 }
198 }
199 }
200 }
201 if (closest_point > 50.) {
202 iso_hits_ += 1;
203 iso_energy_ += hit.getEnergy();
204 }
205 }
206 }
207
208 n_layers_hit_ = layers_hit.size();
209
210 if (summed_det_ > 0.) {
211 x_mean_ /= summed_det_;
212 y_mean_ /= summed_det_;
213 z_mean /= summed_det_;
214 r_mean_ /= summed_det_;
215
216 r_mean_from_photon_track_ /= summed_det_;
217 }
218
220 if (hit.getEnergy() > 0.) {
222 if (det_id.getSection() == 0) {
223 x_std_ += hit.getEnergy() * (hit.getXPos() - x_mean_) *
224 (hit.getXPos() - x_mean_);
225 y_std_ += hit.getEnergy() * (hit.getYPos() - y_mean_) *
226 (hit.getYPos() - y_mean_);
227 z_std_ += hit.getEnergy() * (hit.getZPos() - z_mean) *
228 (hit.getZPos() - z_mean);
229 }
230 }
231 }
232
233 if (summed_det_ > 0.) {
234 x_std_ = sqrt(x_std_ / summed_det_);
235 y_std_ = sqrt(y_std_ / summed_det_);
236 z_std_ = sqrt(z_std_ / summed_det_);
237 }
238
239 result.setVariables(n_layers_hit_, x_std_, y_std_, z_std_, x_mean_, y_mean_,
240 r_mean_, iso_hits_, iso_energy_, n_readout_hits_,
241 summed_det_, r_mean_from_photon_track_);
242
243 buildBDTFeatureVector(result);
244
245 ldmx::ort::FloatArrays inputs({bdt_features_});
246 float pred =
247 rt_->run({feature_list_name_}, inputs, {"probabilities"})[0].at(1);
248 ldmx_log(info) << " Visibles BDT was ran, score is " << pred;
249
250 event.add(collection_name_, result);
251}
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.
Stores reconstructed hit information from the HCAL.
Implements detector ids for HCal subdetector.
Class representing a simulated particle.
Represents a simulated tracker hit in the simulation.
Implementation of a track object.