Jpp
 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 
10 #include "JLang/JComparator.hh"
11 
12 #include "JDB/JToAshort.hh"
13 #include "JDB/JSupport.hh"
14 
15 #include "JDetector/JDetector.hh"
18 #include "JDetector/JTripod.hh"
19 #include "JDetector/JHydrophone.hh"
20 
21 #include "JLang/JPredicate.hh"
22 #include "JLang/JComparator.hh"
23 #include "JTools/JRange.hh"
24 #include "JTools/JQuantile.hh"
25 #include "JTools/JHashMap.hh"
28 #include "JSupport/JMeta.hh"
29 
30 #include "JAcoustics/JEmitterID.hh"
36 #include "JAcoustics/JEvent.hh"
39 #include "JAcoustics/JSupport.hh"
40 
41 #include "Jeep/JContainer.hh"
42 #include "Jeep/JPrint.hh"
43 #include "Jeep/JParser.hh"
44 #include "Jeep/JMessage.hh"
45 
46 
47 /**
48  * \file
49  *
50  * Main program to trigger acoustic data.
51  *
52  * An acoustic event is based on a coincidence of estimated times-of-emission from the same emitter.\n
53  * If the number of coincident times-of-emission exceeds a preset minimum,
54  * the event is triggered and subsequently stored in the output file.\n
55  * Note that an event counter is added which can be used to disambiguate
56  * the different cycles from the same emitter within one "super" event.\n
57  * In this, a "super" event refers to the coincidence of events from a variety of emitters.
58  * \author mdejong
59  */
60 int main(int argc, char **argv)
61 {
62  using namespace std;
63  using namespace JPP;
64 
66 
67  typedef JContainer< vector<JTripod> > tripods_container;
68  typedef JContainer< vector<JHydrophone> > hydrophones_container;
69 
71  JLimit_t& numberOfEvents = inputFile.getLimit();
74  JSoundVelocity V = getSoundVelocity; // default sound velocity
75  string detectorFile;
76  tripods_container tripods; // tripods
77  hydrophones_container hydrophones; // hydrophones
78  double precision;
79  int debug;
80 
81  try {
82 
83  JParser<> zap("Main program to trigger acoustic data.");
84 
85  zap['f'] = make_field(inputFile, "output of JConvertDB -q toashort");
86  zap['n'] = make_field(numberOfEvents) = JLimit::max();
87  zap['@'] = make_field(parameters, "trigger parameters");
88  zap['o'] = make_field(outputFile, "output file") = "event.root";
89  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
90  zap['a'] = make_field(detectorFile, "detector file");
91  zap['T'] = make_field(tripods, "tripod data");
92  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
93  zap['p'] = make_field(precision, "precision time-of-arrival") = 1.0e-6;
94  zap['d'] = make_field(debug) = 1;
95 
96  zap(argc, argv);
97  }
98  catch(const exception &error) {
99  FATAL(error.what() << endl);
100  }
101 
102 
104 
105  try {
106  load(detectorFile, detector);
107  }
108  catch(const JException& error) {
109  FATAL(error);
110  }
111 
113 
114  V.set(detector.getUTMZ());
115 
116  JHashMap<int, JReceiver> receivers;
117  JHashMap<int, JEmitter> emitters;
118 
119  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
120 
121  JPosition3D pos(0.0, 0.0, 0.0);
122 
123  if (i->getFloor() == 0) { // get relative position of hydrophone
124 
125  try {
126  pos = getPosition(hydrophones.begin(),
127  hydrophones.end(),
128  make_predicate(&JHydrophone::getString, i->getString()));
129  }
130  catch(const exception&) {
131  continue; // if no position avialable, discard hydrophone
132  }
133  }
134 
135  receivers[i->getID()] = JReceiver(i->getID(),
136  i->getPosition() + pos,
137  i->getT0() * 1.0e-9);
138  }
139 
140  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
141  emitters[i->getID()] = JEmitter(i->getID(),
142  i->getUTMPosition() - detector.getUTMPosition());
143  }
144 
145 
146  outputFile.open();
147 
148  if (!outputFile.is_open()) {
149  FATAL("Error opening file " << outputFile << endl);
150  }
151 
152  outputFile.put(JMeta(argc, argv));
153  outputFile.put(parameters);
154 
155  // statistics
156 
157  struct map_type :
158  public map<int, JQuantile>
159  {
160  long long int sum() const
161  {
162  long long int count = 0;
163 
164  for (const_iterator i = this->begin(); i != this->end(); ++i) {
165  count += i->second.getCount();
166  }
167 
168  return count;
169  }
170  };
171 
172  map_type number_of_entries;
173  map_type number_of_triggers;
174  map_type number_of_events;
175 
176  // input data
177 
178  string oid; // detector identifier
179 
180  typedef vector<JTransmission> buffer_type; // data type
181 
182  map<int, map< int, buffer_type > > f1; // emitter -> receiver -> data
183 
184  while (inputFile.hasNext()) {
185 
186  if (inputFile.getCounter()%1000 == 0) {
187  STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
188  }
189 
190  JToAshort* parameters = inputFile.next();
191 
192  if (oid == "") {
193  oid = parameters->DETID;
194  }
195 
196  if (oid != parameters->DETID) { // consistency check
197  FATAL("Invalid detector identifier " << parameters->DETID << " != " << oid << endl);
198  }
199 
200  const int id = getEmitterID(parameters->EMITTERID);
201 
202  number_of_entries[parameters->EMITTERID].put(1);
203 
204  if (emitters.has(id) && receivers.has(parameters->DOMID)) {
205 
206  const JTransceiver transceiver(emitters[id], receivers[parameters->DOMID]);
207 
208  f1[transceiver.emitter.getID()][transceiver.receiver.getID()].push_back(transceiver.getTransmission(*parameters, V));
209  }
210  }
211  STATUS(endl);
212 
213 
214  // filter similar hits
215 
216  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
217 
218  for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
219 
220  buffer_type& buffer = receiver->second;
221 
222  sort(buffer.begin(), buffer.end(), JTransmission::compare(precision));
223 
224  buffer.erase(unique(buffer.begin(), buffer.end(), JTransmission::equals(precision)), buffer.end());
225  }
226  }
227 
228 
229  // output events
230 
231  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
232 
233  buffer_type buffer;
234 
235  for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
236 
237  // selection based on quality
238 
239  buffer_type& f2 = receiver->second;
240  double Qmin = parameters.Q;
241 
242  if (parameters.Q > 0.0 && parameters.Q < 1.0) {
243 
244  const JQuantile Q1("quality", make_array(f2.begin(), f2.end(), &JTransmission::getQ), true);
245 
246  Qmin = Q1.getQuantile(parameters.Q);
247  }
248 
249  buffer_type::iterator __end = partition(f2.begin(), f2.end(), make_predicate(&JTransmission::getQ, Qmin, JComparison::gt()));
250 
251  copy(f2.begin(), __end, back_inserter(buffer));
252  }
253 
254 
255  sort(buffer.begin(), buffer.end()); // sort according time-of-emisson
256 
257 
258  JTriggerOutput data;
259 
260  for (buffer_type::const_iterator p = buffer.begin(); p != buffer.end(); ++p) {
261 
262  buffer_type::const_iterator q = p;
263 
264  // basic correlator
265 
266  while (++q != buffer.end() && q->getToE() - p->getToE() <= parameters.TMax_s) {}
267 
268  if (distance(p,q) >= parameters.numberOfHits) {
269 
270  STATUS("trigger: " << setw(8) << number_of_triggers.sum() << '\r' << flush); DEBUG(endl);
271 
272  JEvent event(oid, number_of_triggers.sum(), i->first, p, q);
273 
274  number_of_triggers[i->first].put(event.size());
275 
276  DEBUG("trigger: " << endl << event);
277 
278  data.push_back(event);
279  }
280  }
281 
282  data.merge(JEventOverlap(parameters.TMax_s));
283 
284  for (JTriggerOutput::const_iterator event = data.begin(); event != data.end(); ++event) {
285 
286  number_of_events[event->getID()].put(event->size());
287 
288  outputFile.put(*event);
289  }
290  }
291  STATUS(endl);
292 
293  for (map_type::const_iterator i = number_of_entries.begin(); i != number_of_entries.end(); ++i) {
294  STATUS("number of entries: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
295  }
296 
297  for (map_type::const_iterator i = number_of_triggers.begin(); i != number_of_triggers.end(); ++i) {
298  STATUS("number of triggers: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
299  }
300 
301  for (map_type::const_iterator i = number_of_events.begin(); i != number_of_events.end(); ++i) {
302  STATUS("number of events: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << ' ' << FIXED(7,3) << i->second.getMean() << endl);
303  }
304 
305  JMultipleFileScanner<JMeta> io(inputFile);
306 
307  io >> outputFile;
308 
309  outputFile.close();
310 }
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
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:71
static JEmitterID & getEmitterID
Function object for emitter identification.
Definition: JEmitterID.hh:90
Quantile calculator.
Definition: JQuantile.hh:88
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:80
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.
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:445
string outputFile
Data structure for detector geometry and calibration.
double getQuantile(const double Q, const bool reverse=false) const
Get quantile.
Definition: JQuantile.hh:353
Data structure for hydrophone.
Type list.
Definition: JTypeList.hh:22
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:196
Acoustics toolkit.
Detector toolkit.
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:1961
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition: JVectorize.hh:37
int getID() const
Get identifier.
Definition: JObjectID.hh:50
ROOT I/O of application specific meta data.
JPosition3D getPosition(const Vec &pos)
Get position.
Acoustic transceiver.
Definition: JTransceiver.hh:29
ROOT TTree parameter settings.
int debug
debug level
Definition: JSirene.cc:63
General purpose messaging.
Implementation for 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 configure(const T &value, const JAbstractCollection< JAbscissa_t > &bounds, JBool< false > option)
Configuration of value.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Acoustic trigger parameters.
std::vector< int > count
Definition: JAlgorithm.hh:184
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.
JTransmission getTransmission(const JToAshort &data, const JAbstractSoundVelocity &V=getSoundVelocity) const
Get transmission.
Definition: JTransceiver.hh:62
Acoustic transceiver.
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.
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
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