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