#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #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" /** * \file * * Application to make a global fit of the detector geometry to acoustic data.\n * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; typedef JContainer< vector > tripods_container; typedef JContainer< vector > transmitters_container; typedef JContainer< vector > hydrophones_container; typedef JContainer< set > disable_container; JMultipleFileScanner inputFile; JFileRecorder::typelist> outputFile; string detectorFile; JLimit_t& numberOfEvents = inputFile.getLimit(); JSoundVelocity V = getSoundVelocity; // default sound velocity tripods_container tripods; // tripods transmitters_container transmitters; // transmitters hydrophones_container hydrophones; // hydrophones JFitParameters parameters; // fit parameters disable_container disable; // disable tansmissions size_t threads; // number of parallel threads bool& squash = JFremantle::squash; int debug; try { JParser<> zap("Application to fit position calibration model to acoustic data."); zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]"); zap['o'] = make_field(outputFile) = ""; zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['a'] = make_field(detectorFile); zap['@'] = make_field(parameters) = JPARSER::initialised(); zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised(); zap['T'] = make_field(tripods, "tripod data"); zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised(); zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised(); zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised(); zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised(); zap['N'] = make_field(threads, "number of threads") = 1; zap['s'] = make_field(squash, "squash transmissions in output"); zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } ROOT::EnableThreadSafety(); JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } JHashMap receivers; JHashMap emitters; for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) { receivers[i->getID()] = i->getLocation(); } for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) { emitters[i->getID()] = JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition()); } for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) { try { emitters[i->getID()] = JEmitter(i->getID(), i->getPosition() + detector.getModule(i->getLocation()).getPosition()); } catch(const exception&) {} // if no module available, discard transmitter } V.set(detector.getUTMZ()); // sound velocity at detector depth JGeometry geometry(detector, hydrophones); JGlobalfit katoomba(geometry, V, parameters); JKatoomba::MATRIX_INVERSION = LDU_t; JKatoomba::debug = debug; JKatoomba::MAXIMUM_ITERATIONS = 10000; if (inputFile.size() > 1u) { // sort input files map zmap; for (const string& file_name : inputFile) { STATUS(file_name << '\r'); DEBUG(endl); for (JSingleFileScanner in(file_name, 1); in.hasNext(); ) { const JEvent* evt = in.next(); if (JFremantle::detid != evt->getDetectorID()) { FATAL("Invalid detector identifier " << evt->getDetectorID() << " != " << JFremantle::detid << endl); } if (!evt->empty()) { zmap[evt->begin()->getToE()] = file_name; } } } STATUS(endl); inputFile.clear(); for (map::const_iterator i = zmap.begin(); i != zmap.end(); ++i) { inputFile.push_back(i->second); } } if (outputFile.getFilename() != "" && outputFile.getFilename() != "--") { outputFile.open(); JFremantle::output = &outputFile; } /* outputFile.put(JMeta(argc, argv)); outputFile.put(parameters); */ try { JFremantle fremantle(geometry, V, parameters, threads, 2 * threads); typedef deque buffer_type; for (buffer_type zbuf; inputFile.hasNext(); ) { STATUS(inputFile.getFilename() << '\r'); DEBUG(endl); // read one file at a time for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) { const JEvent* evt = inputFile.next(); if (!evt->empty() && emitters.has(evt->getID())) { zbuf.push_back(*evt); } } sort(zbuf.begin(), zbuf.end()); // sort according first time-of-emission for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) { for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {} if (q == zbuf.end()) { if (inputFile.hasNext()) { zbuf.erase(zbuf.begin(), p); // remove processed data and continue reading break; } } JEvent::overlap(p, q, parameters.deadTime_s); // empty overlapping events if (getNumberOfEmitters(p,q) >= parameters.Nmin) { const JWeight getWeight(p, q); JFremantle::input_type data; for (buffer_type::iterator evt = p; evt != q; ++evt) { sort(evt->begin(), evt->end(), JKatoomba<>::compare); JEvent::iterator __end = unique(evt->begin(), evt->end(), make_comparator(&JTransmission::getID, JComparison::eq())); const JEmitter& emitter = emitters [evt->getID()]; const double weight = getWeight(evt->getID()); for (JEvent::const_iterator i = evt->begin(); i != __end; ++i) { if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 && disable.count(JTransmission_t(-1, i->getID())) == 0) { if (receivers.has(i->getID()) && geometry.hasLocation(receivers[i->getID()]) && i->getQ() >= parameters.Qmin * (parameters.Qmin <= 1.0 ? i->getW() : 1.0)) { data.push_back(JHit(emitter, distance(p,evt), receivers[i->getID()], i->getToA(), parameters.sigma_s, weight)); } } } } if (getMinimumNumberOfEmitters(data.begin(), data.end()) >= parameters.Nmin) { fremantle.enqueue(data); } } } } STATUS(endl); } catch(const exception& error) { FATAL("main " << error.what()); } /* JMultipleFileScanner io(inputFile); io >> outputFile; */ outputFile.close(); JFileOutputStream(3) << SCIENTIFIC(1,10) << JFremantle::Q.getMean(numeric_limits::max()) << endl; }