Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JAcousticsTriggerProcessor.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 #include <vector>
4 #include <map>
5 #include <algorithm>
6 
7 #include "TROOT.h"
8 #include "TFile.h"
9 
10 #include "JLang/JComparator.hh"
11 
12 #include "JDB/JToAshort.hh"
13 #include "JDB/JSupport.hh"
14 
15 #include "JDetector/JDetector.hh"
17 #include "JDetector/JHydrophone.hh"
18 
19 #include "JLang/JPredicate.hh"
20 #include "JLang/JComparator.hh"
21 #include "JTools/JRange.hh"
22 #include "JTools/JQuantile.hh"
23 #include "JTools/JHashMap.hh"
26 #include "JSupport/JMeta.hh"
27 
28 #include "JAcoustics/JEmitterID.hh"
30 #include "JAcoustics/JReceiver.hh"
35 #include "JAcoustics/JEvent.hh"
37 #include "JAcoustics/JSupport.hh"
38 
39 #include "Jeep/JContainer.hh"
40 #include "Jeep/JProperties.hh"
41 #include "Jeep/JPrint.hh"
42 #include "Jeep/JParser.hh"
43 #include "Jeep/JMessage.hh"
44 
45 #include "JTrigger/JMatch.hh"
46 #include "JTrigger/JAlgorithm.hh"
47 
48 
49 namespace JACOUSTICS {
50 
52  using JTRIGGER::JMatch;
53  using JLANG::JClonable;
54  using JTOOLS::JRange;
55 
56 
57  /**
58  * Acoustic hit.
59  */
60  struct hit_type :
61  public JPosition3D,
62  public JTransmission
63  {
64  /**
65  * Default constructor.
66  */
68  {}
69 
70 
71  /**
72  * Constructor.
73  *
74  * \param position position
75  * \param transmission transmission
76  */
77  hit_type(const JPosition3D& position, const JTransmission& transmission) :
78  JPosition3D (position),
79  JTransmission(transmission)
80  {}
81  };
82 
83 
84  /**
85  * 3D match criterion for acoustic signals.
86  */
87  class JMatch3D :
88  public JClonable< JMatch<hit_type>, JMatch3D>
89  {
90  public:
91  /**
92  * Constructor.
93  *
94  * \param V sound velocity
95  * \param Tmax_s maximal extra time [s]
96  */
97  JMatch3D(const JSoundVelocity& V, const double Tmax_s = 0.0) :
98  V(V),
99  Tmax_s(Tmax_s)
100  {}
101 
102 
103  /**
104  * Match operator.
105  *
106  * \param first hit
107  * \param second hit
108  * \return match result
109  */
110  virtual bool operator()(const hit_type& first, const hit_type& second) const override
111  {
112  using namespace JPP;
113 
114  const double t1 = V.getTime(getDistance(first.getPosition(), second.getPosition()),
115  first .getZ(),
116  second.getZ());
117 
118  const double dt = fabs(first.getToA() - second.getToA());
119 
120  return dt <= t1 + Tmax_s;
121  }
122 
123 
125  const double Tmax_s;
126  };
127 
128 
129  /**
130  * Match of two events considering overlap in time.
131  */
132  struct JEventOverlap {
133 
134  typedef JRange<double> time_range; //!< Type definition of time range.
135 
136  /**
137  * Constructor.
138  *
139  * \param Tmax_s maximal time difference between two consecutive events [s]
140  */
141  JEventOverlap(const double Tmax_s = 0.0) :
142  range(-Tmax_s, +Tmax_s)
143  {}
144 
145 
146  /**
147  * Match criterion.
148  *
149  * \param first first event
150  * \param second second event
151  * \return true if two events overlap in time; else false
152  */
153  bool operator()(const JEvent& first, const JEvent& second) const
154  {
155  using namespace JPP;
156 
157  if (first .empty()) return false;
158  if (second.empty()) return false;
159 
160  time_range r1(make_array(first .begin(), first .end(), &JTransmission::getToA));
161  time_range r2(make_array(second.begin(), second.end(), &JTransmission::getToA));
162 
163  r1.add(range);
164  r2.add(range);
165 
166  return overlap(r1, r2);
167  }
168 
170  };
171 }
172 
173 
174 /**
175  * \file
176  *
177  * Main program to trigger acoustic data.
178  *
179  * An acoustic event is based on coincidences between times of arrival.\n
180  * If the number of coincident times of arrival exceeds a preset minimum,
181  * the event is triggered and subsequently stored in the output file.
182  * \author mdejong
183  */
184 int main(int argc, char **argv)
185 {
186  using namespace std;
187  using namespace JPP;
188 
190 
192 
194  JLimit_t& numberOfEvents = inputFile.getLimit();
196  int factoryLimit = 10000;
197  double TMaxExtra_s = 100.0e-6;
199  JSoundVelocity V = getSoundVelocity; // default sound velocity
200  string detectorFile;
201  hydrophones_container hydrophones; // hydrophones
202  double precision;
203  int debug;
204 
205  try {
206 
207  JProperties properties;
208 
209  properties.insert(gmake_property(parameters.Q));
210  properties.insert(gmake_property(parameters.TMax_s));
211  properties.insert(gmake_property(parameters.numberOfHits));
212  properties.insert(gmake_property(factoryLimit));
213  properties.insert(gmake_property(TMaxExtra_s));
214 
215  JParser<> zap("Main program to trigger acoustic data.");
216 
217  zap['f'] = make_field(inputFile, "output of JConvertDB -q toashort");
218  zap['n'] = make_field(numberOfEvents) = JLimit::max();
219  zap['@'] = make_field(properties, "trigger parameters");
220  zap['o'] = make_field(outputFile, "output file") = "event.root";
221  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
222  zap['a'] = make_field(detectorFile, "detector file");
223  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
224  zap['p'] = make_field(precision, "precision time-of-arrival") = 1.0e-6;
225  zap['W'] = make_field(getEmitterID, "waveform identification data") = JPARSER::initialised();
226  zap['d'] = make_field(debug) = 1;
227 
228  zap(argc, argv);
229  }
230  catch(const exception &error) {
231  FATAL(error.what() << endl);
232  }
233 
234 
236 
237  try {
238  load(detectorFile, detector);
239  }
240  catch(const JException& error) {
241  FATAL(error);
242  }
243 
244  V.set(detector.getUTMZ());
245 
246 
247 
248  const JMatch3D match(V, TMaxExtra_s);
249  const JEventOverlap overlap(parameters.TMax_s);
250 
251 
252  JHashMap<int, JReceiver> receivers;
253 
254  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
255 
256  JPosition3D pos(0.0, 0.0, 0.0);
257 
258  if (i->getFloor() == 0) { // get relative position of hydrophone
259 
260  try {
261  pos = getPosition(hydrophones.begin(),
262  hydrophones.end(),
263  make_predicate(&JHydrophone::getString, i->getString()));
264  }
265  catch(const exception&) {
266  continue; // if no position available, discard hydrophone
267  }
268  }
269 
270  receivers[i->getID()] = JReceiver(i->getID(),
271  i->getPosition() + pos,
272  i->getT0() * 1.0e-9);
273  }
274 
275  const JRange<double> unit(0.0, 1.0);
276 
277  outputFile.open();
278 
279  if (!outputFile.is_open()) {
280  FATAL("Error opening file " << outputFile << endl);
281  }
282 
283  outputFile.put(JMeta(argc, argv));
284  outputFile.put(parameters);
285 
286 
287  // input data
288 
289  string oid; // detector identifier
290 
291  typedef vector<hit_type> buffer_type; // data type
292 
293  map<int, map< int, buffer_type > > f1; // emitter -> receiver -> data
294 
295  while (inputFile.hasNext()) {
296 
297  if (inputFile.getCounter()%100000 == 0) {
298  STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
299  }
300 
301  JToAshort* parameters = inputFile.next();
302 
303  if (oid == "") {
304  oid = parameters->DETID;
305  }
306 
307  if (oid != parameters->DETID) { // consistency check
308  FATAL("Invalid detector identifier " << parameters->DETID << " != " << oid << endl);
309  }
310 
311  if (receivers.has(parameters->DOMID)) {
312 
313  const JReceiver& receiver = receivers[parameters->DOMID];
314 
315  const JTransmission transmission(parameters->RUNNUMBER,
316  parameters->DOMID,
317  parameters->QUALITYFACTOR,
318  parameters->QUALITYNORMALISATION,
319  parameters->UNIXTIMEBASE + receiver.getT(parameters->TOA_S),
320  parameters->UNIXTIMEBASE + receiver.getT(parameters->TOA_S));
321  try {
322  f1[getEmitterID(parameters->EMITTERID)][receiver.getID()].push_back(hit_type(receiver.getPosition(), transmission));
323  }
324  catch(const exception&) {}
325  }
326  }
327  STATUS(endl);
328 
329  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
330 
332 
333  for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
334 
335  // filter similar hits
336 
337  sort(receiver->second.begin(), receiver->second.end(), JTransmission::compare(precision));
338 
339  buffer_type::iterator __end = unique(receiver->second.begin(), receiver->second.end(), JTransmission::equals(precision));
340 
341  // selection based on quality
342 
343  for (buffer_type::const_iterator p = receiver->second.begin(); p != __end; ++p) {
344  if (p->getQ() >= parameters.Q * (unit(parameters.Q) ? p->getW() : 1.0)) {
345  data.push_back(*p);
346  }
347  }
348  }
349 
350  sort(data.begin(), data.end(), make_comparator(&hit_type::getToA, JComparison::lt())); // sort according time-of-arrival
351 
352  buffer_type buffer(factoryLimit); // local buffer of hits
353  JEvent out[2]; // FIFO of events
354  int number_of_events = 0;
355 
356  for (buffer_type::const_iterator p = data.begin(); p != data.end(); ++p) {
357 
358  if (distance(data.cbegin(),p)%1000 == 0) {
359  STATUS("processed[" << i->first << "]: " << FIXED(5,1) << (double) (100 * distance(data.cbegin(),p)) / (double) data.size() << "%" << '\r' << flush); DEBUG(endl);
360  }
361 
362  buffer_type::const_iterator q = p;
363 
364  while (++q != data.end() && q->getToA() - p->getToA() <= parameters.TMax_s) {}
365 
366  if (distance(p,q) >= parameters.numberOfHits) {
367 
368  if (distance(p,q) < factoryLimit) {
369 
370  buffer_type::iterator root = buffer.begin();
371  buffer_type::iterator __p = buffer.begin();
372  buffer_type::iterator __q = buffer.begin();
373 
374  *root = *p;
375 
376  ++__p;
377  ++__q;
378 
379  for (buffer_type::const_iterator i = p; ++i != q; ) {
380  if (match(*p,*i)) {
381  *__q = *i;
382  ++__q;
383  }
384  }
385 
386  if (distance(root,__q) >= parameters.numberOfHits) {
387 
388  __q = clusterize(__p, __q, match, parameters.numberOfHits - 1);
389 
390  if (distance(root,__q) >= parameters.numberOfHits) {
391  out[1] = JEvent(oid, number_of_events, i->first, root, __q);
392  }
393  }
394 
395  } else {
396 
397  out[1] = JEvent(oid, number_of_events, i->first, p, q);
398  }
399 
400  if (out[0].empty()) {
401 
402  out[0] = out[1]; // shift
403 
404  } else if (overlap(out[0],out[1])) {
405 
406  out[0].merge(out[1]); // merge
407 
408  } else {
409 
410  number_of_events += 1;
411 
412  outputFile.put(out[0]); // write
413 
414  out[0] = out[1]; // shift
415  }
416 
417  out[1].clear();
418  }
419  }
420 
421  if (!out[0].empty()) {
422 
423  number_of_events += 1;
424 
425  outputFile.put(out[0]); // write
426  }
427  STATUS(endl);
428  STATUS("triggers[" << i->first << "]: " << setw(7) << number_of_events << endl);
429  }
430 
431  JMultipleFileScanner<JMeta> io(inputFile);
432 
433  io >> outputFile;
434 
435  outputFile.close();
436 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
hit_type(const JPosition3D &position, const JTransmission &transmission)
Constructor.
Object writing to file.
double UNIXTIMEBASE
[s]
Definition: JToAshort.hh:25
Utility class to parse command line options.
Definition: JParser.hh:1711
General exception.
Definition: JException.hh:24
bool operator()(const JEvent &first, const JEvent &second) const
Match criterion.
Acoustic receiver.
Definition: JReceiver.hh:27
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
int main(int argc, char *argv[])
Definition: Main.cc:15
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Sound velocity.
void merge(const JEvent &event)
Merge event.
Algorithms for hit clustering and sorting.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
JVector3D getPosition(T __begin, T __end, const JPredicate< JTypename_t, JComparator_t > &predicate)
Get position from element in data which corresponds to given predicate.
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:72
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
General purpose class for hash map of unique elements.
#define STATUS(A)
Definition: JMessage.hh:63
ROOT TTree parameter settings.
Detector data structure.
Definition: JDetector.hh:89
Match of two events considering overlap in time.
Recording of objects on file according a format that follows from the file name extension.
virtual bool operator()(const hit_type &first, const hit_type &second) const override
Match operator.
Utility class to parse parameter values.
Definition: JProperties.hh:497
std::string DETID
constraint
Definition: JToAshort.hh:22
*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
Acoustic event.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:84
Function object interface for hit matching.
Definition: JMatch.hh:25
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
3D match criterion for acoustic signals.
JMatch3D(const JSoundVelocity &V, const double Tmax_s=0.0)
Constructor.
string outputFile
Data structure for detector geometry and calibration.
do JPlot2D f $WORKDIR canberra[${EMITTER}\] root
std::vector< JHitW0 > buffer_type
hits
Definition: JPerth.cc:70
Utility class to parse parameter values.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Data structure for hydrophone.
Type list.
Definition: JTypeList.hh:22
JFIT::JEvent JEvent
Definition: JHistory.hh:353
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
static const JSoundVelocity getSoundVelocity(1541.0,-17.0e-3,-2000.0)
Function object for velocity of sound.
Base class for match operations for cluster and hit-preprocessing methods.
Acoustic transmission.
I/O formatting auxiliaries.
Detector file.
Definition: JHead.hh:226
Acoustics support kit.
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition: JSydney.cc:79
Acoustics toolkit.
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition: JContainer.hh:39
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
range_type & add(argument_type x)
Add offset.
Definition: JRange.hh:446
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition: JVectorize.hh:54
int getID() const
Get identifier.
Definition: JObjectID.hh:50
ROOT I/O of application specific meta data.
Template class for object cloning.
Definition: JClonable.hh:20
Emitter identification.
ROOT TTree parameter settings.
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
Range of values.
Definition: JRange.hh:38
General purpose messaging.
Implementation for depth dependend velocity of sound.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
double getT(const double t_s) const
Get corrected time.
Definition: JReceiver.hh:72
Acoustic receiver.
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
int getString() const
Get string number.
Definition: JLocation.hh:134
JRange< double > time_range
Type definition of time range.
virtual double getTime(const double D_m, const double z1, const double z2) const override
Get propagation time of sound.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Acoustic trigger parameters.
Auxiliary class to define a range between two values.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
hit_type()
Default constructor.
double getToA() const
Get calibrated time of arrival.
Acoustic transmission.
Auxiliary class to compare transmissions.
Acoustic event.
double TOA_S
[s]
Definition: JToAshort.hh:28
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
bool equals(const JFirst_t &first, const JSecond_t &second, const double precision=std::numeric_limits< double >::min())
Check equality.
Definition: JMathToolkit.hh:86
bool overlap(const JRange< T, JComparator_t > &first, const JRange< T, JComparator_t > &second)
Test overlap between ranges.
Definition: JRange.hh:641
JEventOverlap(const double Tmax_s=0.0)
Constructor.
static JEmitterID getEmitterID
Function object for emitter identification.
Definition: JEmitterID.hh:119
do set_variable DETECTOR_TXT $WORKDIR detector
JSoundVelocity & set(const double z0)
Set depth.
double getZ() const
Get z position.
Definition: JVector3D.hh:115
int EMITTERID
waveform identifier
Definition: JToAshort.hh:27
static struct JTRIGGER::clusterize clusterize
Container I/O.
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62