Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JConvertEvt.cc File Reference

Auxiliary program to convert fit results to Evt format. More...

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Auxiliary program to convert fit results to Evt format.

Author
mdejong

Definition in file software/JReconstruction/JConvertEvt.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 66 of file software/JReconstruction/JConvertEvt.cc.

67{
68 using namespace std;
69 using namespace JPP;
70 using namespace KM3NETDAQ;
71
73 JACOUSTICS::JEvt>::typelist calibration_types;
74 typedef JMultipleFileScanner<calibration_types> JCalibration_t;
75
76 JMultipleFileScanner<> inputFile;
78 JLimit_t& numberOfEvents = inputFile.getLimit();
79 string detectorFile;
80 JCalibration_t calibrationFile;
81 double Tmax_s;
82 JPMTParametersMap pmtParameters;
84 size_t numberOfFits;
85 int debug;
86
87 try {
88
89 JParser<> zap("Auxiliary program to convert fit results to Evt format.\
90 \nThe option -L corresponds to the name of a shared library \
91 \nand function so to rearrange the order of fit results.");
92
93 zap['f'] = make_field(inputFile);
94 zap['o'] = make_field(outputFile) = "aanet.root";
95 zap['n'] = make_field(numberOfEvents) = JLimit::max();
96 zap['a'] = make_field(detectorFile) = "";
97 zap['+'] = make_field(calibrationFile) = JPARSER::initialised();
98 zap['T'] = make_field(Tmax_s) = 100.0;
99 zap['P'] = make_field(pmtParameters, "PMT simulation data (or corresponding file name)") = JPARSER::initialised();
101 zap['N'] = make_field(numberOfFits) = 0;
102 zap['d'] = make_field(debug) = 2;
103
104 zap(argc, argv);
105 }
106 catch(const exception& error) {
107 FATAL(error.what() << endl);
108 }
109
110 if (detectorFile == "" && !calibrationFile.empty()) {
111 FATAL("Missing detector file." << endl);
112 }
113
114 typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
115 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
116
118
119 if (detectorFile != "") {
120 try {
121 load(detectorFile, detector);
122 }
123 catch(const JException& error) {
124 FATAL(error);
125 }
126 }
127
128 unique_ptr<JDynamics> dynamics;
129
130 if (!calibrationFile.empty()) {
131
132 try {
133
134 dynamics.reset(new JDynamics(detector, Tmax_s));
135
136 dynamics->load(calibrationFile);
137 }
138 catch(const exception& error) {
139 FATAL(error.what());
140 }
141 }
142
143
144 const JPMTRouter pmt_router(dynamics ? dynamics->getDetector() : detector);
145 const JModuleRouter mod_router(dynamics ? dynamics->getDetector() : detector);
146
147
148 outputFile.open();
149
150 Head header;
151
152 try {
153 header = getHeader(inputFile);
154 } catch(const exception& error) {}
155
156 JAANET::JHead buffer(header);
157
158 // DAQ livetime is not yet set for real data, for Monte Carlo, it is set in JTriggerEfficiency
159
160 if (buffer.pull(&JAANET::JHead::DAQ) == buffer.end()) {
161
162 buffer.DAQ.livetime_s = getLivetime(inputFile.begin(), inputFile.end());
163 buffer.push(&JAANET::JHead::DAQ);
164
165 copy(buffer, header);
166 }
167
168 if (detectorFile != "") {
169
170 buffer.calibration.buffer = (dynamics ? calibration::dynamical() : calibration::statical());
171 buffer.push(&JAANET::JHead::calibration);
172
173 copy(buffer, header);
174 }
175
176
177 outputFile.put(header);
178
179
180 outputFile.put(JMeta(argc, argv));
181
182 map<int, size_t> miss_mod;
183 map<int, size_t> miss_pmt;
184
185 counter_type number_of_events = 0;
186
187 for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
188
189 STATUS("Processing: " << *i << endl);
190
191 JParallelFileScanner_t in(*i);
192 JTreeScanner<Evt> mc(*i);
193
194 in.setLimit(inputFile.getLimit());
195
196 Vec offset(0,0,0);
197
198 int mc_run_id = 0;
199
200 try {
201
202 const JAANET::JHead head(getHeader(*i));
203
204 offset = getOffset(head);
205 mc_run_id = head.start_run.run_id;
206
207 } catch(const exception& error) {}
208
209
210 while (in.hasNext()) {
211
212 STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
213
214 multi_pointer_type ps = in.next();
215
216 JDAQEvent* tev = ps;
217 JFIT::JEvt* evt = ps;
218
219 if (dynamics) {
220 dynamics->update(*tev); // update detector
221 }
222
223 JFIT::JEvt::iterator __end = evt->end();
224
225 if (numberOfFits > 0) {
226 advance(__end = evt->begin(), min(numberOfFits, evt->size()));
227 }
228
229 if (qualitySorter.is_valid()) {
230 partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
231 }
232
233 Evt out; // output event
234
235 if (mc.getEntries() != 0) {
236
237 Evt* event = mc.getEntry(tev->getCounter());
238
239 if (event != NULL) {
240
241 out = *event; // import Monte Carlo true data
242
243 if (out.mc_run_id == 0) {
244 out.mc_run_id = mc_run_id;
245 }
246
247 for (vector<Trk>::iterator i = out.mc_trks.begin(); i != out.mc_trks.end(); ++i) {
248 i->pos += offset; // offset true positions
249 }
250 }
251 }
252
253 read(out, *tev); // import DAQ data
254
255 if (!pmt_router->empty()) { // import calibration data
256
257 for (vector<Hit>::iterator i = out.mc_hits.begin(); i != out.mc_hits.end(); ++i) {
258
259 if (pmt_router.hasPMT(i->pmt_id)) {
260
261 const JPMTAddress address = pmt_router.getAddress(i->pmt_id);
262 const JPMTIdentifier id = pmt_router.getIdentifier(address);
263 const JPMT& pmt = pmt_router.getPMT(address);
264
265 i->dom_id = id.getID();
266 i->channel_id = id.getPMTAddress();
267 i->pos = getPosition (pmt);
268 i->dir = getDirection(pmt);
269
270 } else {
271
272 miss_pmt[i->pmt_id] += 1;
273 }
274 }
275 }
276
277 if (!mod_router->empty()) { // import calibration data
278
279 for (vector<Hit>::iterator i = out.hits.begin(); i != out.hits.end(); ++i) {
280
281 if (mod_router.hasModule(i->dom_id)) {
282
283 const JPMTIdentifier id(i->dom_id, i->channel_id);
284
285 const JPMTParameters& parameters = pmtParameters.getPMTParameters(id);
286 const JPMTAnalogueSignalProcessor cpu(parameters);
287
288 const JPMT& pmt = mod_router.getPMT(id);
289
291
292 const JTRIGGER::JHit hit(getTime(i->tdc, pmt.getCalibration()), i->tot);
293
294 i->pmt_id = pmt.getID();
295 i->pos = getPosition (pmt);
296 i->dir = getDirection(pmt);
297 i->t = hit.getT();
298 i->tot = hit.getToT();
299 i->a = cpu.getNPE(i->tot);
300
301 } else {
302
303 miss_mod[i->dom_id] += 1;
304 }
305 }
306 }
307
308 struct : public set<JUUID> {
309
310 inline int get_index(const JUUID& element) const
311 {
312 const_iterator i = this->find(element);
313
314 if (i != this->end())
315 return std::distance(this->begin(), i);
316 else
317 return -1;
318 }
319 } uuid;
320
321 for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
322 uuid.insert(fit->getUUID());
323 }
324
325 for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) { // import reconstruction data
326
327 Trk trk;
328
329 trk.id = uuid.get_index(fit->getUUID());
330 trk.pos = Vec(fit->getX(), fit->getY(), fit->getZ());
331 trk.dir = Vec(fit->getDX(), fit->getDY(), fit->getDZ());
332 trk.t = fit->getT();
333 trk.E = fit->getE();
334 trk.lik = fit->getQ();
336 trk.status = (fit->getStatus() == COMPLETE_CHAIN ? TRK_ST_FINALSTATE : TRK_ST_UNDEFINED);
337
338 if (fit->hasParentUUID()) {
339 trk.mother_id = uuid.get_index(fit->getParentUUID());
340 }
341
342 for (JHistory::const_iterator i = fit->getHistory().begin(); i != fit->getHistory().end(); ++i) {
343 trk.rec_stages.push_back(i->type);
344 }
345
346 for (int i = 0; i != fit->getN(); ++i) {
347 trk.fitinf.push_back(fit->getW(i));
348 }
349
350 trk.error_matrix = fit->getV();
351
352 out.trks.push_back(trk);
353 }
354
355 out.id = ++number_of_events;
356
357 outputFile.put(out);
358 }
359 }
360
361 STATUS(endl);
362
363 for (const auto& i : miss_pmt) { ERROR("Misses PMT " << setw(8) << i.first << ' ' << setw(8) << i.second << endl); }
364 for (const auto& i : miss_mod) { ERROR("Misses module " << setw(8) << i.first << ' ' << setw(8) << i.second << endl); }
365
367
368 io >> outputFile;
369
370 outputFile.close();
371}
string outputFile
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Monte Carlo run header.
Definition JHead.hh:1236
JAANET::DAQ DAQ
Definition JHead.hh:1607
JAANET::calibration calibration
Definition JHead.hh:1592
const JCalibration & getCalibration() const
Get calibration.
Detector data structure.
Definition JDetector.hh:96
Router for direct addressing of module data in detector data structure.
Address of PMT in detector data structure.
Auxiliary class for map of PMT parameters.
const JPMTParameters & getPMTParameters(const JPMTIdentifier &id) const
Get PMT parameters.
Data structure for PMT parameters.
bool slewing
time slewing of analogue signal
Router for direct addressing of PMT data in detector data structure.
Definition JPMTRouter.hh:37
Data structure for PMT geometry, calibration and status.
Definition JPMT.hh:49
Data structure for set of track fit results.
General exception.
Definition JException.hh:24
int getID() const
Get identifier.
Definition JObjectID.hh:50
Utility class to parse command line options.
Definition JParser.hh:1698
Object writing to file.
General purpose class for object reading from a list of file names.
General purpose class for parallel reading of objects from a single file or multiple files.
Template definition for direct access of elements in ROOT TChain.
Hit data structure.
static void setSlewing(const bool slewing)
Set slewing option.
JTriggerCounter_t getCounter() const
Get trigger counter.
static const int JPP_RECONSTRUCTION_TYPE
KM3NeT Data Definitions v3.6.0 https://git.km3net.de/common/km3net-dataformat.
JDirection3D getDirection(const Vec &dir)
Get direction.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition JHead.cc:163
JPosition3D getPosition(const Vec &pos)
Get position.
Vec getOffset(const JHead &header)
Get offset.
std::istream & read(std::istream &in, JTestSummary &summary, const char delimiter=' ')
Read test summary.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Long64_t counter_type
Type definition for counter.
double getLivetime(const std::string &file_name)
Get data taking live time.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
const char * getTime()
Get current local time conform ISO-8601 standard.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
std::vector< Hit > hits
list of hits
Definition Evt.hh:38
int mc_run_id
MC run identifier.
Definition Evt.hh:27
std::vector< Hit > mc_hits
MC: list of MC truth hits.
Definition Evt.hh:48
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
Definition Evt.hh:49
int id
offline event identifier
Definition Evt.hh:22
std::vector< Trk > trks
list of reconstructed tracks (can be several because of prefits,showers, etc).
Definition Evt.hh:39
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition Head.hh:65
Calibration.
Definition JHead.hh:330
static const std::string dynamical()
Definition JHead.hh:332
Detector file.
Definition JHead.hh:227
Acoustic event fit.
Orientation of module.
Dynamic detector calibration.
Definition JDynamics.hh:86
Auxiliary class for recursive type list generation.
Definition JTypeList.hh:351
Simple wrapper for UUID.
Definition JUUID.hh:24
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
General purpose sorter of fit results.
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
int status
MC status code, see km3net-dataformat/definitions/trkmembers.csv for values.
Definition Trk.hh:28
int id
track identifier
Definition Trk.hh:16
std::vector< double > error_matrix
(NxN) error covariance matrix for fit parameters (stored as linear vector)
Definition Trk.hh:34
std::vector< double > fitinf
place to store additional fit info, see km3net-dataformat/definitions/fitparameters....
Definition Trk.hh:32
Vec dir
track direction
Definition Trk.hh:18
std::vector< int > rec_stages
list of identifyers of succesfull fitting stages resulting in this track
Definition Trk.hh:26
double E
Energy [GeV] (either MC truth or reconstructed)
Definition Trk.hh:20
double t
track time [ns] (when the particle is at pos )
Definition Trk.hh:19
int rec_type
identifier of the fitting algorithm/chain/strategy, see km3net-dataformat/definitions/reconstruction....
Definition Trk.hh:25
int mother_id
MC id of the parent particle.
Definition Trk.hh:29
double lik
likelihood or lambda value (for aafit, lambda)
Definition Trk.hh:23
Vec pos
postion [m] of the track at time t
Definition Trk.hh:17
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition Vec.hh:13
static const int TRK_ST_FINALSTATE
for MC: the particle must be processed by detector simulation ('track_in' tag in evt files)....
Definition trkmembers.hh:15
static const int TRK_ST_UNDEFINED
status was not defined for this MC track (all reco tracks have this value)
Definition trkmembers.hh:14