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/TrigParticle.h"
6
7namespace trigger {
8NtupleWriter::NtupleWriter(const std::string& name, framework::Process& process)
9 : Producer(name, process) {}
10
11void NtupleWriter::configure(framework::config::Parameters& ps) {
12 outPath_ = ps.getParameter<std::string>("outPath");
13}
14
15// precision-limiting function
16// inline float prec(float x, unsigned int nBits=22){ return
17// float(int(x*(1<<nBits)))/(1<<nBits);}
18inline float prec(float x) { return x; }
19
20void NtupleWriter::produce(framework::Event& event) {
21 framework::NtupleManager& n{framework::NtupleManager::getInstance()};
22
23 std::string inTag;
24 inTag = "TargetScoringPlaneHits";
25 if (writeTruth_ && event.exists(inTag)) {
26 const std::vector<ldmx::SimTrackerHit> hits =
27 event.getCollection<ldmx::SimTrackerHit>(inTag);
28 ldmx::SimTrackerHit h, hMaxEle; // the desired truth hits
29 for (const auto& hit : hits) {
30 auto xyz = hit.getPosition();
31 if (xyz[2] > 0 && xyz[2] < 1) {
32 if (hit.getTrackID() == 1) h = hit;
33 if (hit.getPdgID() == 11 && (hit.getEnergy() > hMaxEle.getEnergy()))
34 hMaxEle = hit;
35 } else {
36 continue; // select one sp
37 }
38 }
39 if (h.getPdgID() == 0)
40 h = hMaxEle; // save max energy in case track1 isn't found (A')
41 std::string coll = "Truth";
42 n.setVar(coll + "_e", prec(h.getEnergy()));
43 n.setVar(coll + "_x", prec(h.getPosition()[0]));
44 n.setVar(coll + "_y", prec(h.getPosition()[1]));
45 n.setVar(coll + "_px", prec(h.getMomentum()[0]));
46 n.setVar(coll + "_py", prec(h.getMomentum()[1]));
47 n.setVar(coll + "_pz", prec(h.getMomentum()[2]));
48 n.setVar(coll + "_pdgId", h.getPdgID());
49 }
50 inTag = "EcalScoringPlaneHits";
51 if (writeTruth_ && event.exists(inTag)) {
52 const std::vector<ldmx::SimTrackerHit> hits =
53 event.getCollection<ldmx::SimTrackerHit>(inTag);
54 ldmx::SimTrackerHit h, hMaxEle; // the desired truth hits
55 for (const auto& hit : hits) {
56 auto xyz = hit.getPosition();
57 if (xyz[2] > 239.99 && xyz[2] < 240.01) {
58 if (hit.getTrackID() == 1) h = hit;
59 if (hit.getPdgID() == 11 && (hit.getEnergy() > hMaxEle.getEnergy()))
60 hMaxEle = hit;
61 } else {
62 continue; // select one sp
63 }
64 }
65 if (h.getPdgID() == 0)
66 h = hMaxEle; // save max energy in case track1 isn't found (A')
67 std::string coll = "TruthEcal";
68 n.setVar(coll + "_e", prec(h.getEnergy()));
69 n.setVar(coll + "_x", prec(h.getPosition()[0]));
70 n.setVar(coll + "_y", prec(h.getPosition()[1]));
71 n.setVar(coll + "_px", prec(h.getMomentum()[0]));
72 n.setVar(coll + "_py", prec(h.getMomentum()[1]));
73 n.setVar(coll + "_pz", prec(h.getMomentum()[2]));
74 n.setVar(coll + "_pdgId", h.getPdgID());
75 }
76
77 inTag = "ecalTrigSums";
78 if (writeEcalSums_ && event.exists(inTag)) {
79 const auto sums = event.getCollection<TrigEnergySum>(inTag);
80 // const int nEcalLayers = 34;
81 vector<float> energyAfterLayer; // (nEcalLayers, 0.);
82 for (const auto& sum : sums) {
83 if (!(sum.energy() > 0)) continue;
84 if (sum.layer() >= energyAfterLayer.size())
85 energyAfterLayer.resize(sum.layer() + 1);
86 for (int i = 0; i <= sum.layer(); i++) {
87 energyAfterLayer[i] += sum.energy();
88 }
89 }
90 n.setVar("Ecal_e_afterLayer", energyAfterLayer);
91 n.setVar("Ecal_e_nLayer", int(energyAfterLayer.size()));
92 }
93 inTag = "hcalTrigQuadsBackLayerSums";
94 if (writeHcalSums_ && event.exists(inTag)) {
95 const auto sums = event.getCollection<TrigEnergySum>(inTag);
96 vector<float> energyAfterLayer;
97 for (const auto& sum : sums) {
98 if (!(sum.hwEnergy() > 0)) continue;
99 if (sum.layer() >= energyAfterLayer.size())
100 energyAfterLayer.resize(sum.layer() + 1);
101 for (int i = 0; i <= sum.layer(); i++) {
102 energyAfterLayer[i] += sum.hwEnergy();
103 }
104 }
105 n.setVar("Hcal_e_afterLayer", energyAfterLayer);
106 n.setVar("Hcal_e_nLayer", int(energyAfterLayer.size()));
107 }
108
109 inTag = "trigElectrons";
110 if (writeEle_ && event.exists(inTag)) {
111 const auto eles = event.getCollection<TrigParticle>(inTag);
112 const int nEle = eles.size();
113 int maxE = -1;
114 float maxEVal = 0;
115 int maxPt = -1;
116 float maxPtVal = 0;
117 vector<float> v_e(nEle);
118 vector<float> v_eC(nEle);
119 vector<float> v_zC(nEle);
120 vector<float> v_px(nEle);
121 vector<float> v_py(nEle);
122 vector<float> v_pz(nEle);
123 vector<float> v_dx(nEle);
124 vector<float> v_dy(nEle);
125 vector<float> v_x(nEle);
126 vector<float> v_y(nEle);
127 vector<int> v_tp(nEle);
128 vector<int> v_depth(nEle);
129 for (unsigned int i = 0; i < nEle; i++) {
130 if (eles[i].energy() > maxEVal) {
131 maxEVal = eles[i].energy();
132 maxE = i;
133 }
134 if (eles[i].pt() > maxPtVal) {
135 maxPtVal = eles[i].pt();
136 maxPt = i;
137 }
138 v_e[i] = prec(eles[i].energy());
139 v_eC[i] = prec(eles[i].getClusEnergy());
140 v_zC[i] = prec(eles[i].endz());
141 v_px[i] = prec(eles[i].px());
142 v_py[i] = prec(eles[i].py());
143 v_pz[i] = prec(eles[i].pz());
144 v_dx[i] = prec(eles[i].endx() - eles[i].vx());
145 v_dy[i] = prec(eles[i].endy() - eles[i].vy());
146 v_x[i] = prec(eles[i].vx());
147 v_y[i] = prec(eles[i].vy());
148 v_tp[i] = prec(eles[i].getClusTP());
149 v_depth[i] = prec(eles[i].getClusDepth());
150 }
151 std::string coll = "Electron";
152 n.setVar("n" + coll, nEle);
153 n.setVar("maxE", maxE);
154 n.setVar("maxPt", maxPt);
155 n.setVar(coll + "_e", v_e);
156 n.setVar(coll + "_eClus", v_eC);
157 n.setVar(coll + "_zClus", v_zC);
158 n.setVar(coll + "_px", v_px);
159 n.setVar(coll + "_py", v_py);
160 n.setVar(coll + "_pz", v_pz);
161 n.setVar(coll + "_dx", v_dx);
162 n.setVar(coll + "_dy", v_dy);
163 n.setVar(coll + "_x", v_x);
164 n.setVar(coll + "_y", v_y);
165 n.setVar(coll + "_tp", v_tp);
166 n.setVar(coll + "_depth", v_depth);
167 }
168}
169
170void NtupleWriter::onProcessStart() {
171 // auto hdir = getHistoDirectory();
172 outFile_ = new TFile(outPath_.c_str(), "recreate");
173 outFile_->SetCompressionSettings(209);
174 // 100*alg+level
175 // 2=LZMA, 9 = max compression
176 framework::NtupleManager& n{framework::NtupleManager::getInstance()};
177 n.create(tag_);
178
179 if (writeEle_) {
180 std::string coll = "Electron";
181 n.addVar<int>(tag_, "n" + coll);
182 n.addVar<int>(tag_, "maxE");
183 n.addVar<int>(tag_, "maxPt");
184 n.addVar<vector<float> >(tag_, coll + "_e");
185 n.addVar<vector<float> >(tag_, coll + "_eClus");
186 n.addVar<vector<float> >(tag_, coll + "_zClus");
187 n.addVar<vector<float> >(tag_, coll + "_px");
188 n.addVar<vector<float> >(tag_, coll + "_py");
189 n.addVar<vector<float> >(tag_, coll + "_pz");
190 n.addVar<vector<float> >(tag_, coll + "_dx");
191 n.addVar<vector<float> >(tag_, coll + "_dy");
192 n.addVar<vector<float> >(tag_, coll + "_x"); // at target
193 n.addVar<vector<float> >(tag_, coll + "_y");
194 n.addVar<vector<int> >(tag_, coll + "_tp");
195 n.addVar<vector<int> >(tag_, coll + "_depth");
196 }
197 if (writeTruth_) {
198 n.addVar<float>(tag_, "Truth_x");
199 n.addVar<float>(tag_, "Truth_y");
200 n.addVar<float>(tag_, "Truth_px");
201 n.addVar<float>(tag_, "Truth_py");
202 n.addVar<float>(tag_, "Truth_pz");
203 n.addVar<float>(tag_, "Truth_e");
204 n.addVar<int>(tag_, "Truth_pdgId");
205 n.addVar<float>(tag_, "TruthEcal_x");
206 n.addVar<float>(tag_, "TruthEcal_y");
207 n.addVar<float>(tag_, "TruthEcal_px");
208 n.addVar<float>(tag_, "TruthEcal_py");
209 n.addVar<float>(tag_, "TruthEcal_pz");
210 n.addVar<float>(tag_, "TruthEcal_e");
211 n.addVar<int>(tag_, "TruthEcal_pdgId");
212 }
213 if (writeEcalSums_) {
214 n.addVar<vector<float> >(tag_, "Ecal_e_afterLayer");
215 n.addVar<int>(tag_, "Ecal_e_nLayer");
216 };
217 if (writeHcalSums_) {
218 n.addVar<vector<float> >(tag_, "Hcal_e_afterLayer");
219 n.addVar<int>(tag_, "Hcal_e_nLayer");
220 };
221}
222void NtupleWriter::onProcessEnd() {
223 outFile_->Write();
224 outFile_->Close();
225}
226
227} // namespace trigger
228DECLARE_PRODUCER_NS(trigger, NtupleWriter);
DECLARE_PRODUCER_NS
#define DECLARE_PRODUCER_NS(NS, CLASS)
Macro which allows the framework to construct a producer given its name during configuration.
Definition EventProcessor.h:371
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:41
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:27
framework::config::Parameters::getParameter
T getParameter(const std::string &name) const
Retrieve the parameter of the given name.
Definition Parameters.h:89
ldmx::SimTrackerHit
Represents a simulated tracker hit in the simulation.
Definition SimTrackerHit.h:24
ldmx::SimTrackerHit::getPosition
std::vector< float > getPosition() const
Get the XYZ position of the hit [mm].
Definition SimTrackerHit.h:69
trigger::TrigEnergySum
Contains the trigger output for generic calo objects.
Definition TrigEnergySum.h:24
trigger::TrigParticle
Class for particles reconstructed by the trigger system.
Definition TrigParticle.h:27
Generated by doxygen 1.9.8