#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "TROOT.h" #include "TFile.h" #include "km3net-dataformat/definitions/module_status.hh" #include "JLang/JPredicate.hh" #include "JLang/JComparator.hh" #include "JLang/JComparison.hh" #include "JLang/JSTDObjectWriter.hh" #include "JLang/JSTDObjectReader.hh" #include "JSystem/JStat.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 "JDetector/JLocationRouter.hh" #include "JFit/JGradient.hh" #include "JTools/JHashMap.hh" #include "JTools/JRange.hh" #include "JTools/JQuantile.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JAcoustics/JSoundVelocity.hh" #include "JAcoustics/JEmitter.hh" #include "JAcoustics/JAcousticsToolkit.hh" #include "JAcoustics/JHit.hh" #include "JAcoustics/JFitParameters.hh" #include "JAcoustics/JKatoomba_t.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 "JAcoustics/JPlatypus_t.hh" #include "Jeep/JTimer.hh" #include "Jeep/JeepToolkit.hh" #include "Jeep/JContainer.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace JACOUSTICS { using JLANG::JValueOutOfRange; using JEEP::JContainer; using JFIT::JParameter_t; using namespace JDETECTOR; typedef JContainer< std::vector > tripods_container; typedef JContainer< std::vector > hydrophones_container; typedef JContainer< std::vector > transmitters_container; /** * Script commands. */ static const char skip_t = '#'; //!< skip line static const std::string initialise_t = "initialise"; //!< initialise static const std::string fix_t = "fix"; //!< fix objects static const std::string string_t = "string"; //!< string static const std::string tripod_t = "tripod"; //!< tripod static const std::string stage_t = "stage"; //!< fit stage /** * Auxiliary data structure for handling of file names. */ struct JFilenames { std::string detector; //!< detector std::string tripod; //!< tripod std::string hydrophone; //!< hydrophone std::string transmitter; //!< transmitter }; /** * Auxiliary data structure for setup of complete system. */ struct JSetup { JDetector detector; //!< detector tripods_container tripods; //!< tripods struct : public hydrophones_container { /** * Check if there is a hydrophone on given string. * * \param id string identifier * \return true if hydrophone present; else false */ bool hasString(const int id) const { using namespace std; return (find_if(this->begin(), this->end(), make_predicate(&JHydrophone::getString, id)) != this->end()); } } hydrophones; //!< hydrophones transmitters_container transmitters; //!< transmitters }; /** * Main class for pre-calibration using acoustics data. */ struct JSydney { static constexpr double RADIUS_M = 1.0; //!< maximal horizontal distance between T-bar and emitter/hydrophone /** * List of object identifiers. */ struct ids_t : public std::set { /** * Default constructor. */ ids_t() {} /** * Copy constructor. * * \param buffer list of identifiers */ ids_t(const std::vector& buffer) : std::set(buffer.begin(), buffer.end()) {} /** * Difference constructor. * Make list of all object identifiers in A that are not in B. * * \param A list of identifiers * \param B list of identifiers */ ids_t(const ids_t& A, const ids_t& B) { std::set_difference(A.begin(), A.end(), B.begin(), B.end(), std::inserter(*this, this->begin())); } /** * Fix. * Keep list of object identifiers that are not in B. * * \param B list of identifiers */ void fix(const ids_t& B) { ids_t A; this->swap(A); std::set_difference(A.begin(), A.end(), B.begin(), B.end(), std::inserter(*this, this->begin())); } /** * Read identifiers from input stream * * \param in input stream * \param object identifiers * \return input stream */ friend inline std::istream& operator>>(std::istream& in, ids_t& object) { for (int id; in >> id; ) { object.insert(id); } if (!in.bad()) { in.clear(); } return in; } /** * Write identifiers to output stream * * \param out output stream * \param object identifiers * \return output stream */ friend inline std::ostream& operator<<(std::ostream& out, const ids_t& object) { for (const int id : object) { out << ' ' << id; } return out; } }; /** * Auxiliary data structure for group of lists of identifiers of to-be-fitted objects. */ struct fits_t { /** * Default constructor. */ fits_t() {} /** * Initialise. * * \param setup setup */ void initialise(const JSetup& setup) { using namespace JPP; strings = make_array(setup.detector .begin(), setup.detector .end(), &JModule ::getString); tripods = make_array(setup.tripods .begin(), setup.tripods .end(), &JTripod ::getID); hydrophones = make_array(setup.hydrophones .begin(), setup.hydrophones .end(), &JHydrophone ::getString); transmitters = make_array(setup.transmitters.begin(), setup.transmitters.end(), &JTransmitter::getString); } ids_t strings; //!< identifiers of strings ids_t tripods; //!< identifiers of tripods ids_t hydrophones; //!< identifiers of strings with hydrophone ids_t transmitters; //!< identifiers of strings with transmitter }; /** * Auxiliary class to edit (z) position of module. */ struct JModuleEditor : public JParameter_t { /** * Constructor. * * \param module module */ JModuleEditor(JModule& module) : module(module), direction(JGEOMETRY3D::JVector3Z_t) {} /** * Constructor. * * \param module module * \param direction direction */ JModuleEditor(JModule& module, const JVector3D& direction) : module(module), direction(direction) {} /** * Apply step. * * \param step step */ virtual void apply(const double step) override { using namespace JPP; module.add(direction * step); } private: JModule& module; JVector3D direction; }; /** * Auxiliary class to edit (x,y,z) position of string. */ struct JStringEditor : public JParameter_t { /** * Constructor. * * The option true and false correspond to all modules and optical modules only, respectively. * * \param setup setup * \param id string number * \param direction direction * \param option option */ JStringEditor(JSetup& setup, const int id, const JVector3D& direction, const bool option) : detector (setup.detector), direction(direction) { for (size_t i = 0; i != detector.size(); ++i) { if (detector[i].getString() == id && (detector[i].getFloor() != 0 || option)) { index.push_back(i); } } } /** * Apply step. * * \param step step */ virtual void apply(const double step) override { for (const auto i : index) { detector[i].add(direction * step); } } private: JDetector& detector; JVector3D direction; std::vector index; }; /** * Auxiliary class to edit length of Dyneema ropes. */ struct JDyneemaEditor : public JParameter_t { /** * Constructor. * * \param setup setup * \param id string number * \param z0 reference position */ JDyneemaEditor(JSetup& setup, const int id, const double z0 = 0.0) : detector (setup.detector), z0 (z0) { for (size_t i = 0; i != detector.size(); ++i) { if (detector[i].getString() == id && detector[i].getFloor() != 0) { index.push_back(i); } } } /** * Apply step. * * \param step step */ virtual void apply(const double step) override { for (const auto i : index) { JModule& module = detector[i]; if (step > 0.0) module.set(JVector3D(module.getX(), module.getY(), z0 + (module.getZ() - z0) * (1.0 + step))); else if (step < 0.0) module.set(JVector3D(module.getX(), module.getY(), z0 + (module.getZ() - z0) / (1.0 - step))); } } private: JDetector& detector; double z0; std::vector index; }; /** * Auxiliary class to edit (x,y,z) position of tripod. */ struct JTripodEditor : public JParameter_t { /** * Constructor. * * \param setup setup * \param id tripod identifier * \param direction direction */ JTripodEditor(JSetup& setup, const int id, const JVector3D& direction) : tripods (setup.tripods), direction(direction) { using namespace std; using namespace JPP; index = distance(tripods.begin(), find_if(tripods.begin(), tripods.end(), make_predicate(&JTripod::getID, id))); } /** * Apply step. * * \param step step */ virtual void apply(const double step) override { tripods[index].add(direction * step); } private: std::vector& tripods; JVector3D direction; size_t index; }; /** * Auxiliary class to edit orientation of anchor. */ struct JAnchorEditor : public JParameter_t { /** * Constructor. * * \param setup setup * \param id string identifier */ JAnchorEditor(JSetup& setup, const int id) : hydrophones (setup.hydrophones), transmitters(setup.transmitters) { using namespace std; using namespace JPP; index[0] = distance(hydrophones .begin(), find_if(hydrophones .begin(), hydrophones .end(), make_predicate(&JHydrophone ::getString, id))); index[1] = distance(transmitters.begin(), find_if(transmitters.begin(), transmitters.end(), make_predicate(&JTransmitter::getString, id))); } /** * Apply step. * * \param step step */ virtual void apply(const double step) override { using namespace JPP; const JRotation3Z R(step); if (index[0] != hydrophones .size()) { hydrophones [index[0]].rotate(R); } if (index[1] != transmitters.size()) { transmitters[index[1]].rotate(R); } } private: std::vector& hydrophones; std::vector& transmitters; size_t index[2]; }; /** * Extended data structure for parameters of stage. */ struct JParameters_t : public JFitParameters { /** * Default constuctor. */ JParameters_t() : Nmax (std::numeric_limits::max()), Nextra (0), epsilon(1.0e-4), debug (3) {} /** * Read parameters from input stream * * \param in input stream * \param object parameters * \return input stream */ friend inline std::istream& operator>>(std::istream& in, JParameters_t& object) { object = JParameters_t(); in >> object.option >> object.mestimator >> object.sigma_s >> object.stdev >> object.Nextra; if (in) { for (double value; in >> value; ) { object.steps.push_back(value); } if (!in.bad()) { in.clear(); } } return in; } /** * Write parameters to output stream * * \param out output stream * \param object parameters * \return output stream */ friend inline std::ostream& operator<<(std::ostream& out, const JParameters_t& object) { using namespace std; out << setw(2) << object.option << ' ' << setw(2) << object.mestimator << ' ' << SCIENTIFIC(9,3) << object.sigma_s << ' ' << SCIENTIFIC(9,3) << object.stdev << ' ' << setw(3) << object.Nextra; for (const double value : object.steps) { out << ' ' << FIXED(9,5) << value; } return out; } size_t Nmax; size_t Nextra; double epsilon; int debug; std::vector steps; }; /** * Constructor. * * \param filenames file names * \param V sound velocity * \param threads threads * \param debug debug */ JSydney(const JFilenames& filenames, const JSoundVelocity& V, const size_t threads, const int debug) : filenames(filenames), V(V), threads(threads), debug(debug) { load(filenames.detector, setup.detector); setup.tripods.load(filenames.tripod.c_str()); if (filenames.hydrophone != "") { setup.hydrophones .load(filenames.hydrophone .c_str()); } if (filenames.transmitter != "") { setup.transmitters.load(filenames.transmitter.c_str()); } for (JDetector::const_iterator i = setup.detector.begin(); i != setup.detector.end(); ++i) { receivers[i->getID()] = i->getLocation(); } // detach PMTs detector = setup.detector; for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) { module->clear(); } router.reset(new JLocationRouter(setup.detector)); fits.initialise(setup); this->V.set(setup.detector.getUTMZ()); // sound velocity at detector depth JKatoomba::debug = debug; JKatoomba::MAXIMUM_ITERATIONS = 10000; JKatoomba::MATRIX_INVERSION = LDU_t; ROOT::EnableThreadSafety(); } /** * Auxiliary data structure for decomposed string. */ struct string_type : public std::vector // optical modules { /** * Add module. * * \param module module */ void push_back(const JModule& module) { if (module.getFloor() == 0) this->base = module; else std::vector::push_back(module); } JModule base; // base module }; /** * Auxiliary data structure for detector with decomposed strings. */ struct detector_type : public std::map { /** * Add module. * * \param module module */ void push_back(const JModule& module) { (*this)[module.getString()].push_back(module); } }; /** * Fit procedure to determine the positions of tripods and transmitters using strings that are fixed. * * \param parameters parameters */ void stage_0(const JParameters_t& parameters) { using namespace std; using namespace JPP; this->parameters = parameters; JDetector A; // old strings detector_type B; // new strings with transmitter -> fit only base module JDetector C; // new strings w/o transmitter -> discard completely from fit for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) { if (fits.strings .count(module->getString()) == 0) A.push_back(*module); else if (fits.transmitters.count(module->getString()) != 0) B.push_back(*module); else C.push_back(*module); } setup.detector.swap(A); for (const auto& element : B) { setup.detector.push_back(element.second.base); } router.reset(new JLocationRouter(setup.detector)); JGradient fit(parameters.Nmax, parameters.Nextra, parameters.epsilon, parameters.debug); for (const int i : fits.tripods) { fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "tripod.x: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3X_t), parameters.steps[0])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "tripod.y: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3Y_t), parameters.steps[0])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "tripod.z: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3Z_t), parameters.steps[0])); } for (const int i : fits.transmitters) { fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "transmitter.x: " << RIGHT(4) << i), new JStringEditor(setup, i, JVector3X_t, true), parameters.steps[0])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "transmitter.y: " << RIGHT(4) << i), new JStringEditor(setup, i, JVector3Y_t, true), parameters.steps[0])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "transmitter.z: " << RIGHT(4) << i), new JStringEditor(setup, i, JVector3Z_t, true), parameters.steps[0])); } const double chi2 = fit(*this); for (auto& element : B) { const JModule& base = setup.detector.getModule(router->getAddress(JLocation(element.second.base.getString(),0))); const JVector3D pos = base.getPosition() - element.second.base.getPosition(); for (string_type::iterator module = element.second.begin(); module != element.second.end(); ++module) { module->add(pos); setup.detector.push_back(*module); } } copy(C.begin(), C.end(), back_inserter(setup.detector)); sort(setup.detector.begin(), setup.detector.end(), make_comparator(&JModule::getLocation)); router.reset(new JLocationRouter(setup.detector)); STATUS("detector: " << FIXED(9,4) << chi2 << endl); } /** * Fit procedure to determine the positions of the strings and tripods. * * \param parameters parameters */ void stage_a(const JParameters_t& parameters) { using namespace std; using namespace JPP; this->parameters = parameters; JGradient fit(parameters.Nmax, parameters.Nextra, parameters.epsilon, parameters.debug); for (const int i : fits.strings) { fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "string.x: " << RIGHT(4) << i), new JStringEditor(setup, i, JVector3X_t, true), parameters.steps[0])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "string.y: " << RIGHT(4) << i), new JStringEditor(setup, i, JVector3Y_t, true), parameters.steps[0])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "string.z: " << RIGHT(4) << i), new JStringEditor(setup, i, JVector3Z_t, false), parameters.steps[0])); } for (const int i : fits.tripods) { fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "tripod.x: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3X_t), parameters.steps[1])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "tripod.y: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3Y_t), parameters.steps[1])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "tripod.z: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3Z_t), parameters.steps[1])); } for (const int i : ids_t(fits.hydrophones, fits.transmitters)) { fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "anchor.R: " << RIGHT(4) << i), new JAnchorEditor(setup, i), parameters.steps[0] / RADIUS_M)); } for (const int i : fits.transmitters) { JModule& module = setup.detector.getModule(router->getAddress(JLocation(i,0))); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "anchor.R: " << RIGHT(4) << i), new JAnchorEditor(setup, i), parameters.steps[0] / RADIUS_M)); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "anchor.z: " << RIGHT(4) << i), new JModuleEditor(module), parameters.steps[0])); } const double chi2 = fit(*this); STATUS("detector: " << FIXED(9,4) << chi2 << endl); } /** * Fit procedure to determine the stretching and z-positions of individual strings. * * \param parameters parameters */ void stage_b(const JParameters_t& parameters) { using namespace std; using namespace JPP; this->parameters = parameters; map buffer; double z0 = 0.0; for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) { buffer[module->getString()].push_back(*module); if (module->getZ() > z0) { z0 = module->getZ(); } } JDetector tx; for (transmitters_container::iterator i = setup.transmitters.begin(); i != setup.transmitters.end(); ++i) { try { tx.push_back(router->getModule(i->getLocation())); } catch(const exception&) {} } for (const int i : fits.strings) { setup.detector.swap(buffer[i]); copy(tx.begin(), tx.end(), back_inserter(setup.detector)); router.reset(new JLocationRouter(setup.detector)); JGradient fit(parameters.Nmax, parameters.Nextra, parameters.epsilon, parameters.debug); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "string.M: " << RIGHT(4) << i), new JDyneemaEditor(setup, i, z0), parameters.steps[0])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "string.z: " << RIGHT(4) << i), new JStringEditor (setup, i, JVector3Z_t, false), parameters.steps[1])); const double chi2 = fit(*this); STATUS("string: " << setw(4) << i << ' ' << FIXED(9,4) << chi2 << endl); buffer[i].clear(); copy_if(setup.detector.begin(), setup.detector.end(), back_inserter(buffer[i]), make_predicate(&JModule::getString, i)); } setup.detector.clear(); for (const auto& element : buffer) { copy(element.second.begin(), element.second.end(), back_inserter(setup.detector)); } router.reset(new JLocationRouter(setup.detector)); } /** * Fit procedure to determine the z-positions of the modules. * * \param parameters parameters */ void stage_c(const JParameters_t& parameters) { using namespace std; using namespace JPP; this->parameters = parameters; JGradient fit(parameters.Nmax, parameters.Nextra, parameters.epsilon, parameters.debug); for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) { if (fits.strings.count(module->getString()) != 0 && module->getFloor() != 0) { fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "module.z: " << right << module->getLocation()), new JModuleEditor(*module), parameters.steps[0])); } } const double chi2 = fit(*this); STATUS("detector: " << FIXED(9,4) << chi2 << endl); } /** * Fit procedure to determine the (x,y,z) positions of the modules. * * \param parameters parameters */ void stage_e(const JParameters_t& parameters) { using namespace std; using namespace JPP; this->parameters = parameters; JGradient fit(parameters.Nmax, parameters.Nextra, parameters.epsilon, parameters.debug); for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) { if (fits.strings.count(module->getString()) != 0 && module->getFloor() != 0) { fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "module.x: " << right << module->getLocation()), new JModuleEditor(*module, JVector3X_t), parameters.steps[0])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "module.y: " << right << module->getLocation()), new JModuleEditor(*module, JVector3Y_t), parameters.steps[0])); fit.push_back(JModifier_t(MAKE_STRING(LEFT(16) << "module.z: " << right << module->getLocation()), new JModuleEditor(*module, JVector3Z_t), parameters.steps[0])); } } const double chi2 = fit(*this); STATUS("detector: " << FIXED(9,4) << chi2 << endl); } /** * Get chi2. * * \param option option * \return chi2/NDF */ double operator()(const int option) const { using namespace std; using namespace JPP; const JGeometry geometry(setup.detector, setup.hydrophones); JHashMap emitters; for (tripods_container::const_iterator i = setup.tripods.begin(); i != setup.tripods.end(); ++i) { { emitters[i->getID()] = JEmitter(i->getID(), i->getUTMPosition() - setup.detector.getUTMPosition()); } } for (transmitters_container::const_iterator i = setup.transmitters.begin(); i != setup.transmitters.end(); ++i) { try { emitters[i->getID()] = JEmitter(i->getID(), i->getPosition() + router->getModule(i->getLocation()).getPosition()); } catch(const exception&) {} // if no module available, discard transmitter } if (option == 0 || // step wise improvement of the chi2 option == 1) { // evaluation of the chi2 before the determination of the gradient of the chi2 this->output.clear(); JSTDObjectWriter out(this->output); // write data for subsequent use JFremantle::output = (option == 1 ? &out : NULL); { JFremantle fremantle(geometry, V, parameters, threads, 2 * threads); for (const input_type& superevt : input) { const JWeight getWeight(superevt.begin(), superevt.end()); set counter; vector data; for (input_type::const_iterator evt = superevt.begin(); evt != superevt.end(); ++evt) { if (emitters.has(evt->getID())) { const JEmitter& emitter = emitters [evt->getID()]; const double weight = getWeight(evt->getID()); for (JEvent::const_iterator i = evt->begin(); i != evt->end(); ++i) { if (geometry.hasLocation(receivers[i->getID()])) { data.push_back(JHit(emitter, distance(superevt.begin(), evt), receivers[i->getID()], i->getToA(), parameters.sigma_s, weight)); counter.insert(evt->getID()); } } } } if (counter.size() >= parameters.Nmin) { fremantle.enqueue(data); } } } return JFremantle::Q.getMean(numeric_limits::max()); } else if (option == 2) { // evaluation of the derivative of the chi2 to each fit parameter { JSTDObjectReader in(this->output); JPlatypus platypus(geometry, emitters, V, parameters, in, threads); } return JPlatypus::Q.getMean(numeric_limits::max()); } else { return numeric_limits::max(); } } /** * Run. * * \param script steering script */ void run(const std::string& script) { using namespace std; using namespace JPP; ifstream in(script.c_str()); while (in) { string buffer, key; if (getline(in, buffer)) { if (buffer.empty() || buffer[0] == skip_t) { continue; } istringstream is(buffer); is >> key; if (key == initialise_t) { // set object identifiers fits.initialise(setup); } else if (key == fix_t) { // fix object identifiers string type; // type of object ids_t id; // identifiers if (is >> type >> id) { if (type == string_t) { fits.strings .fix(id); fits.hydrophones .fix(id); fits.transmitters.fix(id); } else if (type == tripod_t) { fits.tripods .fix(id); } else { THROW(JValueOutOfRange, "Invalid type <" << type << ">"); } } } else if (key == stage_t) { // stage string stage; JParameters_t input; if (is >> stage >> input) { STATUS("stage " << setw(3) << stage << " {" << input << "}" << endl); JTimer timer; timer.start(); ofstream out(MAKE_CSTRING("stage-" << stage << ".log")); { JRedirectStream redirect(cout, out); switch (stage[stage.size() - 1]) { case '0': stage_0(input); break; case 'a': case 'A': stage_a(input); break; case 'b': case 'B': stage_b(input); break; case 'c': case 'C': stage_c(input); break; case 'e': case 'E': stage_e(input); break; default: THROW(JValueOutOfRange, "Invalid stage <" << stage << ">"); break; } } out.close(); store(stage); store(); timer.stop(); STATUS("Elapsed time " << FIXED(12,3) << timer.usec_wall * 1.0e-6 << " s." << endl); } } else { THROW(JValueOutOfRange, "Invalid key <" << key << ">"); } } } in.close(); } /** * Store data in given directory. * * \param dir directory */ void store(const std::string& dir = ".") { using namespace JPP; if (getFileStatus(dir.c_str()) || (mkdir(dir.c_str(), S_IRWXU | S_IRWXG) != -1)) { // attach PMTs for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) { module->set(router->getModule(module->getLocation()).getPosition()); } JDETECTOR::store(getFilename(dir, filenames.detector), detector); setup.tripods.store(getFilename(dir, filenames.tripod).c_str()); if (filenames.hydrophone != "") { setup.hydrophones .store(getFilename(dir, filenames.hydrophone) .c_str()); } if (filenames.transmitter != "") { setup.transmitters.store(getFilename(dir, filenames.transmitter).c_str()); } } else { THROW(JValueOutOfRange, "Invalid directory <" << dir << ">"); } } /** * Get list of identifiers of receivers. * * \return list of identifiers */ ids_t getReceivers() const { ids_t buffer; for (JDetector::const_iterator i = setup.detector.begin(); i != setup.detector.end(); ++i) { if ((i->getFloor() != 0 && !i->has(PIEZO_DISABLE)) || (setup.hydrophones.hasString(i->getString()) && !i->has(HYDROPHONE_DISABLE))) { buffer.insert(i->getID()); } } return buffer; } /** * Get list of identifiers of emitters. * * \return list of identifiers */ ids_t getEmitters() const { using namespace std; ids_t buffer; for (tripods_container::const_iterator i = setup.tripods.begin(); i != setup.tripods.end(); ++i) { buffer.insert(i->getID()); } for (transmitters_container::const_iterator i = setup.transmitters.begin(); i != setup.transmitters.end(); ++i) { try { const JModule& module = router->getModule(i->getLocation()); if (!module.has(TRANSMITTER_DISABLE)) { buffer.insert(i->getID()); } } catch(const exception&) {} } return buffer; } JFilenames filenames; JSoundVelocity V; size_t threads; int debug; JSetup setup; fits_t fits; JTOOLS::JHashMap receivers; std::unique_ptr router; typedef std::vector input_type; std::vector input; private: mutable std::vector output; JDetector detector; //!< PMTs JFitParameters parameters; }; } /** * \file * * Application to perform acoustic pre-calibration. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; typedef JContainer< std::set > disable_container; JMultipleFileScanner inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); JFilenames filenames; // file names JFitParameters parameters; // fit parameters string script; // script file JSoundVelocity V = getSoundVelocity; // default sound velocity disable_container disable; // disable tansmissions size_t threads; // number of threads int debug; try { JParser<> zap("Application to perform acoustic pre-calibration."); zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]"); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['a'] = make_field(filenames.detector); zap['@'] = make_field(parameters) = JPARSER::initialised(); zap['s'] = make_field(script, "steering script"); zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised(); zap['T'] = make_field(filenames.tripod, "tripod file"); zap['Y'] = make_field(filenames.transmitter, "transmitter file") = JPARSER::initialised(); zap['H'] = make_field(filenames.hydrophone, "hydrophone file") = 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['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } if (threads == 0) { FATAL("Invalid number of threads " << threads << endl); } JSydney sydney(filenames, V, threads, debug); const JSydney::ids_t receivers = sydney.getReceivers(); const JSydney::ids_t emitters = sydney.getEmitters(); typedef vector buffer_type; buffer_type zbuf; while (inputFile.hasNext()) { const JEvent* evt = inputFile.next(); if (emitters.count(evt->getID())) { zbuf.push_back(*evt); } } const JRange unit(0.0, 1.0); 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; ) {} JEvent::overlap(p, q, parameters.deadTime_s); // empty overlapping events JSydney::input_type buffer; 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())); for (JEvent::iterator i = evt->begin(); i != __end; ) { if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 && disable.count(JTransmission_t(-1, i->getID())) == 0) { if (receivers.count(i->getID()) && i->getQ() >= parameters.Qmin * (unit(parameters.Qmin) ? i->getW() : 1.0)) { ++i; continue; } } iter_swap(i, --__end); } buffer.push_back(JEvent(evt->getOID(), buffer.size(), evt->getID(), evt->begin(), __end)); } if (getNumberOfEmitters(buffer.begin(), buffer.end()) >= parameters.Nmin) { sydney.input.push_back(buffer); } } sydney.run(script); }