Jpp 20.0.0-rc.9-29-gccc23c492-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JKatoomba.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 "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"
#include "JLang/JComparison.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JROOTClassSelector.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 "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/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JTransmission_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 JKatoomba.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 60 of file JKatoomba.cc.

61{
62 using namespace std;
63 using namespace JPP;
64
65 typedef JContainer< vector<JTripod> > tripods_container;
66 typedef JContainer< vector<JTransmitter> > transmitters_container;
67 typedef JContainer< vector<JHydrophone> > hydrophones_container;
68 typedef JContainer< set<JTransmission_t> > disable_container;
69
70 JMultipleFileScanner<JEvent> inputFile;
71 JFileRecorder<JTYPELIST<JEvt, JEvent, JFitParameters, JMeta, TH1D, TH2D>::typelist> outputFile;
72 string detectorFile;
73 JLimit_t& numberOfEvents = inputFile.getLimit();
74 JSoundVelocity V = getSoundVelocity; // default sound velocity
75 tripods_container tripods; // tripods
76 transmitters_container transmitters; // transmitters
77 hydrophones_container hydrophones; // hydrophones
78 JFitParameters parameters; // fit parameters
79 disable_container disable; // disable tansmissions
80 JROOTClassSelection selection = getROOTClassSelection<JTYPELIST<JEvent>::typelist>();
81 JRange<double> utc;
82 bool squash;
83 Long64_t autoflush;
84 int debug;
85
86 try {
87
88 JParser<> zap("Application to fit position calibration model to acoustic data.");
89
90 zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
91 zap['o'] = make_field(outputFile);
92 zap['n'] = make_field(numberOfEvents) = JLimit::max();
93 zap['a'] = make_field(detectorFile);
94 zap['@'] = make_field(parameters) = JPARSER::initialised();
95 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
96 zap['T'] = make_field(tripods, "tripod data");
97 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
98 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
99 zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
100 zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
101 zap['C'] = make_field(selection,
102 "Precede name of data structure by a '+' or '-' "
103 "to add or remove data types in the output, respectively."
104 "\nROOT wildcards are accepted.") = JPARSER::initialised();
105 zap['u'] = make_field(utc, "select UTC time range") = JRange<double>();
106 zap['q'] = make_field(squash, "squash meta data");
107 zap['R'] = make_field(autoflush) = 5000;
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
117 JDetector detector;
118
119 try {
120 load(detectorFile, detector);
121 }
122 catch(const JException& error) {
123 FATAL(error);
124 }
125
126 const floor_range range = getRangeOfFloors(detector);
127
128 JHashMap<int, JLocation> receivers;
129 JHashMap<int, JEmitter> emitters;
130
131 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
132 receivers[i->getID()] = i->getLocation();
133 }
134
135 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
136 emitters[i->getID()] = JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition());
137 }
138
139 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
140 try {
141 emitters[i->getID()] = JEmitter(i->getID(), i->getPosition() + detector.getModule(i->getLocation()).getPosition());
142 }
143 catch(const exception&) {} // if no module available, discard transmitter
144 }
145
146 V.set(detector.getUTMZ()); // sound velocity at detector depth
147
148 JGeometry geometry(detector, hydrophones);
149 JGlobalfit katoomba(geometry, V, parameters);
150
151 JKatoomba<JGandalf>::debug = debug;
152 JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;
153
154 DEBUG(geometry);
155
156 typedef vector<JHit> data_type;
157
158 TH1D h0("chi2/NDF", NULL, 50000, 0.0, 1000.0);
159 TH1D h1("h1", NULL, 51, -0.5, 50.5);
160 TH1D hn("hn", NULL, 100, 0.0, 6.0);
161
162 JManager<int, TH2D> HA(new TH2D("ha[%]", NULL,
163 getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
164 range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
165
166 JManager<int, TH2D> HB(new TH2D("hb[%]", NULL,
167 getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
168 range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
169
170 for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
171 HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
172 HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
173 }
174
175 if (inputFile.size() > 1u) { // sort input files
176
177 map<double, string> zmap;
178
179 for (const string& file_name : inputFile) {
180
181 STATUS(file_name << '\r'); DEBUG(endl);
182
183 for (JSingleFileScanner<JEvent> in(file_name, 1); in.hasNext(); ) {
184
185 const JEvent* evt = in.next();
186
187 if (detector.getID() != evt->getDetectorID()) {
188 FATAL("Invalid detector identifier " << evt->getDetectorID() << " != " << evt->getDetectorID() << endl);
189 }
190
191 if (!evt->empty()) {
192 zmap[evt->begin()->getToE()] = file_name;
193 }
194 }
195 }
196 STATUS(endl);
197
198 inputFile.clear();
199
200 for (map<double, string>::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
201 inputFile.push_back(i->second);
202 }
203 }
204
205 outputFile.open();
206
207 outputFile.put(JMeta(argc, argv));
208 outputFile.put(parameters);
209
210 try { // limit RAM usage
211 dynamic_cast<JTreeWriterObjectOutput<JEvent>&>(*outputFile).getTreeWriter().SetAutoFlush(autoflush);
212 dynamic_cast<JTreeWriterObjectOutput<JEvent>&>(*outputFile).getTreeWriter().SetCacheSize();
213 }
214 catch(const exception&) {}
215
216 int counter[] = { 0, 0 };
217
218 typedef deque<JEvent> buffer_type;
219
220 for (buffer_type zbuf; inputFile.hasNext(); ) {
221
222 STATUS(inputFile.getFilename() << '\r'); DEBUG(endl);
223
224 // read one file at a time
225
226 for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
227
228 const JEvent* evt = inputFile.next();
229
230 if (!evt->empty() && emitters.has(evt->getID())) {
231 if (utc(evt->begin()->getToA()) && utc(evt->rbegin()->getToA())) {
232 zbuf.push_back(*evt);
233 }
234 }
235 }
236
237 sort(zbuf.begin(), zbuf.end()); // sort according first time-of-emission
238
239 for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
240
241 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}
242
243 if (q == zbuf.end()) {
244
245 if (inputFile.hasNext()) {
246
247 zbuf.erase(zbuf.begin(), p); // remove processed data and continue reading
248
249 break;
250 }
251 }
252
253 JEvent::overlap(p, q, parameters.deadTime_s); // empty overlapping events
254
255 if (getNumberOfEmitters(p,q) >= parameters.Nmin) {
256
257 h1.Fill((double) getNumberOfEmitters(p,q));
258
259 const JWeight getWeight(p, q);
260
261 data_type data;
262
263 for (buffer_type::iterator evt = p; evt != q; ++evt) {
264
265 sort(evt->begin(), evt->end(), JKatoomba<>::compare);
266
267 JEvent::iterator __end = unique(evt->begin(), evt->end(), make_comparator(&JTransmission::getID, JComparison::eq()));
268
269 const JEmitter& emitter = emitters [evt->getID()];
270 const double weight = getWeight(evt->getID());
271
272 for (JEvent::const_iterator i = evt->begin(); i != __end; ++i) {
273
274 if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 &&
275 disable.count(JTransmission_t(-1, i->getID())) == 0) {
276
277 if (receivers.has(i->getID()) && geometry.hasLocation(receivers[i->getID()]) && i->getQ() >= parameters.Qmin * (parameters.Qmin <= 1.0 ? i->getW() : 1.0)) {
278
279 data.push_back(JHit(emitter,
280 distance(p,evt),
281 receivers[i->getID()],
282 i->getToA(),
283 parameters.sigma_s,
284 weight));
285 }
286 }
287 }
288 }
289
290 if (getMinimumNumberOfEmitters(data.begin(), data.end()) >= parameters.Nmin) {
291
292 for (data_type::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
293 HA[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
294 }
295
296 // fit
297
298 const auto result = katoomba(data.begin(), data.end());
299
300 for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
301 HB[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
302 }
303
304 if (debug >= debug_t) {
305
306 cout << "result:" << ' '
307 << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl
308 << result.value << endl;
309
310 for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
311 cout << "hit: " << *hit << ' ' << FIXED(9,3) << katoomba.evaluator(result.value, *hit) << endl;
312 }
313 }
314
315 hn.Fill(log10(katoomba.gandalf.numberOfIterations));
316 h0.Fill(result.chi2 / result.ndf);
317
318 // store results
319
320 if ((parameters.chi2perNDF > 0.0 && result.chi2 / result.ndf <= +parameters.chi2perNDF) ||
321 (parameters.chi2perNDF < 0.0 && result.chi2 / result.ndf >= -parameters.chi2perNDF)) {
322
323 const JEvt evt = getEvt(JHead(detector.getID(),
324 result.getTimeRange(),
325 data .size(),
326 result.size(),
327 result.value.getN(),
328 result.ndf,
329 result.chi2,
330 katoomba.gandalf.numberOfIterations),
331 result.value);
332
333 outputFile.put(evt);
334
335 if (selection.is_valid<JEvent>()) {
336
337 for (buffer_type::iterator i = p; i != q; ++i) {
338
339 JEvent out(*i); // event with fitted times of emission
340
341 const double toe = result.value.emission[JEKey(i->getID(), distance(p,i))].t1;
342
343 for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
344 hit->setToE(toe);
345 }
346
347 if (!out.empty()) {
348 outputFile.put(out);
349 }
350 }
351 }
352
353 counter[0] += 1;
354
355 } else {
356
357 WARNING(endl << "Event not written - chi2/NDF " << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl);
358
359 counter[1] += 1;
360 }
361
362 } else {
363
364 WARNING(endl << "Event not accepted - minimum number of emitters " << getMinimumNumberOfEmitters(data.begin(), data.end()) << endl);
365
366 counter[1] += 1;
367 }
368 }
369 }
370 }
371 STATUS(endl);
372
373 STATUS("Number of events written / rejected: " << counter[0] << " / " << counter[1] << endl);
374
375 if (!squash) {
376
377 JMultipleFileScanner<JMeta> io(inputFile);
378
379 io >> outputFile;
380 }
381
382 outputFile.put(h0);
383 outputFile.put(h1);
384 outputFile.put(hn);
385
386 for (JManager<int, TH2D>::const_iterator i = HA.begin(); i != HA.end(); ++i) { outputFile.put(*(i->second)); }
387 for (JManager<int, TH2D>::const_iterator i = HB.begin(); i != HB.end(); ++i) { outputFile.put(*(i->second)); }
388
389 outputFile.close();
390}
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::JEKey
Emitter key.
Definition JEKey.hh:36
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::JEvt
Acoustic event fit.
Definition JEvt.hh:310
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::JFitParameters::chi2perNDF
double chi2perNDF
maximal chi2/NDF to store event
Definition JFitParameters.hh:69
JACOUSTICS::JGlobalfit
Global fit of prameterised detector geometry to acoustics data.
Definition JGlobalfit.hh:29
JACOUSTICS::JHead
Acoustic event header.
Definition JEvt.hh:157
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
JACOUSTICS::getEvt
Auxiliary data structure to convert model to event.
Definition JEvtToolkit.hh:155
Generated by doxygen 1.10.0