Jpp  17.1.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JAcousticsEventBuilder.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 
11 
12 #include "JLang/JComparator.hh"
13 
14 #include "JDB/JToAshort.hh"
15 #include "JDB/JSupport.hh"
16 
17 #include "JDetector/JDetector.hh"
20 #include "JDetector/JTripod.hh"
22 #include "JDetector/JHydrophone.hh"
23 
24 #include "JLang/JPredicate.hh"
25 #include "JLang/JComparator.hh"
26 #include "JTools/JRange.hh"
27 #include "JTools/JQuantile.hh"
28 #include "JTools/JHashMap.hh"
31 #include "JSupport/JMeta.hh"
32 
33 #include "JAcoustics/JEmitterID.hh"
39 #include "JAcoustics/JEvent.hh"
42 #include "JAcoustics/JSupport.hh"
43 
44 #include "Jeep/JContainer.hh"
45 #include "Jeep/JPrint.hh"
46 #include "Jeep/JParser.hh"
47 #include "Jeep/JMessage.hh"
48 
49 
50 /**
51  * \file
52  *
53  * Main program to trigger events in acoustic data.
54  *
55  * An acoustic event is based on a coincidence of estimated times of emission from the same emitter.\n
56  * If the number of coincident times of emission exceeds a preset minimum,
57  * the event is triggered and subsequently stored in the output file.\n
58  * Note that an event counter is added which can be used to disambiguate
59  * the different cycles from the same emitter within one so-called "super" event that is used for the global fit.\n
60  * In this, a super event refers to the coincidence of events from a variety of emitters.
61  * \author mdejong
62  */
63 int main(int argc, char **argv)
64 {
65  using namespace std;
66  using namespace JPP;
67 
69 
70  typedef JContainer< vector<JTripod> > tripods_container;
71  typedef JContainer< vector<JTransmitter> > transmitters_container;
72  typedef JContainer< vector<JHydrophone> > hydrophones_container;
73 
75  JLimit_t& numberOfEvents = inputFile.getLimit();
78  JSoundVelocity V = getSoundVelocity; // default sound velocity
79  string detectorFile;
80  tripods_container tripods; // tripods
81  transmitters_container transmitters; // transmitters
82  hydrophones_container hydrophones; // hydrophones
83  double precision;
84  int debug;
85 
86  try {
87 
88  JParser<> zap("Main program to trigger events in acoustic data.");
89 
90  zap['f'] = make_field(inputFile, "output of JConvertDB -q toashort");
91  zap['n'] = make_field(numberOfEvents) = JLimit::max();
92  zap['@'] = make_field(parameters, "trigger parameters");
93  zap['o'] = make_field(outputFile, "output file") = "event.root";
94  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
95  zap['a'] = make_field(detectorFile, "detector file");
96  zap['T'] = make_field(tripods, "tripod data");
97  zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
98  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
99  zap['W'] = make_field(getEmitterID, "waveform identification data") = JPARSER::initialised();
100  zap['p'] = make_field(precision, "precision time-of-arrival") = 1.0e-6;
101  zap['d'] = make_field(debug) = 1;
102 
103  zap(argc, argv);
104  }
105  catch(const exception &error) {
106  FATAL(error.what() << endl);
107  }
108 
109 
111 
112  try {
113  load(detectorFile, detector);
114  }
115  catch(const JException& error) {
116  FATAL(error);
117  }
118 
119  V.set(detector.getUTMZ());
120 
121  JHashMap<int, JReceiver> receivers;
122  JHashMap<int, JEmitter> emitters;
123 
124  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
125 
126  if ((i->getFloor() == 0 && !i->has(HYDROPHONE_DISABLE)) ||
127  (i->getFloor() != 0 && !i->has(PIEZO_DISABLE))) {
128 
129  JPosition3D pos(0.0, 0.0, 0.0);
130 
131  if (i->getFloor() == 0) { // get relative position of hydrophone
132 
133  try {
134  pos = getPosition(hydrophones.begin(),
135  hydrophones.end(),
136  make_predicate(&JHydrophone::getString, i->getString()));
137  }
138  catch(const exception&) {
139  continue; // if no position available, discard hydrophone
140  }
141  }
142 
143  receivers[i->getID()] = JReceiver(i->getID(),
144  i->getPosition() + pos,
145  i->getT0() * 1.0e-9);
146  }
147  }
148 
149  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
150  emitters[i->getID()] = JEmitter(i->getID(),
151  i->getUTMPosition() - detector.getUTMPosition());
152  }
153 
154  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
155  try {
156  emitters[i->getID()] = JEmitter(i->getID(),
157  i->getPosition() + detector.getModule(i->getLocation()).getPosition());
158  }
159  catch(const exception&) {
160  continue; // if no string available, discard transmitter
161  }
162  }
163 
164 
165  outputFile.open();
166 
167  if (!outputFile.is_open()) {
168  FATAL("Error opening file " << outputFile << endl);
169  }
170 
171  outputFile.put(JMeta(argc, argv));
172  outputFile.put(parameters);
173 
174  // statistics
175 
176  struct map_type :
177  public map<int, JQuantile>
178  {
179  long long int sum() const
180  {
181  long long int count = 0;
182 
183  for (const_iterator i = this->begin(); i != this->end(); ++i) {
184  count += i->second.getCount();
185  }
186 
187  return count;
188  }
189  };
190 
191  map_type number_of_entries;
192  map_type number_of_aliens;
193  map_type number_of_triggers;
194  map_type number_of_events;
195 
196  // input data
197 
198  string oid; // detector identifier
199 
200  typedef vector<JTransmission> buffer_type; // data type
201 
202  map<int, map< int, buffer_type > > f1; // emitter -> receiver -> data
203 
204  while (inputFile.hasNext()) {
205 
206  if (inputFile.getCounter()%1000 == 0) {
207  STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
208  }
209 
210  JToAshort* parameters = inputFile.next();
211 
212  if (oid == "") {
213  oid = parameters->DETID;
214  }
215 
216  if (oid != parameters->DETID) { // consistency check
217  FATAL("Invalid detector identifier " << parameters->DETID << " != " << oid << endl);
218  }
219 
220  try {
221 
222  const int id = getEmitterID(parameters->EMITTERID);
223 
224  number_of_entries[parameters->EMITTERID].put(1);
225 
226  if (emitters.has(id) && receivers.has(parameters->DOMID)) {
227 
228  const JTransceiver transceiver(emitters[id], receivers[parameters->DOMID]);
229 
230  f1[transceiver.emitter.getID()][transceiver.receiver.getID()].push_back(transceiver.getTransmission(*parameters, V));
231  }
232  }
233  catch(const exception&) {
234  number_of_aliens[parameters->EMITTERID].put(1);
235  }
236  }
237  STATUS(endl);
238 
239  for (map_type::const_iterator i = number_of_entries.begin(); i != number_of_entries.end(); ++i) {
240  STATUS("number of entries: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
241  }
242 
243  for (map_type::const_iterator i = number_of_aliens.begin(); i != number_of_aliens.end(); ++i) {
244  STATUS("number of aliens: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
245  }
246 
247 
248  // filter similar hits
249 
250  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
251 
252  for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
253 
254  buffer_type& buffer = receiver->second;
255 
256  sort(buffer.begin(), buffer.end(), JTransmission::compare(precision));
257 
258  buffer.erase(unique(buffer.begin(), buffer.end(), JTransmission::equals(precision)), buffer.end());
259  }
260  }
261 
262 
263  // output events
264 
265  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
266 
267  buffer_type buffer;
268 
269  for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
270 
271  // selection based on quality
272 
273  buffer_type& f2 = receiver->second;
274  double Qmin = parameters.Q;
275 
276  if (parameters.Q > 0.0 && // use quantile to evaluate corresponding quality
277  parameters.Q < 1.0) { //
278 
279  const JQuantile Q1("quality", make_array(f2.begin(), f2.end(), &JTransmission::getQ), true);
280 
281  Qmin = Q1.getQuantile(parameters.Q);
282  }
283 
284  buffer_type::iterator __end = partition(f2.begin(), f2.end(), make_predicate(&JTransmission::getQ, Qmin, JComparison::gt()));
285 
286  copy(f2.begin(), __end, back_inserter(buffer));
287  }
288 
289 
290  sort(buffer.begin(), buffer.end()); // sort according time-of-emisson
291 
292  JQuantile Q;
293 
294  JTriggerOutput data;
295 
296  for (buffer_type::const_iterator p = buffer.begin(); p != buffer.end(); ++p) {
297 
298  buffer_type::const_iterator q = p;
299 
300  // basic correlator
301 
302  while (++q != buffer.end() && q->getToE() - p->getToE() <= parameters.TMax_s) {}
303 
304  Q.put(distance(p,q));
305 
306  if (distance(p,q) >= parameters.numberOfHits) {
307 
308  STATUS("trigger: " << setw(8) << number_of_triggers.sum() << '\r' << flush); DEBUG(endl);
309 
310  JEvent event(oid, number_of_triggers.sum(), i->first, p, q);
311 
312  number_of_triggers[i->first].put(event.size());
313 
314  DEBUG("trigger: " << endl << event);
315 
316  data.push_back(event);
317  }
318  }
319 
320  STATUS("number of toes: " << setw(3) << i->first << ' ' << setw(8) << buffer.size());
321  if (Q.getCount() > 0) {
322  STATUS(' ' << FIXED(6,1) << Q.getMean() <<
323  ' ' << FIXED(6,1) << Q.getXmin() <<
324  ' ' << FIXED(6,1) << Q.getXmax());
325  }
326  STATUS(endl);
327 
328  data.merge(JEventOverlap(parameters.TMax_s));
329 
330  for (JTriggerOutput::const_iterator event = data.begin(); event != data.end(); ++event) {
331 
332  number_of_events[event->getID()].put(event->size());
333 
334  outputFile.put(*event);
335  }
336  }
337  STATUS(endl);
338 
339  for (map_type::const_iterator i = number_of_triggers.begin(); i != number_of_triggers.end(); ++i) {
340  STATUS("number of triggers: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
341  }
342 
343  for (map_type::const_iterator i = number_of_events.begin(); i != number_of_events.end(); ++i) {
344  STATUS("number of events: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << ' ' << FIXED(7,3) << i->second.getMean() << endl);
345  }
346 
347  JMultipleFileScanner<JMeta> io(inputFile);
348 
349  io >> outputFile;
350 
351  outputFile.close();
352 }
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
Acoustic receiver.
Definition: JReceiver.hh:27
Q(UTCMax_s-UTCMin_s)-livetime_s
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
Greater than.
Definition: JComparison.hh:73
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:72
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
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
Recording of objects on file according a format that follows from the file name extension.
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.
then fatal Number of tripods
Definition: JFootprint.sh:45
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
string outputFile
Data structure for detector geometry and calibration.
Data structure for hydrophone.
Type list.
Definition: JTypeList.hh:22
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.
I/O formatting auxiliaries.
Detector file.
Definition: JHead.hh:226
Acoustics support kit.
Detector support kit.
Data structure for transmitter.
Acoustic emitter.
Definition: JEmitter.hh:27
void merge(const JMatch_t &match)
Merge events.
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:1961
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.
JPosition3D getPosition(const Vec &pos)
Get position.
double getQuantile(const double Q, const Quantile_t option=forward_t) const
Get quantile.
Definition: JQuantile.hh:319
JTransmission getTransmission(const JToAshort &data, const JAbstractSoundVelocity &V) const
Get transmission.
Definition: JTransceiver.hh:62
Acoustic transceiver.
Definition: JTransceiver.hh:29
Emitter identification.
ROOT TTree parameter settings.
int debug
debug level
Definition: JSirene.cc:67
double getQ(const double D_m, const double f_kHz, const double d_m)
Get relative quality for given frequency at given distance.
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...
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Acoustic trigger parameters.
std::vector< int > count
Definition: JAlgorithm.hh:180
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.
Acoustic transceiver.
static const int PIEZO_DISABLE
Enable (disable) use of piezo if this status bit is 0 (1);.
Acoustic transmission.
Auxiliary class to compare transmissions.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:139
Acoustic event.
static const int HYDROPHONE_DISABLE
Enable (disable) use of hydrophone if this status bit is 0 (1);.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
bool equals(const JFirst_t &first, const JSecond_t &second, const double precision=std::numeric_limits< double >::min())
Check equality.
Definition: JMathToolkit.hh:86
static JEmitterID getEmitterID
Function object for emitter identification.
Definition: JEmitterID.hh:119
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for tripod.
Match of two events considering overlap in time.
int EMITTERID
waveform identifier
Definition: JToAshort.hh:27
Container I/O.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62