LDMX Software
NtupleWriter.cxx
1#include "Trigger/NtupleWriter.h"
2
3#include "Framework/NtupleManager.h"
4#include "SimCore/Event/SimTrackerHit.h"
5#include "Trigger/Event/TrigEnergySum.h"
6#include "Trigger/Event/TrigMip.h"
7#include "Trigger/Event/TrigParticle.h"
8
9namespace trigger {
10NtupleWriter::NtupleWriter(const std::string& name, framework::Process& process)
11 : Producer(name, process) {}
12
13void NtupleWriter::configure(framework::config::Parameters& ps) {
14 out_path_ = ps.get<std::string>("out_path");
15
16 target_sp_hits_event_passname_ =
17 ps.get<std::string>("target_sp_hits_event_passname");
18 target_sp_passname_ = ps.get<std::string>("target_sp_passname");
19
20 ecal_sp_hits_events_passname_ =
21 ps.get<std::string>("ecal_sp_hits_events_passname");
22 ecal_sp_passname_ = ps.get<std::string>("ecal_sp_passname");
23
24 ecal_trig_sums_event_passname_ =
25 ps.get<std::string>("ecal_trig_sums_event_passname");
26 ecal_trig_sums_passname_ = ps.get<std::string>("ecal_trig_sums_passname");
27
28 trig_electrons_event_passname_ =
29 ps.get<std::string>("trig_electrons_event_passname");
30 trig_electrons_passname_ = ps.get<std::string>("trig_electrons_passname");
31
32 hcal_trig_quads_events_passname_ =
33 ps.get<std::string>("hcal_trig_quads_events_passname");
34 hcal_trig_quads_passname_ = ps.get<std::string>("hcal_trig_quads_passname");
35}
36
37// precision-limiting function
38// inline float prec(float x, unsigned int nBits=22){ return
39// float(int(x*(1<<nBits)))/(1<<nBits);}
40inline float prec(float x) { return x; }
41
42void NtupleWriter::produce(framework::Event& event) {
43 framework::NtupleManager& n{framework::NtupleManager::getInstance()};
44
45 std::string in_tag;
46 in_tag = "TargetScoringPlaneHits";
47 if (write_truth_ && event.exists(in_tag, target_sp_hits_event_passname_)) {
48 const std::vector<ldmx::SimTrackerHit> hits =
49 event.getCollection<ldmx::SimTrackerHit>(in_tag, target_sp_passname_);
50
51 ldmx::SimTrackerHit h, h_max_ele; // the desired truth hits
52 for (const auto& hit : hits) {
53 auto xyz = hit.getPosition();
54 if (xyz[2] > 0 && xyz[2] < 1) {
55 if (hit.getTrackID() == 1) h = hit;
56 if (hit.getPdgID() == 11 && (hit.getEnergy() > h_max_ele.getEnergy()))
57 h_max_ele = hit;
58 } else {
59 continue; // select one sp
60 }
61 }
62 if (h.getPdgID() == 0)
63 h = h_max_ele; // save max energy in case track1 isn't found (A')
64 std::string coll = "Truth";
65 n.setVar(coll + "_e", prec(h.getEnergy()));
66 n.setVar(coll + "_x", prec(h.getPosition()[0]));
67 n.setVar(coll + "_y", prec(h.getPosition()[1]));
68 n.setVar(coll + "_px", prec(h.getMomentum()[0]));
69 n.setVar(coll + "_py", prec(h.getMomentum()[1]));
70 n.setVar(coll + "_pz", prec(h.getMomentum()[2]));
71 n.setVar(coll + "_pdgId", h.getPdgID());
72 }
73 in_tag = "EcalScoringPlaneHits";
74 if (write_truth_ && event.exists(in_tag, ecal_sp_hits_events_passname_)) {
75 const std::vector<ldmx::SimTrackerHit> hits =
76 event.getCollection<ldmx::SimTrackerHit>(in_tag, ecal_sp_passname_);
77 ldmx::SimTrackerHit h, h_max_ele; // the desired truth hits
78 for (const auto& hit : hits) {
79 auto xyz = hit.getPosition();
80 if (xyz[2] > 239.99 && xyz[2] < 240.01) {
81 if (hit.getTrackID() == 1) h = hit;
82 if (hit.getPdgID() == 11 && (hit.getEnergy() > h_max_ele.getEnergy()))
83 h_max_ele = hit;
84 } else {
85 continue; // select one sp
86 }
87 }
88 if (h.getPdgID() == 0)
89 h = h_max_ele; // save max energy in case track1 isn't found (A')
90 std::string coll = "TruthEcal";
91 n.setVar(coll + "_e", prec(h.getEnergy()));
92 n.setVar(coll + "_x", prec(h.getPosition()[0]));
93 n.setVar(coll + "_y", prec(h.getPosition()[1]));
94 n.setVar(coll + "_px", prec(h.getMomentum()[0]));
95 n.setVar(coll + "_py", prec(h.getMomentum()[1]));
96 n.setVar(coll + "_pz", prec(h.getMomentum()[2]));
97 n.setVar(coll + "_pdgId", h.getPdgID());
98 }
99
100 in_tag = "ecalTrigSums";
101 if (write_ecal_sums_ &&
102 event.exists(in_tag, ecal_trig_sums_event_passname_)) {
103 const auto sums =
104 event.getCollection<TrigEnergySum>(in_tag, ecal_trig_sums_passname_);
105 // const int nEcalLayers = 34;
106 vector<float> energy_after_layer; // (nEcalLayers, 0.);
107 for (const auto& sum : sums) {
108 if (!(sum.energy() > 0)) continue;
109 if (sum.layer() >= energy_after_layer.size())
110 energy_after_layer.resize(sum.layer() + 1);
111 for (int i = 0; i <= sum.layer(); i++) {
112 energy_after_layer[i] += sum.energy();
113 }
114 }
115 n.setVar("Ecal_e_afterLayer", energy_after_layer);
116 n.setVar("Ecal_e_nLayer", int(energy_after_layer.size()));
117 }
118 in_tag = "hcalTrigQuadsBackLayerSums";
119 if (write_hcal_sums_ &&
120 event.exists(in_tag, hcal_trig_quads_events_passname_)) {
121 const auto sums =
122 event.getCollection<TrigEnergySum>(in_tag, hcal_trig_quads_passname_);
123 vector<float> energy_after_layer;
124 for (const auto& sum : sums) {
125 if (!(sum.hwEnergy() > 0)) continue;
126 if (sum.layer() >= energy_after_layer.size())
127 energy_after_layer.resize(sum.layer() + 1);
128 for (int i = 0; i <= sum.layer(); i++) {
129 energy_after_layer[i] += sum.hwEnergy();
130 }
131 }
132 n.setVar("Hcal_e_afterLayer", energy_after_layer);
133 n.setVar("Hcal_e_nLayer", int(energy_after_layer.size()));
134 }
135
136 in_tag = "hcalTrigQuadsSideLayerSums";
137 if (write_hcal_sums_ &&
138 event.exists(in_tag, hcal_trig_quads_events_passname_)) {
139 const auto sums =
140 event.getCollection<TrigEnergySum>(in_tag, hcal_trig_quads_passname_);
141 float energy_in_side_hcal = 0.f;
142 for (const auto& sum : sums) {
143 if (!(sum.hwEnergy() > 0)) continue;
144 for (int i = 0; i <= sum.layer(); i++) {
145 energy_in_side_hcal += sum.hwEnergy();
146 }
147 }
148 n.setVar("SideHcal_e", energy_in_side_hcal);
149 }
150
151 in_tag = "ecalTrigMIPs";
152 if (write_ecal_trig_mi_ps_ &&
153 event.exists(in_tag, ecal_trig_sums_event_passname_)) {
154 const auto mips =
155 event.getCollection<TrigMip>(in_tag, ecal_trig_sums_passname_);
156 std::vector<int> lengths;
157 std::vector<int> n_holes;
158 std::vector<float> iso_energies;
159 for (const auto& mip : mips) {
160 lengths.push_back(mip.length());
161 n_holes.push_back(mip.nHoles());
162 iso_energies.push_back(mip.sumEinIsolationRegion());
163 }
164 n.setVar("Ecal_mip_length", lengths);
165 n.setVar("Ecal_mip_nHoles", n_holes);
166 n.setVar("Ecal_mip_isolationEnergy", iso_energies);
167 }
168
169 in_tag = "hcalTrigMIPs";
170 if (write_hcal_trig_mi_ps_ &&
171 event.exists(in_tag, hcal_trig_quads_events_passname_)) {
172 const auto mips =
173 event.getCollection<TrigMip>(in_tag, hcal_trig_quads_passname_);
174 std::vector<int> lengths;
175 std::vector<int> n_holes;
176 // std::vector<float> isoEnergies;
177 for (const auto& mip : mips) {
178 lengths.push_back(mip.length());
179 n_holes.push_back(mip.nHoles());
180 // isoEnergies.push_back(mip.SumEinIsolationRegion());
181 }
182 n.setVar("Hcal_mip_length", lengths);
183 n.setVar("Hcal_mip_nHoles", n_holes);
184 // n.setVar("Hcal_mip_isolationEnergy", isoEnergies);
185 }
186
187 in_tag = "trigElectrons";
188 if (write_ele_ && event.exists(in_tag, trig_electrons_event_passname_)) {
189 const auto eles =
190 event.getCollection<TrigParticle>(in_tag, trig_electrons_passname_);
191 const int n_ele = eles.size();
192 int max_e = -1;
193 float max_e_val = 0;
194 int max_pt = -1;
195 float max_pt_val = 0;
196 vector<float> v_e(n_ele);
197 vector<float> v_e_c(n_ele);
198 vector<float> v_z_c(n_ele);
199 vector<float> v_px(n_ele);
200 vector<float> v_py(n_ele);
201 vector<float> v_pz(n_ele);
202 vector<float> v_dx(n_ele);
203 vector<float> v_dy(n_ele);
204 vector<float> v_x(n_ele);
205 vector<float> v_y(n_ele);
206 vector<int> v_tp(n_ele);
207 vector<int> v_depth(n_ele);
208 for (unsigned int i = 0; i < n_ele; i++) {
209 if (eles[i].energy() > max_e_val) {
210 max_e_val = eles[i].energy();
211 max_e = i;
212 }
213 if (eles[i].pt() > max_pt_val) {
214 max_pt_val = eles[i].pt();
215 max_pt = i;
216 }
217 v_e[i] = prec(eles[i].energy());
218 v_e_c[i] = prec(eles[i].getClusEnergy());
219 v_z_c[i] = prec(eles[i].endz());
220 v_px[i] = prec(eles[i].px());
221 v_py[i] = prec(eles[i].py());
222 v_pz[i] = prec(eles[i].pz());
223 v_dx[i] = prec(eles[i].endx() - eles[i].vx());
224 v_dy[i] = prec(eles[i].endy() - eles[i].vy());
225 v_x[i] = prec(eles[i].vx());
226 v_y[i] = prec(eles[i].vy());
227 v_tp[i] = prec(eles[i].getClusTP());
228 v_depth[i] = prec(eles[i].getClusDepth());
229 }
230 std::string coll = "Electron";
231 n.setVar("n" + coll, n_ele);
232 n.setVar("maxE", max_e);
233 n.setVar("maxPt", max_pt);
234 n.setVar(coll + "_e", v_e);
235 n.setVar(coll + "_eClus", v_e_c);
236 n.setVar(coll + "_zClus", v_z_c);
237 n.setVar(coll + "_px", v_px);
238 n.setVar(coll + "_py", v_py);
239 n.setVar(coll + "_pz", v_pz);
240 n.setVar(coll + "_dx", v_dx);
241 n.setVar(coll + "_dy", v_dy);
242 n.setVar(coll + "_x", v_x);
243 n.setVar(coll + "_y", v_y);
244 n.setVar(coll + "_tp", v_tp);
245 n.setVar(coll + "_depth", v_depth);
246 }
247}
248
249void NtupleWriter::onProcessStart() {
250 // auto hdir = getHistoDirectory();
251 out_file_ = new TFile(out_path_.c_str(), "recreate");
252 out_file_->SetCompressionSettings(209);
253 // 100*alg+level
254 // 2=LZMA, 9 = max compression
255 framework::NtupleManager& n{framework::NtupleManager::getInstance()};
256 n.create(tag_);
257
258 if (write_ele_) {
259 std::string coll = "Electron";
260 n.addVar<int>(tag_, "n" + coll);
261 n.addVar<int>(tag_, "maxE");
262 n.addVar<int>(tag_, "maxPt");
263 n.addVar<vector<float>>(tag_, coll + "_e");
264 n.addVar<vector<float>>(tag_, coll + "_eClus");
265 n.addVar<vector<float>>(tag_, coll + "_zClus");
266 n.addVar<vector<float>>(tag_, coll + "_px");
267 n.addVar<vector<float>>(tag_, coll + "_py");
268 n.addVar<vector<float>>(tag_, coll + "_pz");
269 n.addVar<vector<float>>(tag_, coll + "_dx");
270 n.addVar<vector<float>>(tag_, coll + "_dy");
271 n.addVar<vector<float>>(tag_, coll + "_x"); // at target
272 n.addVar<vector<float>>(tag_, coll + "_y");
273 n.addVar<vector<int>>(tag_, coll + "_tp");
274 n.addVar<vector<int>>(tag_, coll + "_depth");
275 }
276 if (write_truth_) {
277 n.addVar<float>(tag_, "Truth_x");
278 n.addVar<float>(tag_, "Truth_y");
279 n.addVar<float>(tag_, "Truth_px");
280 n.addVar<float>(tag_, "Truth_py");
281 n.addVar<float>(tag_, "Truth_pz");
282 n.addVar<float>(tag_, "Truth_e");
283 n.addVar<int>(tag_, "Truth_pdgId");
284 n.addVar<float>(tag_, "TruthEcal_x");
285 n.addVar<float>(tag_, "TruthEcal_y");
286 n.addVar<float>(tag_, "TruthEcal_px");
287 n.addVar<float>(tag_, "TruthEcal_py");
288 n.addVar<float>(tag_, "TruthEcal_pz");
289 n.addVar<float>(tag_, "TruthEcal_e");
290 n.addVar<int>(tag_, "TruthEcal_pdgId");
291 }
292 if (write_ecal_trig_mi_ps_) {
293 n.addVar<std::vector<int>>(tag_, "Ecal_mip_length");
294 n.addVar<std::vector<int>>(tag_, "Ecal_mip_nHoles");
295 n.addVar<std::vector<float>>(tag_, "Ecal_mip_isolationEnergy");
296 }
297 if (write_hcal_trig_mi_ps_) {
298 n.addVar<std::vector<int>>(tag_, "Hcal_mip_length");
299 n.addVar<std::vector<int>>(tag_, "Hcal_mip_nHoles");
300 // n.addVar<std::vector<float>>(tag_, "Hcal_mip_isolationEnergy");
301 }
302 if (write_ecal_sums_) {
303 n.addVar<vector<float>>(tag_, "Ecal_e_afterLayer");
304 n.addVar<int>(tag_, "Ecal_e_nLayer");
305 };
306 if (write_hcal_sums_) {
307 n.addVar<vector<float>>(tag_, "Hcal_e_afterLayer");
308 n.addVar<int>(tag_, "Hcal_e_nLayer");
309 n.addVar<float>(tag_, "SideHcal_e");
310 };
311}
312void NtupleWriter::onProcessEnd() {
313 out_file_->Write();
314 out_file_->Close();
315}
316
317} // namespace trigger
318DECLARE_PRODUCER(trigger::NtupleWriter);
DECLARE_PRODUCER
#define DECLARE_PRODUCER(CLASS)
Macro which allows the framework to construct a producer given its name during configuration.
Definition EventProcessor.h:346
SimTrackerHit.h
Class which encapsulates information from a hit in a simulated tracking detector.
framework::Event
Implements an event buffer system for storing event data.
Definition Event.h:40
framework::Event::exists
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:107
framework::NtupleManager
Singleton class used to manage the creation and pooling of ntuples.
Definition NtupleManager.h:35
framework::NtupleManager::create
void create(const std::string &tname)
Create a ROOT tree to hold the ntuple variables (ROOT leaves).
Definition NtupleManager.cxx:14
framework::NtupleManager::getInstance
static NtupleManager & getInstance()
Definition NtupleManager.cxx:8
framework::Process
Class which represents the process under execution.
Definition Process.h:34
framework::config::Parameters
Class encapsulating parameters for configuring a processor.
Definition Parameters.h:26
framework::config::Parameters::get
const T & get(const std::string &name) const
Retrieve the parameter of the given name.
Definition Parameters.h:75
ldmx::SimTrackerHit
Represents a simulated tracker hit in the simulation.
Definition SimTrackerHit.h:25
trigger::NtupleWriter
Null algorithm test.
Definition NtupleWriter.h:15
trigger::TrigEnergySum
Contains the trigger output for generic calo objects.
Definition TrigEnergySum.h:24
trigger::TrigMip
Class for clusters built from trigger calo hits.
Definition TrigMip.h:11
trigger::TrigParticle
Class for particles reconstructed by the trigger system.
Definition TrigParticle.h:26
Generated by doxygen 1.12.0