Jpp  18.1.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 
110  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
111  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
112 
113 
115 
116  if (detectorFile != "") {
117  try {
118  load(detectorFile, detector);
119  }
120  catch(const JException& error) {
121  FATAL(error);
122  }
123  }
124 
125  unique_ptr<JDynamics> dynamics;
126 
127  try {
128 
129  dynamics.reset(new JDynamics(detector, Tmax_s));
130 
131  dynamics->load(calibrationFile);
132  }
133  catch(const exception& error) {
134  if (!calibrationFile.empty()) {
135  FATAL(error.what());
136  }
137  }
138 
139 
140  const JPMTRouter pmt_router(dynamics ? dynamics->getDetector() : detector);
141  const JModuleRouter mod_router(dynamics ? dynamics->getDetector() : detector);
142 
143 
144  outputFile.open();
145 
146  Head header;
147 
148  try {
149  header = getHeader(inputFile);
150  } catch(const exception& error) {}
151 
152  JAANET::JHead buffer(header);
153 
154  if (buffer.pull(&JAANET::JHead::DAQ) == buffer.end()) {
155 
156  buffer.DAQ.livetime_s = getLivetime(inputFile.begin(), inputFile.end());
157  buffer.push(&JAANET::JHead::DAQ);
158 
159  copy(buffer, header);
160  }
161 
162  outputFile.put(header);
163 
164 
165  outputFile.put(JMeta(argc, argv));
166 
167 
168  counter_type number_of_events = 0;
169 
170  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
171 
172  STATUS("Processing: " << *i << endl);
173 
174  JParallelFileScanner_t in(*i);
175  JTreeScanner<Evt> mc(*i);
176 
177  in.setLimit(inputFile.getLimit());
178 
179  Vec center(0,0,0);
180 
181  int mc_run_id = 0;
182 
183  try {
184 
185  const JAANET::JHead head(getHeader(*i));
186 
187  center = get<Vec>(head);
188  mc_run_id = head.start_run.run_id;
189 
190  } catch(const exception& error) {}
191 
192 
193  while (in.hasNext()) {
194 
195  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
196 
197  multi_pointer_type ps = in.next();
198 
199  JDAQEvent* tev = ps;
200  JFIT::JEvt* evt = ps;
201 
202  if (dynamics) {
203  dynamics->update(*tev); // update detector
204  }
205 
206  JFIT::JEvt::iterator __end = evt->end();
207 
208  if (numberOfFits > 0) {
209  advance(__end = evt->begin(), min(numberOfFits, evt->size()));
210  }
211 
212  if (qualitySorter.is_valid()) {
213  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
214  }
215 
216  Evt out; // output event
217 
218  if (mc.getEntries() != 0) {
219 
220  Evt* event = mc.getEntry(tev->getCounter());
221 
222  if (event != NULL) {
223 
224  out = *event; // import Monte Carlo true data
225 
226  if (out.mc_run_id == 0) {
227  out.mc_run_id = mc_run_id;
228  }
229 
230  for (vector<Trk>::iterator i = out.mc_trks.begin(); i != out.mc_trks.end(); ++i) {
231  i->pos += center; // offset true positions
232  }
233  }
234  }
235 
236  read(out, *tev); // import DAQ data
237 
238  if (!pmt_router->empty()) { // import calibration data
239 
240  for (vector<Hit>::iterator i = out.mc_hits.begin(); i != out.mc_hits.end(); ++i) {
241 
242  if (pmt_router.hasPMT(i->pmt_id)) {
243 
244  const JPMTAddress address = pmt_router.getAddress(i->pmt_id);
245  const JPMTIdentifier id = pmt_router.getIdentifier(address);
246  const JPMT& pmt = pmt_router.getPMT(address);
247 
248  i->dom_id = id.getID();
249  i->channel_id = id.getPMTAddress();
250  i->pos = getPosition (pmt);
251  i->dir = getDirection(pmt);
252 
253  } else {
254 
255  FATAL("Missing PMT" << i->pmt_id << endl);
256  }
257  }
258  }
259 
260  if (!mod_router->empty()) { // import calibration data
261 
262  for (vector<Hit>::iterator i = out.hits.begin(); i != out.hits.end(); ++i) {
263 
264  if (mod_router.hasModule(i->dom_id)) {
265 
266  const JPMTIdentifier id(i->dom_id, i->channel_id);
267 
268  const JPMTParameters& parameters = pmtParameters.getPMTParameters(id);
269  const JPMTAnalogueSignalProcessor cpu(parameters);
270 
271  const JPMT& pmt = mod_router.getPMT(id);
272 
274 
275  const JTRIGGER::JHit hit(getTime(i->tdc, pmt.getCalibration()), i->tot);
276 
277  i->pmt_id = pmt.getID();
278  i->pos = getPosition (pmt);
279  i->dir = getDirection(pmt);
280  i->t = hit.getT();
281  i->tot = hit.getToT();
282  i->a = cpu.getNPE(i->tot);
283 
284  } else {
285 
286  FATAL("Missing module " << i->dom_id << endl);
287  }
288  }
289  }
290 
291  struct : public set<JUUID> {
292 
293  inline int get_index(const JUUID& element) const
294  {
295  const_iterator i = this->find(element);
296 
297  if (i != this->end())
298  return std::distance(this->begin(), i);
299  else
300  return -1;
301  }
302  } uuid;
303 
304  for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
305  uuid.insert(fit->getUUID());
306  }
307 
308  for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) { // import reconstruction data
309 
310  Trk trk;
311 
312  trk.id = uuid.get_index(fit->getUUID());
313  trk.pos = Vec(fit->getX(), fit->getY(), fit->getZ());
314  trk.dir = Vec(fit->getDX(), fit->getDY(), fit->getDZ());
315  trk.t = fit->getT();
316  trk.E = fit->getE();
317  trk.lik = fit->getQ();
319  trk.status = (fit->getStatus() == COMPLETE_CHAIN ? TRK_ST_FINALSTATE : TRK_ST_UNDEFINED);
320 
321  if (fit->hasParentUUID()) {
322  trk.mother_id = uuid.get_index(fit->getParentUUID());
323  }
324 
325  for (JHistory::const_iterator i = fit->getHistory().begin(); i != fit->getHistory().end(); ++i) {
326  trk.rec_stages.push_back(i->type);
327  }
328 
329  for (int i = 0; i != fit->getN(); ++i) {
330  trk.fitinf.push_back(fit->getW(i));
331  }
332 
333  trk.error_matrix = fit->getV();
334 
335  out.trks.push_back(trk);
336  }
337 
338  out.id = ++number_of_events;
339 
340  outputFile.put(out);
341  }
342  }
343 
344  STATUS(endl);
345 
347 
348  io >> outputFile;
349 
350  outputFile.close();
351 }
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.
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
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.
Monte Carlo run header.
Definition: JHead.hh:1221
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
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.1.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
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
JAANET::DAQ DAQ
Definition: JHead.hh:1589
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:46
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
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