Jpp  18.3.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"
19 #include "JDetector/JTripod.hh"
21 #include "JDetector/JHydrophone.hh"
22 
23 #include "JLang/JPredicate.hh"
24 #include "JLang/JComparator.hh"
25 #include "JTools/JRange.hh"
26 #include "JTools/JQuantile.hh"
27 #include "JTools/JHashMap.hh"
30 #include "JSupport/JMeta.hh"
31 
32 #include "JAcoustics/JEmitterID.hh"
38 #include "JAcoustics/JEvent.hh"
41 #include "JAcoustics/JSupport.hh"
42 
43 #include "Jeep/JContainer.hh"
44 #include "Jeep/JPrint.hh"
45 #include "Jeep/JParser.hh"
46 #include "Jeep/JMessage.hh"
47 
48 
49 /**
50  * \file
51  *
52  * Main program to trigger events in acoustic data.
53  *
54  * An acoustic event is based on a coincidence of estimated times of emission from the same emitter.\n
55  * If the number of coincident times of emission exceeds a preset minimum,
56  * the event is triggered and subsequently stored in the output file.\n
57  * Note that an event counter is added which can be used to disambiguate
58  * the different cycles from the same emitter within one "super" event that is used for the global fit.\n
59  * In this, a super event refers to the coincidence of events from a variety of emitters.
60  * \author mdejong
61  */
62 int main(int argc, char **argv)
63 {
64  using namespace std;
65  using namespace JPP;
66 
68 
72 
74  JLimit_t& numberOfEvents = inputFile.getLimit();
77  JSoundVelocity V = getSoundVelocity; // default sound velocity
78  string detectorFile;
79  tripods_container tripods; // tripods
80  transmitters_container transmitters; // transmitters
81  hydrophones_container hydrophones; // hydrophones
82  double precision;
83  int debug;
84 
85  try {
86 
87  JParser<> zap("Main program to trigger events in acoustic data.");
88 
89  zap['f'] = make_field(inputFile, "output of JConvertDB -q toashort");
90  zap['n'] = make_field(numberOfEvents) = JLimit::max();
91  zap['@'] = make_field(parameters, "trigger parameters");
92  zap['o'] = make_field(outputFile, "output file") = "event.root";
93  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
94  zap['a'] = make_field(detectorFile, "detector file");
95  zap['T'] = make_field(tripods, "tripod data");
96  zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
97  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
98  zap['W'] = make_field(getEmitterID, "waveform identification data") = JPARSER::initialised();
99  zap['p'] = make_field(precision, "precision time-of-arrival") = 1.0e-6;
100  zap['d'] = make_field(debug) = 1;
101 
102  zap(argc, argv);
103  }
104  catch(const exception &error) {
105  FATAL(error.what() << endl);
106  }
107 
108 
110 
111  try {
112  load(detectorFile, detector);
113  }
114  catch(const JException& error) {
115  FATAL(error);
116  }
117 
118  V.set(detector.getUTMZ());
119 
120  JHashMap<int, JReceiver> receivers;
121  JHashMap<int, JEmitter> emitters;
122 
123  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
124 
125  if ((i->getFloor() == 0 && !i->has(HYDROPHONE_DISABLE)) ||
126  (i->getFloor() != 0 && !i->has(PIEZO_DISABLE))) {
127 
128  JPosition3D pos(0.0, 0.0, 0.0);
129 
130  if (i->getFloor() == 0) { // get relative position of hydrophone
131 
132  try {
133  pos = getPosition(hydrophones.begin(),
134  hydrophones.end(),
135  make_predicate(&JHydrophone::getString, i->getString()));
136  }
137  catch(const exception&) {
138  continue; // if no position available, discard hydrophone
139  }
140  }
141 
142  receivers[i->getID()] = JReceiver(i->getID(),
143  i->getPosition() + pos,
144  i->getT0() * 1.0e-9);
145  }
146  }
147 
148  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
149  emitters[i->getID()] = JEmitter(i->getID(),
150  i->getUTMPosition() - detector.getUTMPosition());
151  }
152 
153  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
154  try {
155  emitters[i->getID()] = JEmitter(i->getID(),
156  i->getPosition() + detector.getModule(i->getLocation()).getPosition());
157  }
158  catch(const exception&) {
159  continue; // if no string available, discard transmitter
160  }
161  }
162 
163  const JRange<double> unit(0.0, 1.0);
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()%10000 == 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  const JEventOverlap match(parameters.TMax_s);
248 
249  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
250 
251  buffer_type buffer;
252 
253  for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
254 
255  // filter similar hits
256 
257  sort(receiver->second.begin(), receiver->second.end(), JTransmission::compare(precision));
258 
259  buffer_type::iterator __end = unique(receiver->second.begin(), receiver->second.end(), JTransmission::equals(precision));
260 
261  // selection based on quality
262 
263  for (buffer_type::const_iterator p = receiver->second.begin(); p != __end; ++p) {
264  if (p->getQ() >= parameters.Q * (unit(parameters.Q) ? p->getW() : 1.0)) {
265  buffer.push_back(*p);
266  }
267  }
268  }
269 
270  sort(buffer.begin(), buffer.end()); // sort according time-of-emisson
271 
272  JQuantile Q; // monitor of trigger
273  JEvent out; // temporary storage of event(s)
274  JTriggerOutput data; //
275 
276  for (buffer_type::const_iterator p = buffer.begin(); p != buffer.end(); ++p) {
277 
278  buffer_type::const_iterator q = p;
279 
280  // basic correlator
281 
282  while (++q != buffer.end() && q->getToE() - p->getToE() <= parameters.TMax_s) {}
283 
284  Q.put(distance(p,q));
285 
286  if (distance(p,q) >= parameters.numberOfHits) {
287 
288  JEvent event(oid, number_of_triggers.sum(), i->first, p, q);
289 
290  number_of_triggers[i->first].put(event.size());
291 
292  DEBUG("trigger: " << endl << event);
293 
294  if (out.empty()) {
295 
296  out = event;
297 
298  } else if (match(out, event)) {
299 
300  out.merge(event);
301 
302  } else {
303 
304  data.push_back(out);
305 
306  out = event;
307  }
308  }
309  }
310 
311  if (!out.empty()) {
312  data.push_back(out); // pending event
313  }
314 
315  STATUS("number of toes: " << setw(3) << i->first << ' ' << setw(8) << buffer.size());
316  if (Q.getCount() > 0) {
317  STATUS(' ' << FIXED(6,1) << Q.getMean() <<
318  ' ' << FIXED(6,1) << Q.getXmin() <<
319  ' ' << FIXED(6,1) << Q.getXmax());
320  }
321  STATUS(endl);
322 
323  for (JTriggerOutput::const_iterator event = data.begin(); event != data.end(); ++event) {
324 
325  number_of_events[event->getID()].put(event->size());
326 
327  outputFile.put(*event);
328  }
329  }
330  STATUS(endl);
331 
332  for (map_type::const_iterator i = number_of_triggers.begin(); i != number_of_triggers.end(); ++i) {
333  STATUS("number of triggers: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
334  }
335 
336  for (map_type::const_iterator i = number_of_events.begin(); i != number_of_events.end(); ++i) {
337  STATUS("number of events: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << ' ' << FIXED(7,3) << i->second.getMean() << endl);
338  }
339 
340  JMultipleFileScanner<JMeta> io(inputFile);
341 
342  io >> outputFile;
343 
344  outputFile.close();
345 
346  return 0;
347 }
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
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
JContainer< std::vector< JTransmitter > > transmitters_container
Definition: JSydney.cc:79
void merge(const JEvent &event)
Merge event.
std::map< int, buffer_type > map_type
string -&gt; hits
Definition: JPerth.cc:68
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
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:83
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.
JMODEL::JString getString(const JFit &fit)
Get model parameters of string.
std::vector< JHitW0 > buffer_type
hits
Definition: JPerth.cc:67
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.
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition: JSydney.cc:78
Data structure for transmitter.
Acoustic emitter.
Definition: JEmitter.hh:27
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:1989
int getID() const
Get identifier.
Definition: JObjectID.hh:50
ROOT I/O of application specific meta data.
JPosition3D getPosition(const Vec &pos)
Get position.
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.
JContainer< std::vector< JTripod > > tripods_container
Definition: JSydney.cc:77
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...
Auxiliary data structure for average.
Definition: JKatoomba_t.hh:76
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.
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.
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:84
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.
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62