Jpp  17.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/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 64 of file software/JReconstruction/JConvertEvt.cc.

65 {
66  using namespace std;
67  using namespace JPP;
68  using namespace KM3NETDAQ;
69 
71  JACOUSTICS::JEvt>::typelist calibration_types;
72  typedef JMultipleFileScanner<calibration_types> JCalibration_t;
73 
74  JMultipleFileScanner<> inputFile;
76  JLimit_t& numberOfEvents = inputFile.getLimit();
77  string detectorFile;
78  JCalibration_t calibrationFile;
79  double Tmax_s;
80  JPMTParametersMap pmtParameters;
82  size_t numberOfFits;
83  int debug;
84 
85  try {
86 
87  JParser<> zap("Auxiliary program to convert fit results to Evt format.\
88  \nThe option -L corresponds to the name of a shared library \
89  \nand function so to rearrange the order of fit results.");
90 
91  zap['f'] = make_field(inputFile);
92  zap['o'] = make_field(outputFile) = "aanet.root";
93  zap['n'] = make_field(numberOfEvents) = JLimit::max();
94  zap['a'] = make_field(detectorFile) = "";
95  zap['+'] = make_field(calibrationFile) = JPARSER::initialised();
96  zap['T'] = make_field(Tmax_s) = 100.0;
97  zap['P'] = make_field(pmtParameters, "PMT simulation data (or corresponding file name)") = JPARSER::initialised();
99  zap['N'] = make_field(numberOfFits) = 0;
100  zap['d'] = make_field(debug) = 2;
101 
102  zap(argc, argv);
103  }
104  catch(const exception& error) {
105  FATAL(error.what() << endl);
106  }
107 
108 
109  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
110  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
111 
112 
114 
115  if (detectorFile != "") {
116  try {
117  load(detectorFile, detector);
118  }
119  catch(const JException& error) {
120  FATAL(error);
121  }
122  }
123 
124  getMechanics.load(detector.getID());
125 
126  unique_ptr<JDynamics> dynamics;
127 
128  try {
129 
130  dynamics.reset(new JDynamics(detector, Tmax_s));
131 
132  dynamics->load(calibrationFile);
133  }
134  catch(const exception& error) {
135  if (!calibrationFile.empty()) {
136  FATAL(error.what());
137  }
138  }
139 
140 
141  const JPMTRouter pmt_router(dynamics ? dynamics->getDetector() : detector);
142  const JModuleRouter mod_router(dynamics ? dynamics->getDetector() : detector);
143 
144 
145  outputFile.open();
146 
147  Head header;
148 
149  try {
150  header = getHeader(inputFile);
151  } catch(const exception& error) {}
152 
153  JAANET::JHead buffer(header);
154 
155  if (buffer.pull(&JAANET::JHead::DAQ) == buffer.end()) {
156 
157  buffer.DAQ.livetime_s = getLivetime(inputFile.begin(), inputFile.end());
158  buffer.push(&JAANET::JHead::DAQ);
159 
160  copy(buffer, header);
161  }
162 
163  outputFile.put(header);
164 
165 
166  outputFile.put(JMeta(argc, argv));
167 
168 
169  counter_type number_of_events = 0;
170 
171  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
172 
173  STATUS("Processing: " << *i << endl);
174 
175  JParallelFileScanner_t in(*i);
176  JTreeScanner<Evt> mc(*i);
177 
178  in.setLimit(inputFile.getLimit());
179 
180  Vec center(0,0,0);
181 
182  int mc_run_id = 0;
183 
184  try {
185 
186  const JAANET::JHead head(getHeader(*i));
187 
188  center = get<Vec>(head);
189  mc_run_id = head.start_run.run_id;
190 
191  } catch(const exception& error) {}
192 
193 
194  while (in.hasNext()) {
195 
196  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
197 
198  multi_pointer_type ps = in.next();
199 
200  JDAQEvent* tev = ps;
201  JFIT::JEvt* evt = ps;
202 
203  if (dynamics) {
204  dynamics->update(*tev); // update detector
205  }
206 
207  JFIT::JEvt::iterator __end = evt->end();
208 
209  if (numberOfFits > 0) {
210  advance(__end = evt->begin(), min(numberOfFits, evt->size()));
211  }
212 
213  if (qualitySorter.is_valid()) {
214  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
215  }
216 
217  Evt out; // output event
218 
219  if (mc.getEntries() != 0) {
220 
221  Evt* event = mc.getEntry(tev->getCounter());
222 
223  if (event != NULL) {
224 
225  out = *event; // import Monte Carlo true data
226 
227  if (out.mc_run_id == 0) {
228  out.mc_run_id = mc_run_id;
229  }
230 
231  for (vector<Trk>::iterator i = out.mc_trks.begin(); i != out.mc_trks.end(); ++i) {
232  i->pos += center; // offset true positions
233  }
234  }
235  }
236 
237  read(out, *tev); // import DAQ data
238 
239  if (!pmt_router->empty()) { // import calibration data
240 
241  for (vector<Hit>::iterator i = out.mc_hits.begin(); i != out.mc_hits.end(); ++i) {
242 
243  if (pmt_router.hasPMT(i->pmt_id)) {
244 
245  const JPMTAddress address = pmt_router.getAddress(i->pmt_id);
246  const JPMTIdentifier id = pmt_router.getIdentifier(address);
247  const JPMT& pmt = pmt_router.getPMT(address);
248 
249  i->dom_id = id.getID();
250  i->channel_id = id.getPMTAddress();
251  i->pos = getPosition (pmt);
252  i->dir = getDirection(pmt);
253 
254  } else {
255 
256  FATAL("Missing PMT" << i->pmt_id << endl);
257  }
258  }
259  }
260 
261  if (!mod_router->empty()) { // import calibration data
262 
263  for (vector<Hit>::iterator i = out.hits.begin(); i != out.hits.end(); ++i) {
264 
265  if (mod_router.hasModule(i->dom_id)) {
266 
267  const JPMTIdentifier id(i->dom_id, i->channel_id);
268 
269  const JPMTParameters& parameters = pmtParameters.getPMTParameters(id);
270  const JPMTAnalogueSignalProcessor cpu(parameters);
271 
272  const JPMT& pmt = mod_router.getPMT(id);
273 
275 
276  const JTRIGGER::JHit hit(getTime(i->tdc, pmt.getCalibration()), i->tot);
277 
278  i->pmt_id = pmt.getID();
279  i->pos = getPosition (pmt);
280  i->dir = getDirection(pmt);
281  i->t = hit.getT();
282  i->tot = hit.getToT();
283  i->a = cpu.getNPE(i->tot);
284 
285  } else {
286 
287  FATAL("Missing module " << i->dom_id << endl);
288  }
289  }
290  }
291 
292  struct : public set<JUUID> {
293 
294  inline int get_index(const JUUID& element) const
295  {
296  const_iterator i = this->find(element);
297 
298  if (i != this->end())
299  return std::distance(this->begin(), i);
300  else
301  return -1;
302  }
303  } uuid;
304 
305  for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
306  uuid.insert(fit->getUUID());
307  }
308 
309  for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) { // import reconstruction data
310 
311  Trk trk;
312 
313  trk.id = uuid.get_index(fit->getUUID());
314  trk.pos = Vec(fit->getX(), fit->getY(), fit->getZ());
315  trk.dir = Vec(fit->getDX(), fit->getDY(), fit->getDZ());
316  trk.t = fit->getT();
317  trk.E = fit->getE();
318  trk.lik = fit->getQ();
320  trk.status = (fit->getStatus() == COMPLETE_CHAIN ? TRK_ST_FINALSTATE : TRK_ST_UNDEFINED);
321 
322  if (fit->hasParentUUID()) {
323  trk.mother_id = uuid.get_index(fit->getParentUUID());
324  }
325 
326  for (JHistory::const_iterator i = fit->getHistory().begin(); i != fit->getHistory().end(); ++i) {
327  trk.rec_stages.push_back(i->type);
328  }
329 
330  for (int i = 0; i != fit->getN(); ++i) {
331  trk.fitinf.push_back(fit->getW(i));
332  }
333 
334  trk.error_matrix = fit->getV();
335 
336  out.trks.push_back(trk);
337  }
338 
339  out.id = ++number_of_events;
340 
341  outputFile.put(out);
342  }
343  }
344 
345  STATUS(endl);
346 
348 
349  io >> outputFile;
350 
351  outputFile.close();
352 }
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:1500
General exception.
Definition: JException.hh:23
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
static JDetectorMechanics getMechanics
Function object to get string mechanics.
Definition: JMechanics.hh:243
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:66
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:1961
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, for jgandalf, see JFitParameters.hh
Definition: Trk.hh:32
Auxiliary class for map of PMT parameters.
void load(const std::string &file_name)
Load mechanical model parameters from file.
Definition: JMechanics.hh:142
int debug
debug level
Definition: JSirene.cc:67
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:1167
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:63
static const int TRK_ST_FINALSTATE
particle to be tracked by detector-level MC (&#39;track_in&#39; tag in evt files from gseagen ...
Definition: trkmembers.hh:15
Vec pos
postion of the track at time t [m]
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 v2.2.0-12-ge740342 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:139
int id
offline event identifier
Definition: Evt.hh:22
Orientation of module.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
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:1518
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
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
identifyer for the overall fitting algorithm/chain/strategy
Definition: Trk.hh:25