Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
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
17#include "JDetector/JTripod.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/JToA.hh"
37#include "JAcoustics/JEvent.hh"
41
42#include "Jeep/JContainer.hh"
43#include "Jeep/JPrint.hh"
44#include "Jeep/JParser.hh"
45#include "Jeep/JMessage.hh"
46
47
48/**
49 * \file
50 *
51 * Main program to trigger events in acoustic data.
52 *
53 * An acoustic event is based on a coincidence of estimated times of emission from the same emitter.\n
54 * If the number of coincident times of emission exceeds a preset minimum,
55 * the event is triggered and subsequently stored in the output file.\n
56 * Note that an event counter is added which can be used to disambiguate
57 * the different cycles from the same emitter within one "super" event that is used for the global fit.\n
58 * In this, a super event refers to the coincidence of events from a variety of emitters.
59 * \author mdejong
60 */
61int main(int argc, char **argv)
62{
63 using namespace std;
64 using namespace JPP;
65
67
71
73 JLimit_t& numberOfEvents = inputFile.getLimit();
74 JTriggerParameters parameters;
76 JSoundVelocity V = getSoundVelocity; // default sound velocity
77 string detectorFile;
78 tripods_container tripods; // tripods
79 transmitters_container transmitters; // transmitters
80 hydrophones_container hydrophones; // hydrophones
81 double precision;
82 int debug;
83
84 try {
85
86 JParser<> zap("Main program to trigger events in acoustic data.");
87
88 zap['f'] = make_field(inputFile, "output of JToA");
89 zap['n'] = make_field(numberOfEvents) = JLimit::max();
90 zap['@'] = make_field(parameters, "trigger parameters");
91 zap['o'] = make_field(outputFile, "output file") = "event.root";
92 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
93 zap['a'] = make_field(detectorFile, "detector file");
94 zap['T'] = make_field(tripods, "tripod data");
95 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
96 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
97 zap['W'] = make_field(getEmitterID, "waveform identification data") = JPARSER::initialised();
98 zap['p'] = make_field(precision, "precision time-of-arrival") = 1.0e-6;
99 zap['d'] = make_field(debug) = 1;
100
101 zap(argc, argv);
102 }
103 catch(const exception &error) {
104 FATAL(error.what() << endl);
105 }
106
107
109
110 try {
111 load(detectorFile, detector);
112 }
113 catch(const JException& error) {
114 FATAL(error);
115 }
116
117 V.set(detector.getUTMZ());
118
119 JHashMap<int, JReceiver> receivers;
121
122 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
123
124 if ((i->getFloor() == 0 && !i->has(HYDROPHONE_DISABLE)) ||
125 (i->getFloor() != 0 && !i->has(PIEZO_DISABLE))) {
126
128
129 if (i->getFloor() == 0) { // get relative position of hydrophone
130
131 try {
132 pos = getPosition(hydrophones.begin(),
133 hydrophones.end(),
134 make_predicate(&JHydrophone::getString, i->getString()));
135 }
136 catch(const exception&) {
137 continue; // if no position available, discard hydrophone
138 }
139 }
140
141 receivers[i->getID()] = JReceiver(i->getID(),
142 i->getPosition() + pos,
143 i->getT0() * 1.0e-9);
144 }
145 }
146
147 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
148 emitters[i->getID()] = JEmitter(i->getID(),
149 i->getUTMPosition() - detector.getUTMPosition());
150 }
151
152 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
153 try {
154 emitters[i->getID()] = JEmitter(i->getID(),
155 i->getPosition() + detector.getModule(i->getLocation()).getPosition());
156 }
157 catch(const exception&) {
158 continue; // if no string available, discard transmitter
159 }
160 }
161
162 outputFile.open();
163
164 if (!outputFile.is_open()) {
165 FATAL("Error opening file " << outputFile << endl);
166 }
167
168 outputFile.put(JMeta(argc, argv));
169 outputFile.put(parameters);
170
171 // statistics
172
173 struct map_type :
174 public map<int, JQuantile>
175 {
176 long long int sum() const
177 {
178 long long int count = 0;
179
180 for (const_iterator i = this->begin(); i != this->end(); ++i) {
181 count += i->second.getCount();
182 }
183
184 return count;
185 }
186 };
187
188 map_type number_of_entries;
189 map_type number_of_aliens;
190 map_type number_of_misses;
191 map_type number_of_triggers;
192 map_type number_of_events;
193
194 // input data
195
196 typedef vector<JTransmission> buffer_type; // data type
197
198 map<int, map< int, buffer_type > > f1; // emitter -> receiver -> data
199
200 while (inputFile.hasNext()) {
201
202 if (inputFile.getCounter()%10000 == 0) {
203 STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
204 }
205
206 JToA* parameters = inputFile.next();
207
208 if (detector.getID() != parameters->DETID) { // consistency check
209 FATAL("Invalid detector identifier " << parameters->DETID << " != " << detector.getID() << endl);
210 }
211
212 try {
213
214 const int id = getEmitterID(parameters->WAVEFORMID);
215
216 number_of_entries[parameters->WAVEFORMID].put(1);
217
218 if (emitters.has(id)) {
219
220 if (receivers.has(parameters->DOMID)) {
221
222 const JTransceiver transceiver(emitters[id], receivers[parameters->DOMID]);
223
224 f1[transceiver.emitter.getID()][transceiver.receiver.getID()].push_back(transceiver.getTransmission(*parameters, V));
225 }
226
227 } else {
228 number_of_misses[id].put(1);
229 }
230 }
231 catch(const exception&) {
232 number_of_aliens[parameters->WAVEFORMID].put(1);
233 }
234 }
235 STATUS(endl);
236
237 for (map_type::const_iterator i = number_of_entries.begin(); i != number_of_entries.end(); ++i) {
238 STATUS("number of entries: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
239 }
240
241 for (map_type::const_iterator i = number_of_aliens.begin(); i != number_of_aliens.end(); ++i) {
242 STATUS("number of aliens: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
243 }
244
245 for (map_type::const_iterator i = number_of_misses.begin(); i != number_of_misses.end(); ++i) {
246 STATUS("number of misses: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
247 }
248
249 const JEventOverlap match(parameters.TMax_s);
250
251 for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
252
253 buffer_type buffer;
254
255 for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
256
257 // filter similar hits
258
259 sort(receiver->second.begin(), receiver->second.end(), JTransmission::compare(precision));
260
261 buffer_type::iterator __end = unique(receiver->second.begin(), receiver->second.end(), JTransmission::equals(precision));
262
263 // selection based on quality
264
265 for (buffer_type::const_iterator p = receiver->second.begin(); p != __end; ++p) {
266 if (p->getQ() >= parameters.Q * (parameters.Q <= 1.0 ? p->getW() : 1.0)) {
267 buffer.push_back(*p);
268 }
269 }
270 }
271
272 sort(buffer.begin(), buffer.end()); // sort according time-of-emisson
273
274 JQuantile Q; // monitor of trigger
275 JEvent out; // temporary storage of event(s)
276
277 for (buffer_type::const_iterator p = buffer.begin(); p != buffer.end(); ++p) {
278
279 buffer_type::const_iterator q = p;
280
281 // basic correlator
282
283 while (++q != buffer.end() && q->getToE() - p->getToE() <= parameters.TMax_s) {}
284
285 Q.put(distance(p,q));
286
287 if (distance(p,q) >= parameters.numberOfHits) {
288
289 JEvent event(detector.getID(), number_of_triggers.sum(), i->first, p, q);
290
291 number_of_triggers[i->first].put(event.size());
292
293 DEBUG("trigger: " << endl << event);
294
295 if (out.empty()) {
296
297 out = event;
298
299 } else if (match(out, event)) {
300
301 out.merge(event);
302
303 } else {
304
305 outputFile.put(out);
306
307 number_of_events[out.getID()].put(out.size());
308
309 out = event;
310 }
311 }
312 }
313
314 if (!out.empty()) {
315
316 outputFile.put(out); // pending event
317
318 number_of_events[out.getID()].put(out.size());
319 }
320
321 STATUS("number of toes: " << setw(3) << i->first << ' ' << setw(8) << buffer.size());
322 if (Q.getCount() > 0) {
323 STATUS(' ' << FIXED(6,1) << Q.getMean() <<
324 ' ' << FIXED(6,1) << Q.getXmin() <<
325 ' ' << FIXED(6,1) << Q.getXmax());
326 }
327 STATUS(endl);
328 }
329 STATUS(endl);
330
331 for (map_type::const_iterator i = number_of_triggers.begin(); i != number_of_triggers.end(); ++i) {
332 STATUS("number of triggers: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
333 }
334
335 for (map_type::const_iterator i = number_of_events.begin(); i != number_of_events.end(); ++i) {
336 STATUS("number of events: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << ' ' << FIXED(7,3) << i->second.getMean() << endl);
337 }
338
339 JMultipleFileScanner<JMeta> io(inputFile);
340
341 io >> outputFile;
342
343 outputFile.close();
344
345 return 0;
346}
int main(int argc, char **argv)
Acoustics support kit.
Acoustics toolkit.
Acoustic event.
ROOT TTree parameter settings.
Acoustic trigger parameters.
Container I/O.
string outputFile
Data structure for detector geometry and calibration.
Emitter identification.
Recording of objects on file according a format that follows from the file name extension.
General purpose class for hash map of unique elements.
Data structure for hydrophone.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
ROOT I/O of application specific meta data.
Module support kit.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
I/O formatting auxiliaries.
Auxiliary class to define a range between two values.
Acoustic event.
Acoustic transceiver.
Acoustic transmission.
Data structure for transmitter.
Data structure for tripod.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Detector data structure.
Definition JDetector.hh:96
int getString() const
Get string number.
Definition JLocation.hh:135
Data structure for position in three dimensions.
General exception.
Definition JException.hh:24
int getID() const
Get identifier.
Definition JObjectID.hh:50
Utility class to parse command line options.
Definition JParser.hh:1698
Object writing to file.
General purpose class for object reading from a list of file names.
virtual bool hasNext() override
Check availability of next element.
counter_type getCounter() const
Get counter.
virtual const pointer_type & next() override
Get next element.
bool has(const T &value) const
Test whether given value is present.
static const int HYDROPHONE_DISABLE
Enable (disable) use of hydrophone if this status bit is 0 (1);.
static const int PIEZO_DISABLE
Enable (disable) use of piezo if this status bit is 0 (1);.
JPosition3D getPosition(const Vec &pos)
Get position.
static JEmitterID getEmitterID
Function object for emitter identification.
JContainer< std::vector< JTripod > > tripods_container
Definition JSydney.cc:79
JContainer< std::vector< JTransmitter > > transmitters_container
Definition JSydney.cc:81
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition JSydney.cc:80
static const JSoundVelocity getSoundVelocity(1541.0, -17.0e-3, -2000.0)
Function object for velocity of sound.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
JPosition3D getPiezoPosition()
Get relative position of piezo in optical module.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::vector< JHitW0 > buffer_type
hits
Definition JPerth.cc:70
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Detector file.
Definition JHead.hh:227
Acoustic emitter.
Definition JEmitter.hh:30
Match of two events considering overlap in time and position.
void merge(const JEvent &event)
Merge event.
int getID() const
Get emitter identifier.
Acoustic receiver.
Definition JReceiver.hh:30
Implementation for depth dependend velocity of sound.
JSoundVelocity & set(const double z0)
Set depth.
Time-of-arrival data from acoustic piezo sensor or hydrophone.
Definition JToA.hh:28
uint32_t DOMID
DAQ run number.
Definition JToA.hh:34
int32_t WAVEFORMID
DOM unique identifeir.
Definition JToA.hh:35
int32_t DETID
Definition JToA.hh:32
Acoustic transceiver.
JTransmission getTransmission(const JToA &data, const JAbstractSoundVelocity &V) const
Get transmission.
Auxiliary class to compare transmissions.
Auxiliary class to compare transmissions.
double Q
minimal quality if larger than one; else minimal normalised quality
double TMax_s
maximal difference between times of emission [s]
int numberOfHits
minimal number of hits to trigger event
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition JContainer.hh:42
Type list.
Definition JTypeList.hh:23
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
General purpose class for hash map of unique keys.
Definition JHashMap.hh:75
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
double getXmin() const
Get minimum value.
Definition JQuantile.hh:208
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133
long long int getCount() const
Get total count.
Definition JQuantile.hh:186
double getMean() const
Get mean value.
Definition JQuantile.hh:252
double getXmax() const
Get maximum value.
Definition JQuantile.hh:219