Jpp 20.0.0-rc.9-29-gccc23c492-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JFremantle.cc File Reference

Application to make a global fit of the detector geometry to acoustic data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <deque>
#include <algorithm>
#include <limits>
#include <type_traits>
#include <functional>
#include <future>
#include <mutex>
#include <thread>
#include <queue>
#include "TROOT.h"
#include "TFile.h"
#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"
#include "JLang/JComparison.hh"
#include "JLang/JFileStream.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JTransmitter.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JHydrophone.hh"
#include "JTools/JHashMap.hh"
#include "JTools/JRange.hh"
#include "JTools/JQuantile.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JGlobalfit.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JSuperEvt.hh"
#include "JAcoustics/JSuperEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JTransmission_t.hh"
#include "JAcoustics/JFremantle_t.hh"
#include "Jeep/JContainer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Application to make a global fit of the detector geometry to acoustic data.


Author
mdejong

Definition in file JFremantle.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 67 of file JFremantle.cc.

68{
69 using namespace std;
70 using namespace JPP;
71
72 typedef JContainer< vector<JTripod> > tripods_container;
73 typedef JContainer< vector<JTransmitter> > transmitters_container;
74 typedef JContainer< vector<JHydrophone> > hydrophones_container;
75 typedef JContainer< set<JTransmission_t> > disable_container;
76
77 JMultipleFileScanner<JEvent> inputFile;
78 JFileRecorder<JTYPELIST<JSuperEvt, JFitParameters, JMeta>::typelist> outputFile;
79 string detectorFile;
80 JLimit_t& numberOfEvents = inputFile.getLimit();
81 JSoundVelocity V = getSoundVelocity; // default sound velocity
82 tripods_container tripods; // tripods
83 transmitters_container transmitters; // transmitters
84 hydrophones_container hydrophones; // hydrophones
85 JFitParameters parameters; // fit parameters
86 disable_container disable; // disable tansmissions
87 size_t threads; // number of parallel threads
88 bool& squash = JFremantle::squash;
89 int debug;
90
91 try {
92
93 JParser<> zap("Application to fit position calibration model to acoustic data.");
94
95 zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
96 zap['o'] = make_field(outputFile) = "";
97 zap['n'] = make_field(numberOfEvents) = JLimit::max();
98 zap['a'] = make_field(detectorFile);
99 zap['@'] = make_field(parameters) = JPARSER::initialised();
100 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
101 zap['T'] = make_field(tripods, "tripod data");
102 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
103 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
104 zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
105 zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
106 zap['N'] = make_field(threads, "number of threads") = 1;
107 zap['s'] = make_field(squash, "squash transmissions in output");
108 zap['d'] = make_field(debug) = 1;
109
110 zap(argc, argv);
111 }
112 catch(const exception &error) {
113 FATAL(error.what() << endl);
114 }
115
116 ROOT::EnableThreadSafety();
117
118 JDetector detector;
119
120 try {
121 load(detectorFile, detector);
122 }
123 catch(const JException& error) {
124 FATAL(error);
125 }
126
127 JHashMap<int, JLocation> receivers;
128 JHashMap<int, JEmitter> emitters;
129
130 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
131 receivers[i->getID()] = i->getLocation();
132 }
133
134 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
135 emitters[i->getID()] = JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition());
136 }
137
138 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
139 try {
140 emitters[i->getID()] = JEmitter(i->getID(), i->getPosition() + detector.getModule(i->getLocation()).getPosition());
141 }
142 catch(const exception&) {} // if no module available, discard transmitter
143 }
144
145 V.set(detector.getUTMZ()); // sound velocity at detector depth
146
147 JGeometry geometry(detector, hydrophones);
148 JGlobalfit katoomba(geometry, V, parameters);
149
150 JKatoomba<JEstimator>::MATRIX_INVERSION = LDU_t;
151
152 JKatoomba<JGandalf>::debug = debug;
153 JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;
154
155 if (inputFile.size() > 1u) { // sort input files
156
157 map<double, string> zmap;
158
159 for (const string& file_name : inputFile) {
160
161 STATUS(file_name << '\r'); DEBUG(endl);
162
163 for (JSingleFileScanner<JEvent> in(file_name, 1); in.hasNext(); ) {
164
165 const JEvent* evt = in.next();
166
167 if (JFremantle::detid != evt->getDetectorID()) {
168 FATAL("Invalid detector identifier " << evt->getDetectorID() << " != " << JFremantle::detid << endl);
169 }
170
171 if (!evt->empty()) {
172 zmap[evt->begin()->getToE()] = file_name;
173 }
174 }
175 }
176 STATUS(endl);
177
178 inputFile.clear();
179
180 for (map<double, string>::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
181 inputFile.push_back(i->second);
182 }
183 }
184
185 if (outputFile.getFilename() != "" &&
186 outputFile.getFilename() != "--") {
187
188 outputFile.open();
189
190 JFremantle::output = &outputFile;
191 }
192 /*
193 outputFile.put(JMeta(argc, argv));
194 outputFile.put(parameters);
195 */
196
197 try {
198
199 JFremantle fremantle(geometry, V, parameters, threads, 2 * threads);
200
201 typedef deque<JEvent> buffer_type;
202
203 for (buffer_type zbuf; inputFile.hasNext(); ) {
204
205 STATUS(inputFile.getFilename() << '\r'); DEBUG(endl);
206
207 // read one file at a time
208
209 for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
210
211 const JEvent* evt = inputFile.next();
212
213 if (!evt->empty() && emitters.has(evt->getID())) {
214 zbuf.push_back(*evt);
215 }
216 }
217
218 sort(zbuf.begin(), zbuf.end()); // sort according first time-of-emission
219
220 for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
221
222 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}
223
224 if (q == zbuf.end()) {
225
226 if (inputFile.hasNext()) {
227
228 zbuf.erase(zbuf.begin(), p); // remove processed data and continue reading
229
230 break;
231 }
232 }
233
234 JEvent::overlap(p, q, parameters.deadTime_s); // empty overlapping events
235
236 if (getNumberOfEmitters(p,q) >= parameters.Nmin) {
237
238 const JWeight getWeight(p, q);
239
240 JFremantle::input_type data;
241
242 for (buffer_type::iterator evt = p; evt != q; ++evt) {
243
244 sort(evt->begin(), evt->end(), JKatoomba<>::compare);
245
246 JEvent::iterator __end = unique(evt->begin(), evt->end(), make_comparator(&JTransmission::getID, JComparison::eq()));
247
248 const JEmitter& emitter = emitters [evt->getID()];
249 const double weight = getWeight(evt->getID());
250
251 for (JEvent::const_iterator i = evt->begin(); i != __end; ++i) {
252
253 if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 &&
254 disable.count(JTransmission_t(-1, i->getID())) == 0) {
255
256 if (receivers.has(i->getID()) && geometry.hasLocation(receivers[i->getID()]) && i->getQ() >= parameters.Qmin * (parameters.Qmin <= 1.0 ? i->getW() : 1.0)) {
257
258 data.push_back(JHit(emitter,
259 distance(p,evt),
260 receivers[i->getID()],
261 i->getToA(),
262 parameters.sigma_s,
263 weight));
264 }
265 }
266 }
267 }
268
269 if (getMinimumNumberOfEmitters(data.begin(), data.end()) >= parameters.Nmin) {
270 fremantle.enqueue(data);
271 }
272 }
273 }
274 }
275 STATUS(endl);
276 }
277 catch(const exception& error) {
278 FATAL("main " << error.what());
279 }
280 /*
281 JMultipleFileScanner<JMeta> io(inputFile);
282
283 io >> outputFile;
284 */
285 outputFile.close();
286
287 JFileOutputStream(3) << SCIENTIFIC(1,10) << JFremantle::Q.getMean(numeric_limits<float>::max()) << endl;
288}
JACOUSTICS::JFremantle
Thread pool for global fits.
Definition JFremantle_t.hh:33
JACOUSTICS::JFremantle::input_type
std::vector< JHit > input_type
Definition JFremantle_t.hh:36
JACOUSTICS::getMinimumNumberOfEmitters
size_t getMinimumNumberOfEmitters(T __begin, T __end)
Get minimum number of emitters for any string in data.
Definition JEvtToolkit.hh:69
JACOUSTICS::getWeight
double getWeight(T __begin, T __end)
Get total weight of data points.
Definition JKatoomba_t.hh:60
JACOUSTICS::tripods_container
JContainer< std::vector< JTripod > > tripods_container
Definition JSydney.cc:79
JACOUSTICS::transmitters_container
JContainer< std::vector< JTransmitter > > transmitters_container
Definition JSydney.cc:81
JACOUSTICS::hydrophones_container
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition JSydney.cc:80
JACOUSTICS::getNumberOfEmitters
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
Definition JEvtToolkit.hh:45
JACOUSTICS::JEmitter
Acoustic emitter.
Definition JEmitter.hh:30
JACOUSTICS::JEvent
Acoustic event.
Definition JEvent.hh:43
JACOUSTICS::JEvent::getID
int getID() const
Get emitter identifier.
Definition JEvent.hh:151
JACOUSTICS::JEvent::getDetectorID
const int getDetectorID() const
Get detector identifier.
Definition JEvent.hh:129
JACOUSTICS::JFitParameters::Qmin
double Qmin
minimal quality transmission
Definition JFitParameters.hh:62
JACOUSTICS::JFitParameters::deadTime_s
double deadTime_s
dead time between events [s]
Definition JFitParameters.hh:70
JACOUSTICS::JFitParameters::Nmin
size_t Nmin
minimum number of emitters
Definition JFitParameters.hh:66
JACOUSTICS::JFitParameters::sigma_s
double sigma_s
time-of-arrival resolution [s]
Definition JFitParameters.hh:64
JACOUSTICS::JFitParameters::Tmax_s
double Tmax_s
time window to combine events [s]
Definition JFitParameters.hh:65
JACOUSTICS::JGlobalfit
Global fit of prameterised detector geometry to acoustics data.
Definition JGlobalfit.hh:29
JACOUSTICS::JHit
Acoustics hit.
Definition JHit.hh:35
JACOUSTICS::JKatoomba
Template definition of fit function of acoustic model.
Definition JKatoomba_t.hh:76
JACOUSTICS::JSoundVelocity
Implementation for depth dependend velocity of sound.
Definition JSoundVelocity.hh:37
JACOUSTICS::JSoundVelocity::set
JSoundVelocity & set(const double z0)
Set depth.
Definition JSoundVelocity.hh:62
JACOUSTICS::JTransmission_t
Acoustic transmission identifier.
Definition JTransmission_t.hh:28
JACOUSTICS::JWeight
Auxiliary data structure to unify weights of acoustics data according to the number of pings per emit...
Definition JAcousticsToolkit.hh:55
Generated by doxygen 1.10.0