LDMX Software
dqm::VisiblesFeatureProducer Class Reference

Public Member Functions

 VisiblesFeatureProducer (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

bool training_ {false}
 
std::string training_file_
 
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_
 

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 VisiblesFeatureProducer.h.

Constructor & Destructor Documentation

◆ VisiblesFeatureProducer()

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

Definition at line 19 of file VisiblesFeatureProducer.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::VisiblesFeatureProducer::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 53 of file VisiblesFeatureProducer.cxx.

53 {
54 std::vector<double> bdt_features;
55
56 const auto &particle_map{event.getMap<int, ldmx::SimParticle>(
57 sim_particles_coll_name_, sim_particles_pass_name_)};
58
59 // Get target scoring plane hits for recoil electron
60 // Use this to calculate the projected photon line vector
61 // This currently uses truth-level information, but it should be replaced
62 // by reconstructed tracker information, when available
63 std::vector<double> gamma_p(3);
64 std::vector<double> gamma_x0(3);
65
66 std::vector<double> recoil_p(3);
67 bool found_recoil_e = false;
68
69 if (recoil_from_tracking_) {
70 const auto &recoil_tracks{
71 event.getCollection<ldmx::Track>(track_collection_, track_pass_name_)};
72 // Fill this in later when you know how to use it
73 for (auto &track : recoil_tracks) {
74 // need to figure out how to best isolate candidate electron track
75 auto trk_pos = track.getPositionAtTarget();
76 auto trk_mom = track.getMomentumAtTarget();
77 if (track.getCharge() == 1 && track.getNhits() == 5 &&
78 trk_pos.size() == 3 && trk_mom.size() == 3) {
79 gamma_x0 = trk_pos;
80 gamma_p[0] = -1. * trk_mom[0];
81 gamma_p[1] = -1. * trk_mom[1];
82 gamma_p[2] = beam_energy_mev_ - trk_mom[2];
83 }
84 }
85 } else {
86 if (event.exists(sp_collection_, sp_pass_name_)) {
87 const auto &target_sp_hits = event.getCollection<ldmx::SimTrackerHit>(
88 sp_collection_, sp_pass_name_);
89 bool found_rec = false;
90 for (auto const &it : particle_map) {
91 for (auto const &sphit : target_sp_hits) {
92 if (sphit.getPosition()[2] > 0) {
93 if (it.first == sphit.getTrackID()) {
94 if (it.second.getPdgID() == 622) {
95 std::vector<float> x0_float = sphit.getPosition();
96 std::vector<double> x0_double(x0_float.begin(), x0_float.end());
97 gamma_x0 = x0_double;
98 gamma_p = sphit.getMomentum();
99 found_rec = true;
100 }
101 if (it.second.getPdgID() == 11 &&
102 inList(it.second.getParents(), 0)) {
103 if (!found_rec) {
104 std::vector<float> x0_float = sphit.getPosition();
105 std::vector<double> x0_double(x0_float.begin(),
106 x0_float.end());
107 gamma_x0 = x0_double;
108 gamma_p[0] = -1. * sphit.getMomentum()[0];
109 gamma_p[1] = -1. * sphit.getMomentum()[1];
110 gamma_p[2] = beam_energy_mev_ - sphit.getMomentum()[2];
111 found_rec = true;
112 }
113 recoil_p = sphit.getMomentum();
114 found_recoil_e = true;
115 }
116 }
117 }
118 }
119 }
120 }
121 }
122
123 double p_mag = 0.;
124 if (found_recoil_e) {
125 p_mag = std::sqrt(recoil_p[0] * recoil_p[0] + recoil_p[1] * recoil_p[1] +
126 recoil_p[2] * recoil_p[2]);
127 }
128
129 const auto &ecal_rec_hits = event.getCollection<ldmx::EcalHit>(
130 ecal_rec_collection_, ecal_rec_pass_name_);
131 const auto &hcal_rec_hits = event.getCollection<ldmx::HcalHit>(
132 hcal_rec_collection_, hcal_rec_pass_name_);
133
134 double ecal_energy = 0.;
135 double hcal_energy = 0.;
136 bool hcal_containment = true;
137
138 for (const ldmx::EcalHit &hit : ecal_rec_hits) {
139 if (hit.getEnergy() > 0.) {
140 ecal_energy += hit.getEnergy();
141 }
142 }
143 for (const ldmx::HcalHit &hit : hcal_rec_hits) {
144 if (hit.getEnergy() > 0.) {
145 ldmx::HcalID det_id(hit.getID());
146 if (det_id.getSection() != 0) {
147 continue;
148 }
149 if (det_id.getLayerID() == 1 && hit.getPE() > 5) {
150 hcal_containment = false;
151 }
152 hcal_energy += 12. * hit.getEnergy();
153 }
154 }
155
156 if (ecal_energy < 3160 && hcal_energy > 4840 && hcal_containment &&
157 p_mag < 2400) {
158 // initialize all of the features
159 int n_layers_hit = 0;
160 double x_std = 0.;
161 double y_std = 0.;
162 double z_std = 0.;
163 double x_mean = 0.;
164 double y_mean = 0.;
165 double r_mean = 0.;
166 int iso_hits = 0;
167 double iso_energy = 0.;
168 int n_readout_hits = 0;
169 double summed_det = 0.;
170 double r_mean_from_photon_track = 0.;
171
172 double z_mean = 0.; // need this when calculating z_std
173 std::vector<int> layers_hit;
174
175 for (const ldmx::HcalHit &hit : hcal_rec_hits) {
176 if (hit.getEnergy() > 0.) {
177 ldmx::HcalID det_id(hit.getID());
178 if (det_id.getSection() != 0) { // skip hits that aren't in main Hcal
179 continue;
180 }
181 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
182 continue;
183 }
184
185 n_readout_hits += 1;
186 double hit_x = hit.getXPos();
187 double hit_y = hit.getYPos();
188 double hit_z = hit.getZPos();
189 double hit_r = sqrt(hit_x * hit_x + hit_y * hit_y);
190
191 summed_det += hit.getEnergy();
192
193 x_mean += hit_x * hit.getEnergy();
194 y_mean += hit_y * hit.getEnergy();
195 z_mean += hit_z * hit.getEnergy();
196 r_mean += hit_r * hit.getEnergy();
197
198 // check if this is a new layer in the collection
199 if (!(std::find(layers_hit.begin(), layers_hit.end(),
200 det_id.getLayerID()) != layers_hit.end())) {
201 layers_hit.push_back(det_id.getLayerID());
202 }
203
204 double proj_x =
205 gamma_x0[0] + (hit_z - gamma_x0[2]) * gamma_p[0] / gamma_p[2];
206 double proj_y =
207 gamma_x0[1] + (hit_z - gamma_x0[2]) * gamma_p[1] / gamma_p[2];
208
209 r_mean_from_photon_track +=
210 hit.getEnergy() * sqrt((hit_x - proj_x) * (hit_x - proj_x) +
211 (hit_y - proj_y) * (hit_y - proj_y));
212
213 // Calculate isolated hits
214 double closest_point = 9999.;
215 for (const ldmx::HcalHit &hit2 : hcal_rec_hits) {
216 if (hit2.getEnergy() > 0.) {
217 ldmx::HcalID det_i_d2(hit2.getID());
218 if (fabs(hit2.getXPos()) > 1000 || fabs(hit2.getYPos()) > 1000) {
219 continue;
220 }
221 if (det_i_d2.getLayerID() == det_id.getLayerID()) {
222 // Determine if a bar is vertical (along y-axis) or horizontal
223 // (along x-axis) Odd layers have horizontal strips Even layers
224 // have vertical strips
225 if (hit.isOrientationX()) {
226 if (fabs(hit2.getYPos() - hit_y) > 0) {
227 if (fabs(hit2.getYPos() - hit_y) < closest_point) {
228 closest_point = fabs(hit2.getYPos() - hit_y);
229 }
230 }
231 }
232 if (hit.isOrientationY()) {
233 if (fabs(hit2.getXPos() - hit_x) > 0) {
234 if (fabs(hit2.getXPos() - hit_x) < closest_point) {
235 closest_point = fabs(hit2.getXPos() - hit_x);
236 }
237 }
238 }
239 }
240 }
241 }
242 if (closest_point > 50.) {
243 iso_hits += 1;
244 iso_energy += hit.getEnergy();
245 }
246 }
247 }
248
249 n_layers_hit = layers_hit.size();
250
251 if (summed_det > 0.) {
252 x_mean /= summed_det;
253 y_mean /= summed_det;
254 z_mean /= summed_det;
255 r_mean /= summed_det;
256
257 r_mean_from_photon_track /= summed_det;
258 }
259
260 for (const ldmx::HcalHit &hit : hcal_rec_hits) {
261 if (hit.getEnergy() > 0.) {
262 if (fabs(hit.getXPos()) > 1000 || fabs(hit.getYPos()) > 1000) {
263 continue;
264 }
265 ldmx::HcalID det_id(hit.getID());
266 if (det_id.getSection() == 0) {
267 x_std += hit.getEnergy() * (hit.getXPos() - x_mean) *
268 (hit.getXPos() - x_mean);
269 y_std += hit.getEnergy() * (hit.getYPos() - y_mean) *
270 (hit.getYPos() - y_mean);
271 z_std += hit.getEnergy() * (hit.getZPos() - z_mean) *
272 (hit.getZPos() - z_mean);
273 }
274 }
275 }
276
277 if (summed_det > 0.) {
278 x_std = sqrt(x_std / summed_det);
279 y_std = sqrt(y_std / summed_det);
280 z_std = sqrt(z_std / summed_det);
281 }
282
283 // Fill histograms
284 histograms_.fill("layers_hit", n_layers_hit);
285 histograms_.fill("x_std", x_std);
286 histograms_.fill("y_std", y_std);
287 histograms_.fill("z_std", z_std);
288 histograms_.fill("x_mean", x_mean);
289 histograms_.fill("y_mean", y_mean);
290 histograms_.fill("r_mean", r_mean);
291 histograms_.fill("iso_hits", iso_hits);
292 histograms_.fill("iso_energy", iso_energy);
293 histograms_.fill("n_hits", n_readout_hits);
294 histograms_.fill("total_energy", hcal_energy);
295 histograms_.fill("photon_track", r_mean_from_photon_track);
296
297 bdt_features.push_back(n_layers_hit);
298 bdt_features.push_back(x_std);
299 bdt_features.push_back(y_std);
300 bdt_features.push_back(z_std);
301 bdt_features.push_back(x_mean);
302 bdt_features.push_back(y_mean);
303 bdt_features.push_back(r_mean);
304 bdt_features.push_back(iso_hits);
305 bdt_features.push_back(iso_energy);
306 bdt_features.push_back(n_readout_hits);
307 bdt_features.push_back(hcal_energy);
308 bdt_features.push_back(r_mean_from_photon_track);
309
310 if (training_) {
311 std::ofstream file(training_file_, std::ios::app);
312 if (!file.is_open()) {
313 ldmx_log(fatal) << "Error: Could not open file " << training_file_;
314 return;
315 }
316 for (int i = 0; i < bdt_features.size(); ++i) {
317 file << bdt_features[i] << (i + 1 == bdt_features.size() ? "\n" : ", ");
318 }
319 }
320 }
321
322 return;
323}
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:105
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:24
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::VisiblesFeatureProducer::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 22 of file VisiblesFeatureProducer.cxx.

23 {
24 training_ = parameters.get<bool>("training");
25 training_file_ = parameters.get<std::string>("training_file");
26
27 beam_energy_mev_ = parameters.get<double>("beam_energy");
28
29 // collection names
30 hcal_rec_collection_ = parameters.get<std::string>("hcal_rec_coll_name");
31 hcal_rec_pass_name_ = parameters.get<std::string>("hcal_rec_pass_name");
32
33 ecal_rec_collection_ = parameters.get<std::string>("ecal_rec_coll_name");
34 ecal_rec_pass_name_ = parameters.get<std::string>("ecal_rec_pass_name");
35
36 recoil_from_tracking_ = parameters.get<bool>("recoil_from_tracking");
37 track_collection_ = parameters.get<std::string>("track_collection");
38 track_pass_name_ = parameters.get<std::string>("track_pass_name");
39
40 sp_collection_ = parameters.get<std::string>("sp_coll_name");
41 sp_pass_name_ = parameters.get<std::string>("sp_pass_name");
42
43 sim_particles_coll_name_ =
44 parameters.get<std::string>("sim_particles_coll_name");
45 sim_particles_pass_name_ =
46 parameters.get<std::string>("sim_particles_pass_name");
47}
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::VisiblesFeatureProducer::inList ( std::vector< int > parents,
int track_id )
private

Definition at line 49 of file VisiblesFeatureProducer.cxx.

49 {
50 return std::find(parents.begin(), parents.end(), track_id) != parents.end();
51}

Member Data Documentation

◆ beam_energy_mev_

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

Definition at line 32 of file VisiblesFeatureProducer.h.

32{0.};

◆ ecal_rec_collection_

std::string dqm::VisiblesFeatureProducer::ecal_rec_collection_
private

Definition at line 36 of file VisiblesFeatureProducer.h.

◆ ecal_rec_pass_name_

std::string dqm::VisiblesFeatureProducer::ecal_rec_pass_name_
private

Definition at line 37 of file VisiblesFeatureProducer.h.

◆ hcal_rec_collection_

std::string dqm::VisiblesFeatureProducer::hcal_rec_collection_
private

Definition at line 34 of file VisiblesFeatureProducer.h.

◆ hcal_rec_pass_name_

std::string dqm::VisiblesFeatureProducer::hcal_rec_pass_name_
private

Definition at line 35 of file VisiblesFeatureProducer.h.

◆ recoil_from_tracking_

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

Definition at line 38 of file VisiblesFeatureProducer.h.

38{false};

◆ sim_particles_coll_name_

std::string dqm::VisiblesFeatureProducer::sim_particles_coll_name_
private

Definition at line 43 of file VisiblesFeatureProducer.h.

◆ sim_particles_pass_name_

std::string dqm::VisiblesFeatureProducer::sim_particles_pass_name_
private

Definition at line 44 of file VisiblesFeatureProducer.h.

◆ sp_collection_

std::string dqm::VisiblesFeatureProducer::sp_collection_
private

Definition at line 41 of file VisiblesFeatureProducer.h.

◆ sp_pass_name_

std::string dqm::VisiblesFeatureProducer::sp_pass_name_
private

Definition at line 42 of file VisiblesFeatureProducer.h.

◆ track_collection_

std::string dqm::VisiblesFeatureProducer::track_collection_
private

Definition at line 39 of file VisiblesFeatureProducer.h.

◆ track_pass_name_

std::string dqm::VisiblesFeatureProducer::track_pass_name_
private

Definition at line 40 of file VisiblesFeatureProducer.h.

◆ training_

bool dqm::VisiblesFeatureProducer::training_ {false}
private

Definition at line 29 of file VisiblesFeatureProducer.h.

29{false};

◆ training_file_

std::string dqm::VisiblesFeatureProducer::training_file_
private

Definition at line 30 of file VisiblesFeatureProducer.h.


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