Jpp  15.0.5
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"
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 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 so-called "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 
69  typedef JContainer< vector<JTripod> > tripods_container;
70  typedef JContainer< vector<JHydrophone> > hydrophones_container;
71 
73  JLimit_t& numberOfEvents = inputFile.getLimit();
76  JSoundVelocity V = getSoundVelocity; // default sound velocity
77  string detectorFile;
78  tripods_container tripods; // tripods
79  hydrophones_container hydrophones; // hydrophones
80  double precision;
81  int debug;
82 
83  try {
84 
85  JParser<> zap("Main program to trigger acoustic data.");
86 
87  zap['f'] = make_field(inputFile, "output of JConvertDB -q toashort");
88  zap['n'] = make_field(numberOfEvents) = JLimit::max();
89  zap['@'] = make_field(parameters, "trigger parameters");
90  zap['o'] = make_field(outputFile, "output file") = "event.root";
91  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
92  zap['a'] = make_field(detectorFile, "detector file");
93  zap['T'] = make_field(tripods, "tripod data");
94  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
95  zap['W'] = make_field(getEmitterID, "waveform identification data") = JPARSER::initialised();
96  zap['p'] = make_field(precision, "precision time-of-arrival") = 1.0e-6;
97  zap['d'] = make_field(debug) = 1;
98 
99  zap(argc, argv);
100  }
101  catch(const exception &error) {
102  FATAL(error.what() << endl);
103  }
104 
105 
107 
108  try {
109  load(detectorFile, detector);
110  }
111  catch(const JException& error) {
112  FATAL(error);
113  }
114 
115  V.set(detector.getUTMZ());
116 
117  JHashMap<int, JReceiver> receivers;
118  JHashMap<int, JEmitter> emitters;
119 
120  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
121 
122  if ((i->getFloor() == 0 && !i->has(HYDROPHONE_DISABLE)) ||
123  (i->getFloor() != 0 && !i->has(PIEZO_DISABLE))) {
124 
125  JPosition3D pos(0.0, 0.0, 0.0);
126 
127  if (i->getFloor() == 0) { // get relative position of hydrophone
128 
129  try {
130  pos = getPosition(hydrophones.begin(),
131  hydrophones.end(),
132  make_predicate(&JHydrophone::getString, i->getString()));
133  }
134  catch(const exception&) {
135  continue; // if no position avialable, discard hydrophone
136  }
137  }
138 
139  receivers[i->getID()] = JReceiver(i->getID(),
140  i->getPosition() + pos,
141  i->getT0() * 1.0e-9);
142  }
143  }
144 
145  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
146  emitters[i->getID()] = JEmitter(i->getID(),
147  i->getUTMPosition() - detector.getUTMPosition());
148  }
149 
150 
151  outputFile.open();
152 
153  if (!outputFile.is_open()) {
154  FATAL("Error opening file " << outputFile << endl);
155  }
156 
157  outputFile.put(JMeta(argc, argv));
158  outputFile.put(parameters);
159 
160  // statistics
161 
162  struct map_type :
163  public map<int, JQuantile>
164  {
165  long long int sum() const
166  {
167  long long int count = 0;
168 
169  for (const_iterator i = this->begin(); i != this->end(); ++i) {
170  count += i->second.getCount();
171  }
172 
173  return count;
174  }
175  };
176 
177  map_type number_of_entries;
178  map_type number_of_triggers;
179  map_type number_of_events;
180 
181  // input data
182 
183  string oid; // detector identifier
184 
185  typedef vector<JTransmission> buffer_type; // data type
186 
187  map<int, map< int, buffer_type > > f1; // emitter -> receiver -> data
188 
189  while (inputFile.hasNext()) {
190 
191  if (inputFile.getCounter()%1000 == 0) {
192  STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
193  }
194 
195  JToAshort* parameters = inputFile.next();
196 
197  if (oid == "") {
198  oid = parameters->DETID;
199  }
200 
201  if (oid != parameters->DETID) { // consistency check
202  FATAL("Invalid detector identifier " << parameters->DETID << " != " << oid << endl);
203  }
204 
205  const int id = getEmitterID(parameters->EMITTERID);
206 
207  number_of_entries[parameters->EMITTERID].put(1);
208 
209  if (emitters.has(id) && receivers.has(parameters->DOMID)) {
210 
211  const JTransceiver transceiver(emitters[id], receivers[parameters->DOMID]);
212 
213  f1[transceiver.emitter.getID()][transceiver.receiver.getID()].push_back(transceiver.getTransmission(*parameters, V));
214  }
215  }
216  STATUS(endl);
217 
218  for (map_type::const_iterator i = number_of_entries.begin(); i != number_of_entries.end(); ++i) {
219  STATUS("number of entries: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
220  }
221 
222 
223  // filter similar hits
224 
225  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
226 
227  for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
228 
229  buffer_type& buffer = receiver->second;
230 
231  sort(buffer.begin(), buffer.end(), JTransmission::compare(precision));
232 
233  buffer.erase(unique(buffer.begin(), buffer.end(), JTransmission::equals(precision)), buffer.end());
234  }
235  }
236 
237 
238  // output events
239 
240  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
241 
242  buffer_type buffer;
243 
244  for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
245 
246  // selection based on quality
247 
248  buffer_type& f2 = receiver->second;
249  double Qmin = parameters.Q;
250 
251  if (parameters.Q > 0.0 && // use quantile to evaluate corresponding quality
252  parameters.Q < 1.0) { //
253 
254  const JQuantile Q1("quality", make_array(f2.begin(), f2.end(), &JTransmission::getQ), true);
255 
256  Qmin = Q1.getQuantile(parameters.Q);
257  }
258 
259  buffer_type::iterator __end = partition(f2.begin(), f2.end(), make_predicate(&JTransmission::getQ, Qmin, JComparison::gt()));
260 
261  copy(f2.begin(), __end, back_inserter(buffer));
262  }
263 
264 
265  sort(buffer.begin(), buffer.end()); // sort according time-of-emisson
266 
267  STATUS("number of toes: " << setw(3) << i->first << ' ' << setw(8) << buffer.size() << endl);
268 
269  JTriggerOutput data;
270 
271  for (buffer_type::const_iterator p = buffer.begin(); p != buffer.end(); ++p) {
272 
273  buffer_type::const_iterator q = p;
274 
275  // basic correlator
276 
277  while (++q != buffer.end() && q->getToE() - p->getToE() <= parameters.TMax_s) {}
278 
279  if (distance(p,q) >= parameters.numberOfHits) {
280 
281  STATUS("trigger: " << setw(8) << number_of_triggers.sum() << '\r' << flush); DEBUG(endl);
282 
283  JEvent event(oid, number_of_triggers.sum(), i->first, p, q);
284 
285  number_of_triggers[i->first].put(event.size());
286 
287  DEBUG("trigger: " << endl << event);
288 
289  data.push_back(event);
290  }
291  }
292 
293  data.merge(JEventOverlap(parameters.TMax_s));
294 
295  for (JTriggerOutput::const_iterator event = data.begin(); event != data.end(); ++event) {
296 
297  number_of_events[event->getID()].put(event->size());
298 
299  outputFile.put(*event);
300  }
301  }
302  STATUS(endl);
303 
304  for (map_type::const_iterator i = number_of_triggers.begin(); i != number_of_triggers.end(); ++i) {
305  STATUS("number of triggers: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
306  }
307 
308  for (map_type::const_iterator i = number_of_events.begin(); i != number_of_events.end(); ++i) {
309  STATUS("number of events: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << ' ' << FIXED(7,3) << i->second.getMean() << endl);
310  }
311 
312  JMultipleFileScanner<JMeta> io(inputFile);
313 
314  io >> outputFile;
315 
316  outputFile.close();
317 }
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: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.
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
string outputFile
Data structure for detector geometry and calibration.
double getQuantile(const double Q, const bool reverse=false) const
Get quantile.
Definition: JQuantile.hh:309
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 support kit.
Detector support kit.
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: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.
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:63
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 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);.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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.