Jpp  debug
the software that should make you happy
software/JReconstruction/JConvertEvt.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <algorithm>
6 #include <memory>
7 #include <set>
8 #include <map>
9 
10 #include "TROOT.h"
11 #include "TFile.h"
12 
18 
19 #include "JDetector/JDetector.hh"
21 #include "JDetector/JPMTRouter.hh"
25 
26 #include "JDynamics/JDynamics.hh"
27 
28 #include "JCompass/JEvt.hh"
29 #include "JCompass/JSupport.hh"
30 
31 #include "JAcoustics/JEvt.hh"
32 #include "JAcoustics/JSupport.hh"
33 
34 #include "JTrigger/JHit.hh"
35 #include "JTrigger/JHitToolkit.hh"
36 
37 #include "JAAnet/JHead.hh"
38 #include "JAAnet/JHeadToolkit.hh"
39 #include "JAAnet/JAAnetToolkit.hh"
40 
41 #include "JDAQ/JDAQEventIO.hh"
43 
45 #include "JSupport/JTreeScanner.hh"
49 #include "JSupport/JSupport.hh"
50 #include "JSupport/JMeta.hh"
51 
52 #include "JReconstruction/JEvt.hh"
54 
55 #include "Jeep/JeepToolkit.hh"
56 #include "Jeep/JParser.hh"
57 #include "Jeep/JMessage.hh"
58 
59 
60 /**
61  * \file
62  * Auxiliary program to convert fit results to Evt format.
63  *
64  * \author mdejong
65  */
66 int main(int argc, char **argv)
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());
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 }
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
Acoustic event fit.
ROOT TTree parameter settings.
Compass event data types.
ROOT TTree parameter settings.
string outputFile
Data structure for detector geometry and calibration.
Dynamic detector calibration.
Recording of objects on file according a format that follows from the file name extension.
Tools for handling different hit types.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define STATUS(A)
Definition: JMessage.hh:63
#define ERROR(A)
Definition: JMessage.hh:66
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
ROOT I/O of application specific meta data.
Direct access to module in detector data structure.
PMT analogue signal processor.
Direct access to PMT in detector data structure.
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
Support methods.
ROOT TTree parameter settings of various packages.
Basic data structure for time and time over threshold information of hit.
Auxiliary methods for handling file names, type names and environment.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Monte Carlo run header.
Definition: JHead.hh:1236
JAANET::start_run start_run
Definition: JHead.hh:1583
JAANET::DAQ DAQ
Definition: JHead.hh:1607
void push(T JHead::*pd)
Push given data member to Head.
Definition: JHead.hh:1374
const_iterator pull(T JHead::*pd) const
Pull given data member from Head.
Definition: JHead.hh:1349
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.
bool hasModule(const JObjectID &id) const
Has module.
const JPMT & getPMT(const JPMTIdentifier &id) const
Get PMT parameters.
Address of PMT in detector data structure.
Definition: JPMTAddress.hh:35
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
bool hasPMT(const JObjectID &id) const
Has PMT.
Definition: JPMTRouter.hh:116
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Definition: JPMTRouter.hh:92
JPMTIdentifier getIdentifier(const JPMTAddress &address) const
Get identifier of PMT.
Definition: JPMTRouter.hh:128
const JPMTAddress & getAddress(const JObjectID &id) const
Get address of PMT.
Definition: JPMTRouter.hh:80
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:1714
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.
Hit data structure.
static void setSlewing(const bool slewing)
Set slewing option.
double getToT() const
Get calibrated time over threshold of hit.
double getT() const
Get calibrated time of hit.
JTriggerCounter_t getCounter() const
Get trigger counter.
static const int JPP_RECONSTRUCTION_TYPE
KM3NeT Data Definitions v3.4.0-8-ge14cb17 https://git.km3net.de/common/km3net-dataformat.
JDirection3D getDirection(const Vec &dir)
Get direction.
double getTime(const Hit &hit)
Get true time of hit.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:162
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.
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
Definition: JSTDTypes.hh:14
int main(int argc, char **argv)
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
double livetime_s
Live time [s].
Definition: JHead.hh:1053
std::string buffer
General purpose name.
Definition: JHead.hh:217
Detector file.
Definition: JHead.hh:227
int run_id
MC run number.
Definition: JHead.hh:143
Acoustic event fit.
Orientation of module.
virtual double getNPE(const double tot_ns) const override
Get number of photo-electrons.
Dynamic detector calibration.
Definition: JDynamics.hh:84
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:84
General purpose sorter of fit results.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:45
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