LDMX Software
dqm::VisiblesCutflow Class Reference

Public Member Functions

 VisiblesCutflow (const std::string &name, framework::Process &process)
 
void configure (framework::config::Parameters &parameters) override
 Callback for the EventProcessor to configure itself from the given set of parameters.
 
void analyze (const framework::Event &event) override
 Process the event and make histograms or summaries.
 
- Public Member Functions inherited from framework::Analyzer
 Analyzer (const std::string &name, Process &process)
 Class constructor.
 
virtual void process (Event &event) final
 Processing an event for an Analyzer is calling analyze.
 
virtual void beforeNewRun (ldmx::RunHeader &run_header) final
 Don't allow Analyzers to add parameters to the run header.
 
- Public Member Functions inherited from framework::EventProcessor
 DECLARE_FACTORY (EventProcessor, EventProcessor *, const std::string &, Process &)
 declare that we have a factory for this class
 
 EventProcessor (const std::string &name, Process &process)
 Class constructor.
 
virtual ~EventProcessor ()=default
 Class destructor.
 
virtual void onNewRun (const ldmx::RunHeader &run_header)
 Callback for the EventProcessor to take any necessary action when the run being processed changes.
 
virtual void onFileOpen (EventFile &event_file)
 Callback for the EventProcessor to take any necessary action when a new event input ROOT file is opened.
 
virtual void onFileClose (EventFile &event_file)
 Callback for the EventProcessor to take any necessary action when a event input ROOT file is closed.
 
virtual void onProcessStart ()
 Callback for the EventProcessor to take any necessary action when the processing of events starts, such as creating histograms.
 
virtual void onProcessEnd ()
 Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.
 
template<class T >
const T & getCondition (const std::string &condition_name)
 Access a conditions object for the current event.
 
TDirectory * getHistoDirectory ()
 Access/create a directory in the histogram file for this event processor to create histograms and analysis tuples.
 
void setStorageHint (framework::StorageControl::Hint hint)
 Mark the current event as having the given storage control hint from this module_.
 
void setStorageHint (framework::StorageControl::Hint hint, const std::string &purposeString)
 Mark the current event as having the given storage control hint from this module and the given purpose string.
 
int getLogFrequency () const
 Get the current logging frequency from the process.
 
int getRunNumber () const
 Get the run number from the process.
 
std::string getName () const
 Get the processor name.
 
void createHistograms (const std::vector< framework::config::Parameters > &histos)
 Internal function which is used to create histograms passed from the python configuration @parma histos vector of Parameters that configure histograms to create.
 

Private Member Functions

bool inList (std::vector< int > parents, int track_id)
 

Private Attributes

std::unique_ptr< ldmx::ort::ONNXRuntimert_
 
double bdt_cut_val_ {0.}
 
std::string feature_list_name_
 
bool all_cuts_ {true}
 
double beam_energy_mev_ {0.}
 
std::string hcal_rec_collection_
 
std::string hcal_rec_pass_name_
 
std::string ecal_rec_collection_
 
std::string ecal_rec_pass_name_
 
bool recoil_from_tracking_ {false}
 
std::string track_collection_
 
std::string track_pass_name_
 
std::string sp_collection_
 
std::string sp_pass_name_
 
std::string sim_particles_coll_name_
 
std::string sim_particles_pass_name_
 
std::string ecal_veto_collection_
 
std::string ecal_veto_pass_
 
double ecal_bdt_cut_val_ {0.}
 

Additional Inherited Members

- Protected Member Functions inherited from framework::EventProcessor
void abortEvent ()
 Abort the event immediately.
 
- Protected Attributes inherited from framework::EventProcessor
HistogramPool histograms_
 helper object for making and filling histograms
 
NtupleManagerntuple_ {NtupleManager::getInstance()}
 Manager for any ntuples.
 
logging::logger the_log_
 The logger for this EventProcessor.
 

Detailed Description

Definition at line 17 of file VisiblesCutflow.h.

Constructor & Destructor Documentation

◆ VisiblesCutflow()

dqm::VisiblesCutflow::VisiblesCutflow ( const std::string & name,
framework::Process & process )
inline

Definition at line 19 of file VisiblesCutflow.h.

20 : Analyzer(name, process) {}
virtual void process(Event &event) final
Processing an event for an Analyzer is calling analyze.
Analyzer(const std::string &name, Process &process)
Class constructor.

Member Function Documentation

◆ analyze()

void dqm::VisiblesCutflow::analyze ( const framework::Event & event)
overridevirtual

Process the event and make histograms or summaries.

Parameters
eventThe Event to analyze

Implements framework::Analyzer.

Definition at line 61 of file VisiblesCutflow.cxx.

61 {
62 std::vector<float> bdt_features;
63
64 double decay_z;
65
66 const auto &particle_map{event.getMap<int, ldmx::SimParticle>(
67 sim_particles_coll_name_, sim_particles_pass_name_)};
68
69 for (auto const &it : particle_map) {
70 std::vector<int> parents = it.second.getParents();
71 for (int track_id : parents) {
72 if (track_id == 0 && it.second.getPdgID() == -11) {
73 decay_z = it.second.getVertex()[2];
74 }
75 }
76 }
77
78 histograms_.fill("total_events", decay_z);
79
80 // Get target scoring plane hits for recoil electron
81 // Use this to calculate the projected photon line vector
82 // This currently uses truth-level information, but it should be replaced
83 // by reconstructed tracker information, when available
84 std::vector<double> gamma_p(3);
85 double gamma_energy = 90000.;
86 std::vector<double> gamma_x0(3);
87
88 std::vector<double> recoil_p(3);
89 bool found_recoil_e = false;
90
91 if (recoil_from_tracking_) {
92 const auto &recoil_tracks{
93 event.getCollection<ldmx::Track>(track_collection_, track_pass_name_)};
94 for (auto &track : recoil_tracks) {
95 // need to figure out how to best isolate candidate electron track
96 if (track.q() == 1 && track.getNhits() == 5) {
97 gamma_x0 = track.getPosition();
98 gamma_p[0] = -1. * track.getMomentum()[0];
99 gamma_p[1] = -1. * track.getMomentum()[1];
100 gamma_p[2] = 8000. - track.getMomentum()[2];
101 }
102 }
103 } else {
104 if (event.exists(sp_collection_, sp_pass_name_)) {
105 const auto &target_sp_hits = event.getCollection<ldmx::SimTrackerHit>(
106 sp_collection_, sp_pass_name_);
107 bool found_rec = false;
108 for (auto const &it : particle_map) {
109 for (auto const &sphit : target_sp_hits) {
110 if (sphit.getPosition()[2] > 0) {
111 if (it.first == sphit.getTrackID()) {
112 if (it.second.getPdgID() == 622) {
113 std::vector<float> x0_float = sphit.getPosition();
114 std::vector<double> x0_double(x0_float.begin(), x0_float.end());
115 gamma_x0 = x0_double;
116 gamma_p = sphit.getMomentum();
117 gamma_energy = sphit.getEnergy();
118 found_rec = true;
119 }
120 if (it.second.getPdgID() == 11 &&
121 inList(it.second.getParents(), 0)) {
122 if (!found_rec) {
123 std::vector<float> x0_float = sphit.getPosition();
124 std::vector<double> x0_double(x0_float.begin(),
125 x0_float.end());
126 gamma_x0 = x0_double;
127 gamma_p[0] = -1. * sphit.getMomentum()[0];
128 gamma_p[1] = -1. * sphit.getMomentum()[1];
129 gamma_p[2] = beam_energy_mev_ - sphit.getMomentum()[2];
130 gamma_energy = beam_energy_mev_ - sphit.getEnergy();
131 found_rec = true;
132 }
133 recoil_p = sphit.getMomentum();
134 found_recoil_e = true;
135 }
136 }
137 }
138 }
139 }
140 }
141 }
142 histograms_.fill("beam_energy_frac", gamma_energy / 8000.);
143 histograms_.fill("beam_angle",
144 std::acos(gamma_p[2] / std::sqrt(gamma_p[0] * gamma_p[0] +
145 gamma_p[1] * gamma_p[1] +
146 gamma_p[2] * gamma_p[2])));
147
148 double p_mag = 0.;
149 if (found_recoil_e) {
150 p_mag = std::sqrt(recoil_p[0] * recoil_p[0] + recoil_p[1] * recoil_p[1] +
151 recoil_p[2] * recoil_p[2]);
152 }
153
154 if (p_mag < 50. && all_cuts_) {
155 return;
156 }
157 histograms_.fill("pass_acceptance", decay_z);
158
159 // Get EcalRecHits, check that trigger is passed
160 const auto &ecal_rec_hits = event.getCollection<ldmx::EcalHit>(
161 ecal_rec_collection_, ecal_rec_pass_name_);
162 const auto &hcal_rec_hits = event.getCollection<ldmx::HcalHit>(
163 hcal_rec_collection_, hcal_rec_pass_name_);
164
165 double trigger = 0.;
166 double ecal_energy = 0.;
167 double hcal_energy = 0.;
168 bool hcal_containment = true;
169
170 for (const ldmx::EcalHit &hit : ecal_rec_hits) {
171 if (hit.getEnergy() > 0.) {
172 ecal_energy += hit.getEnergy();
173 if (hit.getZPos() <= 541.722) {
174 trigger += hit.getEnergy();
175 }
176 }
177 }
178 for (const ldmx::HcalHit &hit : hcal_rec_hits) {
179 if (hit.getEnergy() > 0.) {
180 ldmx::HcalID det_id(hit.getID());
181 if (det_id.getSection() != 0) {
182 continue;
183 }
184 if (det_id.getLayerID() == 1 && hit.getPE() > 5) {
185 hcal_containment = false;
186 }
187 hcal_energy += 12. * hit.getEnergy();
188 }
189 }
190
191 if (trigger > 3160. && all_cuts_) {
192 return;
193 }
194 histograms_.fill("pass_trigger", decay_z);
195
196 if (ecal_energy > 3160. && all_cuts_) {
197 return;
198 }
199 histograms_.fill("pass_ecal_energy", decay_z);
200
201 if (p_mag > 2400. && all_cuts_) {
202 return;
203 }
204 histograms_.fill("pass_tracker_veto", decay_z);
205
206 const auto &ecal_veto{event.getObject<ldmx::EcalVetoResult>(
207 ecal_veto_collection_, ecal_veto_pass_)};
208 if (ecal_veto.getDisc() < ecal_bdt_cut_val_ && all_cuts_) {
209 return;
210 }
211 histograms_.fill("pass_ecal_bdt", decay_z);
212
213 if (hcal_energy < 4840 && all_cuts_) {
214 return;
215 }
216 histograms_.fill("pass_hcal_energy", decay_z);
217
218 if (!hcal_containment && all_cuts_) {
219 return;
220 }
221 histograms_.fill("pass_hcal_containment", decay_z);
222
223 // initialize all of the features
224 int n_layers_hit = 0;
225 double x_std = 0.;
226 double y_std = 0.;
227 double z_std = 0.;
228 double x_mean = 0.;
229 double y_mean = 0.;
230 double r_mean = 0.;
231 int iso_hits = 0;
232 double iso_energy = 0.;
233 int n_readout_hits = 0;
234 double summed_det = 0.;
235 double r_mean_from_photon_track = 0.;
236
237 double z_mean = 0.; // need this when calculating z_std
238 std::vector<int> layers_hit;
239 for (const ldmx::HcalHit &hit : hcal_rec_hits) {
240 if (hit.getEnergy() > 0.) {
241 ldmx::HcalID det_id(hit.getID());
242 if (det_id.getSection() != 0) { // skip hits that aren't in main Hcal
243 continue;
244 }
245 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
246 continue;
247 }
248 n_readout_hits += 1;
249 double hit_x = hit.getXPos();
250 double hit_y = hit.getYPos();
251 double hit_z = hit.getZPos();
252 double hit_r = sqrt(hit_x * hit_x + hit_y * hit_y);
253
254 summed_det += hit.getEnergy();
255
256 x_mean += hit_x * hit.getEnergy();
257 y_mean += hit_y * hit.getEnergy();
258 z_mean += hit_z * hit.getEnergy();
259 r_mean += hit_r * hit.getEnergy();
260
261 // check if this is a new layer in the collection
262 if (!(std::find(layers_hit.begin(), layers_hit.end(),
263 det_id.getLayerID()) != layers_hit.end())) {
264 layers_hit.push_back(det_id.getLayerID());
265 }
266
267 double proj_x =
268 gamma_x0[0] + (hit_z - gamma_x0[2]) * gamma_p[0] / gamma_p[2];
269 double proj_y =
270 gamma_x0[1] + (hit_z - gamma_x0[2]) * gamma_p[1] / gamma_p[2];
271
272 r_mean_from_photon_track +=
273 hit.getEnergy() * sqrt((hit_x - proj_x) * (hit_x - proj_x) +
274 (hit_y - proj_y) * (hit_y - proj_y));
275
276 // Calculate isolated hits
277 double closest_point = 9999.;
278 for (const ldmx::HcalHit &hit2 : hcal_rec_hits) {
279 if (hit2.getEnergy() > 0.) {
280 ldmx::HcalID det_i_d2(hit2.getID());
281 if (fabs(hit2.getXPos()) > 1000 || fabs(hit2.getYPos()) > 1000) {
282 continue;
283 }
284 if (det_i_d2.getLayerID() == det_id.getLayerID()) {
285 // Determine if a bar is vertical (along y-axis) or horizontal
286 // (along x-axis) Odd layers have horizontal strips Even layers have
287 // vertical strips
288 if (hit.isOrientationX()) {
289 if (fabs(hit2.getYPos() - hit_y) > 0) {
290 if (fabs(hit2.getYPos() - hit_y) < closest_point) {
291 closest_point = abs(hit2.getYPos() - hit_y);
292 }
293 }
294 }
295 if (hit.isOrientationY()) {
296 if (fabs(hit2.getXPos() - hit_x) > 0) {
297 if (fabs(hit2.getXPos() - hit_x) < closest_point) {
298 closest_point = fabs(hit2.getXPos() - hit_x);
299 }
300 }
301 }
302 }
303 }
304 }
305 if (closest_point > 50.) {
306 iso_hits += 1;
307 iso_energy += hit.getEnergy();
308 }
309 }
310 }
311
312 n_layers_hit = layers_hit.size();
313
314 if (summed_det > 0.) {
315 x_mean /= summed_det;
316 y_mean /= summed_det;
317 z_mean /= summed_det;
318 r_mean /= summed_det;
319
320 r_mean_from_photon_track /= summed_det;
321 }
322
323 for (const ldmx::HcalHit &hit : hcal_rec_hits) {
324 if (hit.getEnergy() > 0.) {
325 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
326 continue;
327 }
328 ldmx::HcalID det_id(hit.getID());
329 if (det_id.getSection() == 0) {
330 x_std += hit.getEnergy() * (hit.getXPos() - x_mean) *
331 (hit.getXPos() - x_mean);
332 y_std += hit.getEnergy() * (hit.getYPos() - y_mean) *
333 (hit.getYPos() - y_mean);
334 z_std += hit.getEnergy() * (hit.getZPos() - z_mean) *
335 (hit.getZPos() - z_mean);
336 }
337 }
338 }
339
340 if (summed_det > 0.) {
341 x_std = sqrt(x_std / summed_det);
342 y_std = sqrt(y_std / summed_det);
343 z_std = sqrt(z_std / summed_det);
344 }
345
346 // Fill histograms
347 histograms_.fill("layers_hit", n_layers_hit);
348 histograms_.fill("x_std", x_std);
349 histograms_.fill("y_std", y_std);
350 histograms_.fill("z_std", z_std);
351 histograms_.fill("x_mean", x_mean);
352 histograms_.fill("y_mean", y_mean);
353 histograms_.fill("r_mean", r_mean);
354 histograms_.fill("iso_hits", iso_hits);
355 histograms_.fill("iso_energy", iso_energy);
356 histograms_.fill("n_hits", n_readout_hits);
357 histograms_.fill("total_energy", hcal_energy);
358 histograms_.fill("photon_track", r_mean_from_photon_track);
359
360 bdt_features.push_back(n_layers_hit);
361 bdt_features.push_back(x_std);
362 bdt_features.push_back(y_std);
363 bdt_features.push_back(z_std);
364 bdt_features.push_back(x_mean);
365 bdt_features.push_back(y_mean);
366 bdt_features.push_back(r_mean);
367 bdt_features.push_back(iso_hits);
368 bdt_features.push_back(iso_energy);
369 bdt_features.push_back(n_readout_hits);
370 bdt_features.push_back(hcal_energy);
371 bdt_features.push_back(r_mean_from_photon_track);
372
373 ldmx::ort::FloatArrays inputs({bdt_features});
374 float pred =
375 rt_->run({feature_list_name_}, inputs, {"probabilities"})[0].at(1);
376
377 histograms_.fill("visibles_disc", pred);
378
379 for (int i = 0; i < 10000; i++) {
380 double disc = 0.999 + (double(i) / 10000.) * (1. - 0.999);
381 histograms_.fill("visibles_disc_high_norm",
382 0.999 + ((double(i) + 0.5) / 10000.) * (1. - 0.999));
383 double disc_roc = 0.99 + (double(i) / 10000.) * (1. - 0.99);
384 if (pred >= disc) {
385 histograms_.fill("visibles_disc_high",
386 0.999 + ((double(i) + 0.5) / 10000.) * (1. - 0.999));
387 }
388 if (pred >= disc_roc) {
389 histograms_.fill("roc", disc_roc);
390 }
391 }
392
393 histograms_.fill("ecal_disc_vs_vis_disc", pred, ecal_veto.getDisc());
394
395 if (pred < bdt_cut_val_) {
396 return;
397 }
398 histograms_.fill("pass_visibles_bdt", decay_z);
399
400 return;
401}
HistogramPool histograms_
helper object for making and filling histograms
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.
Definition Event.cxx:92
void fill(const std::string &name, const T &val)
Fill a 1D histogram.
Stores reconstructed hit information from the ECAL.
Definition EcalHit.h:19
Stores reconstructed hit information from the HCAL.
Definition HcalHit.h:24
Implements detector ids for HCal subdetector.
Definition HcalID.h:19
Class representing a simulated particle.
Definition SimParticle.h:23
Represents a simulated tracker hit in the simulation.
Implementation of a track object.
Definition Track.h:53

References framework::Event::exists(), framework::HistogramPool::fill(), ldmx::HcalID::getLayerID(), and framework::EventProcessor::histograms_.

◆ configure()

void dqm::VisiblesCutflow::configure ( framework::config::Parameters & parameters)
overridevirtual

Callback for the EventProcessor to configure itself from the given set of parameters.

The parameters a processor has access to are the member variables of the python class in the sequence that has class_name equal to the EventProcessor class name.

For an example, look at MyProcessor.

Parameters
parametersParameters for configuration.

Reimplemented from framework::EventProcessor.

Definition at line 24 of file VisiblesCutflow.cxx.

24 {
25 rt_ = std::make_unique<ldmx::ort::ONNXRuntime>(
26 parameters.get<std::string>("bdt_file"));
27 feature_list_name_ = parameters.get<std::string>("feature_list_name");
28 bdt_cut_val_ = parameters.get<double>("disc_cut");
29 ecal_bdt_cut_val_ = parameters.get<double>("ecal_disc_cut");
30 all_cuts_ = parameters.get<bool>("all_cuts");
31
32 beam_energy_mev_ = parameters.get<double>("beam_energy");
33
34 // collection names
35 hcal_rec_collection_ = parameters.get<std::string>("hcal_rec_coll_name");
36 hcal_rec_pass_name_ = parameters.get<std::string>("hcal_rec_pass_name");
37
38 ecal_rec_collection_ = parameters.get<std::string>("ecal_rec_coll_name");
39 ecal_rec_pass_name_ = parameters.get<std::string>("ecal_rec_pass_name");
40
41 recoil_from_tracking_ = parameters.get<bool>("recoil_from_tracking");
42 track_collection_ = parameters.get<std::string>("track_collection");
43 track_pass_name_ = parameters.get<std::string>("track_pass_name");
44
45 sp_collection_ = parameters.get<std::string>("sp_coll_name");
46 sp_pass_name_ = parameters.get<std::string>("sp_pass_name");
47
48 sim_particles_coll_name_ =
49 parameters.get<std::string>("sim_particles_coll_name");
50 sim_particles_pass_name_ =
51 parameters.get<std::string>("sim_particles_pass_name");
52
53 ecal_veto_collection_ = parameters.get<std::string>("ecal_veto_coll_name");
54 ecal_veto_pass_ = parameters.get<std::string>("ecal_veto_pass_name");
55}
const T & get(const std::string &name) const
Retrieve the parameter of the given name.
Definition Parameters.h:78

References framework::config::Parameters::get().

◆ inList()

bool dqm::VisiblesCutflow::inList ( std::vector< int > parents,
int track_id )
private

Definition at line 57 of file VisiblesCutflow.cxx.

57 {
58 return std::find(parents.begin(), parents.end(), track_id) != parents.end();
59}

Member Data Documentation

◆ all_cuts_

bool dqm::VisiblesCutflow::all_cuts_ {true}
private

Definition at line 33 of file VisiblesCutflow.h.

33{true};

◆ bdt_cut_val_

double dqm::VisiblesCutflow::bdt_cut_val_ {0.}
private

Definition at line 30 of file VisiblesCutflow.h.

30{0.};

◆ beam_energy_mev_

double dqm::VisiblesCutflow::beam_energy_mev_ {0.}
private

Definition at line 35 of file VisiblesCutflow.h.

35{0.};

◆ ecal_bdt_cut_val_

double dqm::VisiblesCutflow::ecal_bdt_cut_val_ {0.}
private

Definition at line 50 of file VisiblesCutflow.h.

50{0.};

◆ ecal_rec_collection_

std::string dqm::VisiblesCutflow::ecal_rec_collection_
private

Definition at line 39 of file VisiblesCutflow.h.

◆ ecal_rec_pass_name_

std::string dqm::VisiblesCutflow::ecal_rec_pass_name_
private

Definition at line 40 of file VisiblesCutflow.h.

◆ ecal_veto_collection_

std::string dqm::VisiblesCutflow::ecal_veto_collection_
private

Definition at line 48 of file VisiblesCutflow.h.

◆ ecal_veto_pass_

std::string dqm::VisiblesCutflow::ecal_veto_pass_
private

Definition at line 49 of file VisiblesCutflow.h.

◆ feature_list_name_

std::string dqm::VisiblesCutflow::feature_list_name_
private

Definition at line 31 of file VisiblesCutflow.h.

◆ hcal_rec_collection_

std::string dqm::VisiblesCutflow::hcal_rec_collection_
private

Definition at line 37 of file VisiblesCutflow.h.

◆ hcal_rec_pass_name_

std::string dqm::VisiblesCutflow::hcal_rec_pass_name_
private

Definition at line 38 of file VisiblesCutflow.h.

◆ recoil_from_tracking_

bool dqm::VisiblesCutflow::recoil_from_tracking_ {false}
private

Definition at line 41 of file VisiblesCutflow.h.

41{false};

◆ rt_

std::unique_ptr<ldmx::ort::ONNXRuntime> dqm::VisiblesCutflow::rt_
private

Definition at line 29 of file VisiblesCutflow.h.

◆ sim_particles_coll_name_

std::string dqm::VisiblesCutflow::sim_particles_coll_name_
private

Definition at line 46 of file VisiblesCutflow.h.

◆ sim_particles_pass_name_

std::string dqm::VisiblesCutflow::sim_particles_pass_name_
private

Definition at line 47 of file VisiblesCutflow.h.

◆ sp_collection_

std::string dqm::VisiblesCutflow::sp_collection_
private

Definition at line 44 of file VisiblesCutflow.h.

◆ sp_pass_name_

std::string dqm::VisiblesCutflow::sp_pass_name_
private

Definition at line 45 of file VisiblesCutflow.h.

◆ track_collection_

std::string dqm::VisiblesCutflow::track_collection_
private

Definition at line 42 of file VisiblesCutflow.h.

◆ track_pass_name_

std::string dqm::VisiblesCutflow::track_pass_name_
private

Definition at line 43 of file VisiblesCutflow.h.


The documentation for this class was generated from the following files: