Jpp  18.3.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
software/JReconstruction/JConvertEvt.cc File Reference

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

#include <string>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <memory>
#include <set>
#include "TROOT.h"
#include "TFile.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/io_online.hh"
#include "km3net-dataformat/definitions/trkmembers.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JDynamics/JDynamics.hh"
#include "JCompass/JEvt.hh"
#include "JCompass/JSupport.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JSupport.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JHitToolkit.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupportToolkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "Jeep/JeepToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

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

int main ( int  argc,
char **  argv 
)

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

66 {
67  using namespace std;
68  using namespace JPP;
69  using namespace KM3NETDAQ;
70 
72  JACOUSTICS::JEvt>::typelist calibration_types;
73  typedef JMultipleFileScanner<calibration_types> JCalibration_t;
74 
75  JMultipleFileScanner<> inputFile;
77  JLimit_t& numberOfEvents = inputFile.getLimit();
78  string detectorFile;
79  JCalibration_t calibrationFile;
80  double Tmax_s;
81  JPMTParametersMap pmtParameters;
83  size_t numberOfFits;
84  int debug;
85 
86  try {
87 
88  JParser<> zap("Auxiliary program to convert fit results to Evt format.\
89  \nThe option -L corresponds to the name of a shared library \
90  \nand function so to rearrange the order of fit results.");
91 
92  zap['f'] = make_field(inputFile);
93  zap['o'] = make_field(outputFile) = "aanet.root";
94  zap['n'] = make_field(numberOfEvents) = JLimit::max();
95  zap['a'] = make_field(detectorFile) = "";
96  zap['+'] = make_field(calibrationFile) = JPARSER::initialised();
97  zap['T'] = make_field(Tmax_s) = 100.0;
98  zap['P'] = make_field(pmtParameters, "PMT simulation data (or corresponding file name)") = JPARSER::initialised();
100  zap['N'] = make_field(numberOfFits) = 0;
101  zap['d'] = make_field(debug) = 2;
102 
103  zap(argc, argv);
104  }
105  catch(const exception& error) {
106  FATAL(error.what() << endl);
107  }
108 
109  if (detectorFile == "" && !calibrationFile.empty()) {
110  FATAL("Missing detector file." << endl);
111  }
112 
113  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
114  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
115 
117 
118  if (detectorFile != "") {
119  try {
120  load(detectorFile, detector);
121  }
122  catch(const JException& error) {
123  FATAL(error);
124  }
125  }
126 
127  unique_ptr<JDynamics> dynamics;
128 
129  try {
130 
131  dynamics.reset(new JDynamics(detector, Tmax_s));
132 
133  dynamics->load(calibrationFile);
134  }
135  catch(const exception& error) {
136  if (!calibrationFile.empty()) {
137  FATAL(error.what());
138  }
139  }
140 
141 
142  const JPMTRouter pmt_router(dynamics ? dynamics->getDetector() : detector);
143  const JModuleRouter mod_router(dynamics ? dynamics->getDetector() : detector);
144 
145 
146  outputFile.open();
147 
148  Head header;
149 
150  try {
151  header = getHeader(inputFile);
152  } catch(const exception& error) {}
153 
154  JAANET::JHead buffer(header);
155 
156  // DAQ livetime is not yet set for real data, for Monte Carlo, it is set in JTriggerEfficiency
157 
158  if (buffer.pull(&JAANET::JHead::DAQ) == buffer.end()) {
159 
160  buffer.DAQ.livetime_s = getLivetime(inputFile.begin(), inputFile.end());
161  buffer.push(&JAANET::JHead::DAQ);
162 
163  copy(buffer, header);
164  }
165 
166  if (detectorFile != "") {
167 
168  buffer.calibration.buffer = (dynamics ? calibration::dynamical() : calibration::statical());
169  buffer.push(&JAANET::JHead::calibration);
170 
171  copy(buffer, header);
172  }
173 
174 
175  outputFile.put(header);
176 
177 
178  outputFile.put(JMeta(argc, argv));
179 
180 
181  counter_type number_of_events = 0;
182 
183  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
184 
185  STATUS("Processing: " << *i << endl);
186 
187  JParallelFileScanner_t in(*i);
188  JTreeScanner<Evt> mc(*i);
189 
190  in.setLimit(inputFile.getLimit());
191 
192  Vec center(0,0,0);
193 
194  int mc_run_id = 0;
195 
196  try {
197 
198  const JAANET::JHead head(getHeader(*i));
199 
200  center = get<Vec>(head);
201  mc_run_id = head.start_run.run_id;
202 
203  } catch(const exception& error) {}
204 
205 
206  while (in.hasNext()) {
207 
208  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
209 
210  multi_pointer_type ps = in.next();
211 
212  JDAQEvent* tev = ps;
213  JFIT::JEvt* evt = ps;
214 
215  if (dynamics) {
216  dynamics->update(*tev); // update detector
217  }
218 
219  JFIT::JEvt::iterator __end = evt->end();
220 
221  if (numberOfFits > 0) {
222  advance(__end = evt->begin(), min(numberOfFits, evt->size()));
223  }
224 
225  if (qualitySorter.is_valid()) {
226  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
227  }
228 
229  Evt out; // output event
230 
231  if (mc.getEntries() != 0) {
232 
233  Evt* event = mc.getEntry(tev->getCounter());
234 
235  if (event != NULL) {
236 
237  out = *event; // import Monte Carlo true data
238 
239  if (out.mc_run_id == 0) {
240  out.mc_run_id = mc_run_id;
241  }
242 
243  for (vector<Trk>::iterator i = out.mc_trks.begin(); i != out.mc_trks.end(); ++i) {
244  i->pos += center; // offset true positions
245  }
246  }
247  }
248 
249  read(out, *tev); // import DAQ data
250 
251  if (!pmt_router->empty()) { // import calibration data
252 
253  for (vector<Hit>::iterator i = out.mc_hits.begin(); i != out.mc_hits.end(); ++i) {
254 
255  if (pmt_router.hasPMT(i->pmt_id)) {
256 
257  const JPMTAddress address = pmt_router.getAddress(i->pmt_id);
258  const JPMTIdentifier id = pmt_router.getIdentifier(address);
259  const JPMT& pmt = pmt_router.getPMT(address);
260 
261  i->dom_id = id.getID();
262  i->channel_id = id.getPMTAddress();
263  i->pos = getPosition (pmt);
264  i->dir = getDirection(pmt);
265 
266  } else {
267 
268  FATAL("Missing PMT" << i->pmt_id << endl);
269  }
270  }
271  }
272 
273  if (!mod_router->empty()) { // import calibration data
274 
275  for (vector<Hit>::iterator i = out.hits.begin(); i != out.hits.end(); ++i) {
276 
277  if (mod_router.hasModule(i->dom_id)) {
278 
279  const JPMTIdentifier id(i->dom_id, i->channel_id);
280 
281  const JPMTParameters& parameters = pmtParameters.getPMTParameters(id);
282  const JPMTAnalogueSignalProcessor cpu(parameters);
283 
284  const JPMT& pmt = mod_router.getPMT(id);
285 
287 
288  const JTRIGGER::JHit hit(getTime(i->tdc, pmt.getCalibration()), i->tot);
289 
290  i->pmt_id = pmt.getID();
291  i->pos = getPosition (pmt);
292  i->dir = getDirection(pmt);
293  i->t = hit.getT();
294  i->tot = hit.getToT();
295  i->a = cpu.getNPE(i->tot);
296 
297  } else {
298 
299  FATAL("Missing module " << i->dom_id << endl);
300  }
301  }
302  }
303 
304  struct : public set<JUUID> {
305 
306  inline int get_index(const JUUID& element) const
307  {
308  const_iterator i = this->find(element);
309 
310  if (i != this->end())
311  return std::distance(this->begin(), i);
312  else
313  return -1;
314  }
315  } uuid;
316 
317  for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
318  uuid.insert(fit->getUUID());
319  }
320 
321  for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) { // import reconstruction data
322 
323  Trk trk;
324 
325  trk.id = uuid.get_index(fit->getUUID());
326  trk.pos = Vec(fit->getX(), fit->getY(), fit->getZ());
327  trk.dir = Vec(fit->getDX(), fit->getDY(), fit->getDZ());
328  trk.t = fit->getT();
329  trk.E = fit->getE();
330  trk.lik = fit->getQ();
332  trk.status = (fit->getStatus() == COMPLETE_CHAIN ? TRK_ST_FINALSTATE : TRK_ST_UNDEFINED);
333 
334  if (fit->hasParentUUID()) {
335  trk.mother_id = uuid.get_index(fit->getParentUUID());
336  }
337 
338  for (JHistory::const_iterator i = fit->getHistory().begin(); i != fit->getHistory().end(); ++i) {
339  trk.rec_stages.push_back(i->type);
340  }
341 
342  for (int i = 0; i != fit->getN(); ++i) {
343  trk.fitinf.push_back(fit->getW(i));
344  }
345 
346  trk.error_matrix = fit->getV();
347 
348  out.trks.push_back(trk);
349  }
350 
351  out.id = ++number_of_events;
352 
353  outputFile.put(out);
354  }
355  }
356 
357  STATUS(endl);
358 
360 
361  io >> outputFile;
362 
363  outputFile.close();
364 }
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
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:24
std::istream & read(std::istream &in, JTestSummary &summary, const char delimiter= ' ')
Read test summary.
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
Detector data structure.
Definition: JDetector.hh:89
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.
*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 getTime(const Hit &hit)
Get true time of hit.
string outputFile
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
Calibration.
Definition: JHead.hh:328
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
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
bool slewing
time slewing of analogue signal
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
int getID() const
Get identifier.
Definition: JObjectID.hh:50
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:43
JPosition3D getPosition(const Vec &pos)
Get position.
std::vector< double > fitinf
place to store additional fit info, see km3net-dataformat/definitions/fitparameters.csv
Definition: Trk.hh:32
Auxiliary class for map of PMT parameters.
static const int TRK_ST_UNDEFINED
status was not defined for this MC track (all reco tracks have this value)
Definition: trkmembers.hh:14
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:65
#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
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
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.2.0-6-ge33aafe 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
int id
offline event identifier
Definition: Evt.hh:22
Orientation of module.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
Data structure for PMT parameters.
std::vector< Hit > hits
list of hits
Definition: Evt.hh:38
do set_variable DETECTOR_TXT $WORKDIR detector
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