Jpp  18.5.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  try {
131 
132  dynamics.reset(new JDynamics(detector, Tmax_s));
133 
134  dynamics->load(calibrationFile);
135  }
136  catch(const exception& error) {
137  if (!calibrationFile.empty()) {
138  FATAL(error.what());
139  }
140  }
141 
142 
143  const JPMTRouter pmt_router(dynamics ? dynamics->getDetector() : detector);
144  const JModuleRouter mod_router(dynamics ? dynamics->getDetector() : detector);
145 
146 
147  outputFile.open();
148 
149  Head header;
150 
151  try {
152  header = getHeader(inputFile);
153  } catch(const exception& error) {}
154 
155  JAANET::JHead buffer(header);
156 
157  // DAQ livetime is not yet set for real data, for Monte Carlo, it is set in JTriggerEfficiency
158 
159  if (buffer.pull(&JAANET::JHead::DAQ) == buffer.end()) {
160 
161  buffer.DAQ.livetime_s = getLivetime(inputFile.begin(), inputFile.end());
162  buffer.push(&JAANET::JHead::DAQ);
163 
164  copy(buffer, header);
165  }
166 
167  if (detectorFile != "") {
168 
169  buffer.calibration.buffer = (dynamics ? calibration::dynamical() : calibration::statical());
171 
172  copy(buffer, header);
173  }
174 
175 
176  outputFile.put(header);
177 
178 
179  outputFile.put(JMeta(argc, argv));
180 
181  map<int, size_t> miss_mod;
182  map<int, size_t> miss_pmt;
183 
184  counter_type number_of_events = 0;
185 
186  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
187 
188  STATUS("Processing: " << *i << endl);
189 
190  JParallelFileScanner_t in(*i);
191  JTreeScanner<Evt> mc(*i);
192 
193  in.setLimit(inputFile.getLimit());
194 
195  Vec center(0,0,0);
196 
197  int mc_run_id = 0;
198 
199  try {
200 
201  const JAANET::JHead head(getHeader(*i));
202 
203  center = get<Vec>(head);
204  mc_run_id = head.start_run.run_id;
205 
206  } catch(const exception& error) {}
207 
208 
209  while (in.hasNext()) {
210 
211  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
212 
213  multi_pointer_type ps = in.next();
214 
215  JDAQEvent* tev = ps;
216  JFIT::JEvt* evt = ps;
217 
218  if (dynamics) {
219  dynamics->update(*tev); // update detector
220  }
221 
222  JFIT::JEvt::iterator __end = evt->end();
223 
224  if (numberOfFits > 0) {
225  advance(__end = evt->begin(), min(numberOfFits, evt->size()));
226  }
227 
228  if (qualitySorter.is_valid()) {
229  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
230  }
231 
232  Evt out; // output event
233 
234  if (mc.getEntries() != 0) {
235 
236  Evt* event = mc.getEntry(tev->getCounter());
237 
238  if (event != NULL) {
239 
240  out = *event; // import Monte Carlo true data
241 
242  if (out.mc_run_id == 0) {
243  out.mc_run_id = mc_run_id;
244  }
245 
246  for (vector<Trk>::iterator i = out.mc_trks.begin(); i != out.mc_trks.end(); ++i) {
247  i->pos += center; // offset true positions
248  }
249  }
250  }
251 
252  read(out, *tev); // import DAQ data
253 
254  if (!pmt_router->empty()) { // import calibration data
255 
256  for (vector<Hit>::iterator i = out.mc_hits.begin(); i != out.mc_hits.end(); ++i) {
257 
258  if (pmt_router.hasPMT(i->pmt_id)) {
259 
260  const JPMTAddress address = pmt_router.getAddress(i->pmt_id);
261  const JPMTIdentifier id = pmt_router.getIdentifier(address);
262  const JPMT& pmt = pmt_router.getPMT(address);
263 
264  i->dom_id = id.getID();
265  i->channel_id = id.getPMTAddress();
266  i->pos = getPosition (pmt);
267  i->dir = getDirection(pmt);
268 
269  } else {
270 
271  miss_pmt[i->pmt_id] += 1;
272  }
273  }
274  }
275 
276  if (!mod_router->empty()) { // import calibration data
277 
278  for (vector<Hit>::iterator i = out.hits.begin(); i != out.hits.end(); ++i) {
279 
280  if (mod_router.hasModule(i->dom_id)) {
281 
282  const JPMTIdentifier id(i->dom_id, i->channel_id);
283 
284  const JPMTParameters& parameters = pmtParameters.getPMTParameters(id);
285  const JPMTAnalogueSignalProcessor cpu(parameters);
286 
287  const JPMT& pmt = mod_router.getPMT(id);
288 
290 
291  const JTRIGGER::JHit hit(getTime(i->tdc, pmt.getCalibration()), i->tot);
292 
293  i->pmt_id = pmt.getID();
294  i->pos = getPosition (pmt);
295  i->dir = getDirection(pmt);
296  i->t = hit.getT();
297  i->tot = hit.getToT();
298  i->a = cpu.getNPE(i->tot);
299 
300  } else {
301 
302  miss_mod[i->dom_id] += 1;
303  }
304  }
305  }
306 
307  struct : public set<JUUID> {
308 
309  inline int get_index(const JUUID& element) const
310  {
311  const_iterator i = this->find(element);
312 
313  if (i != this->end())
314  return std::distance(this->begin(), i);
315  else
316  return -1;
317  }
318  } uuid;
319 
320  for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
321  uuid.insert(fit->getUUID());
322  }
323 
324  for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) { // import reconstruction data
325 
326  Trk trk;
327 
328  trk.id = uuid.get_index(fit->getUUID());
329  trk.pos = Vec(fit->getX(), fit->getY(), fit->getZ());
330  trk.dir = Vec(fit->getDX(), fit->getDY(), fit->getDZ());
331  trk.t = fit->getT();
332  trk.E = fit->getE();
333  trk.lik = fit->getQ();
335  trk.status = (fit->getStatus() == COMPLETE_CHAIN ? TRK_ST_FINALSTATE : TRK_ST_UNDEFINED);
336 
337  if (fit->hasParentUUID()) {
338  trk.mother_id = uuid.get_index(fit->getParentUUID());
339  }
340 
341  for (JHistory::const_iterator i = fit->getHistory().begin(); i != fit->getHistory().end(); ++i) {
342  trk.rec_stages.push_back(i->type);
343  }
344 
345  for (int i = 0; i != fit->getN(); ++i) {
346  trk.fitinf.push_back(fit->getW(i));
347  }
348 
349  trk.error_matrix = fit->getV();
350 
351  out.trks.push_back(trk);
352  }
353 
354  out.id = ++number_of_events;
355 
356  outputFile.put(out);
357  }
358  }
359 
360  STATUS(endl);
361 
362  for (const auto& i : miss_pmt) { ERROR("Misses PMT " << setw(8) << i.first << ' ' << setw(8) << i.second << endl); }
363  for (const auto& i : miss_mod) { ERROR("Misses module " << setw(8) << i.first << ' ' << setw(8) << i.second << endl); }
364 
366 
367  io >> outputFile;
368 
369  outputFile.close();
370 }
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:35
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
double getT() const
Get calibrated time of hit.
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:24
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
std::istream & read(std::istream &in, JTestSummary &summary, const char delimiter= ' ')
Read test summary.
const_iterator pull(T JHead::*pd) const
Pull given data member from Head.
Definition: JHead.hh:1349
static void setSlewing(const bool slewing)
Set slewing option.
double t
track time [ns] (when the particle is at pos )
Definition: Trk.hh:19
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
const JCalibration & getCalibration() const
Get calibration.
#define STATUS(A)
Definition: JMessage.hh:63
ROOT TTree parameter settings.
Detector data structure.
Definition: JDetector.hh:89
Recording of objects on file according a format that follows from the file name extension.
Router for direct addressing of module data in detector data structure.
Vec dir
track direction
Definition: Trk.hh:18
General purpose sorter of fit results.
JAANET::start_run start_run
Definition: JHead.hh:1580
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
Hit data structure.
General purpose class for parallel reading of objects from a single file or multiple files...
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
Long64_t counter_type
Type definition for counter.
double livetime_s
Live time [s].
Definition: JHead.hh:1053
double getTime(const Hit &hit)
Get true time of hit.
Basic data structure for time and time over threshold information of hit.
string outputFile
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
Data structure for detector geometry and calibration.
Tools for handling different hit types.
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
ROOT TTree parameter settings.
JAANET::calibration calibration
Definition: JHead.hh:1589
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
int mother_id
MC id of the parent particle.
Definition: Trk.hh:29
double getToT() const
Get calibrated time over threshold of hit.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
bool slewing
time slewing of analogue signal
std::string buffer
General purpose name.
Definition: JHead.hh:217
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:12
Detector file.
Definition: JHead.hh:226
JDirection3D getDirection(const Vec &dir)
Get direction.
Acoustic event fit.
int mc_run_id
MC run identifier.
Definition: Evt.hh:27
Monte Carlo run header.
Definition: JHead.hh:1234
Auxiliary class for recursive type list generation.
Definition: JTypeList.hh:351
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
Auxiliary methods for handling file names, type names and environment.
int getID() const
Get identifier.
Definition: JObjectID.hh:50
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:43
ROOT I/O of application specific meta data.
JPosition3D getPosition(const Vec &pos)
Get position.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
std::vector< double > fitinf
place to store additional fit info, see km3net-dataformat/definitions/fitparameters.csv
Definition: Trk.hh:32
#define ERROR(A)
Definition: JMessage.hh:66
Support methods.
virtual double getNPE(const double tot_ns) const override
Get number of photo-electrons.
Auxiliary class for map of PMT parameters.
Dynamic detector calibration.
Direct access to PMT in detector data structure.
static const int TRK_ST_UNDEFINED
status was not defined for this MC track (all reco tracks have this value)
Definition: trkmembers.hh:14
General purpose messaging.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
int run_id
MC run number.
Definition: JHead.hh:143
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:65
void push(T JHead::*pd)
Push given data member to Head.
Definition: JHead.hh:1374
#define FATAL(A)
Definition: JMessage.hh:67
int id
track identifier
Definition: Trk.hh:16
int status
MC status code, see km3net-dataformat/definitions/trkmembers.csv for values.
Definition: Trk.hh:28
Direct access to module in detector data structure.
Dynamic detector calibration.
Definition: JDynamics.hh:81
static const int TRK_ST_FINALSTATE
for MC: the particle must be processed by detector simulation (&#39;track_in&#39; tag in evt files)...
Definition: trkmembers.hh:15
Vec pos
postion [m] of the track at time t
Definition: Trk.hh:17
Address of PMT in detector data structure.
Definition: JPMTAddress.hh:32
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Data structure for set of track fit results.
std::vector< Trk > trks
list of reconstructed tracks (can be several because of prefits,showers, etc).
Definition: Evt.hh:39
double getLivetime(const std::string &file_name)
Get data taking live time.
std::vector< int > rec_stages
list of identifyers of succesfull fitting stages resulting in this track
Definition: Trk.hh:26
General purpose class for object reading from a list of file names.
std::vector< double > error_matrix
(NxN) error covariance matrix for fit parameters (stored as linear vector)
Definition: Trk.hh:34
double lik
likelihood or lambda value (for aafit, lambda)
Definition: Trk.hh:23
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:48
Utility class to parse command line options.
std::vector< Hit > mc_hits
MC: list of MC truth hits.
Definition: Evt.hh:48
static const int JPP_RECONSTRUCTION_TYPE
KM3NeT Data Definitions v3.3.0 https://git.km3net.de/common/km3net-dataformat.
Simple wrapper for UUID.
Definition: JUUID.hh:22
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:162
PMT analogue signal processor.
int id
offline event identifier
Definition: Evt.hh:22
Orientation of module.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
Compass event data types.
Data structure for PMT parameters.
std::vector< Hit > hits
list of hits
Definition: Evt.hh:38
do set_variable DETECTOR_TXT $WORKDIR detector
Acoustic event fit.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
JAANET::DAQ DAQ
Definition: JHead.hh:1604
JTriggerCounter_t getCounter() const
Get trigger counter.
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
Definition: Evt.hh:49
int debug
debug level
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
int rec_type
identifier of the fitting algorithm/chain/strategy, see km3net-dataformat/definitions/reconstruction...
Definition: Trk.hh:25