#include #include #include #include #include #include "TROOT.h" #include "TFile.h" #include "km3net-dataformat/definitions/module_status.hh" #include "JLang/JComparator.hh" #include "JDB/JToAshort.hh" #include "JDB/JSupport.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JDetectorSupportkit.hh" #include "JDetector/JTripod.hh" #include "JDetector/JTransmitter.hh" #include "JDetector/JHydrophone.hh" #include "JLang/JPredicate.hh" #include "JLang/JComparator.hh" #include "JTools/JRange.hh" #include "JTools/JQuantile.hh" #include "JTools/JHashMap.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JFileRecorder.hh" #include "JSupport/JMeta.hh" #include "JAcoustics/JEmitterID.hh" #include "JAcoustics/JTransceiver.hh" #include "JAcoustics/JTransmission.hh" #include "JAcoustics/JAcousticsToolkit.hh" #include "JAcoustics/JAcousticsSupportkit.hh" #include "JAcoustics/JTriggerParameters.hh" #include "JAcoustics/JEvent.hh" #include "JAcoustics/JEventOverlap.hh" #include "JAcoustics/JTriggerOutput.hh" #include "JAcoustics/JSupport.hh" #include "Jeep/JContainer.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Main program to trigger events in acoustic data. * * An acoustic event is based on a coincidence of estimated times of emission from the same emitter.\n * If the number of coincident times of emission exceeds a preset minimum, * the event is triggered and subsequently stored in the output file.\n * Note that an event counter is added which can be used to disambiguate * the different cycles from the same emitter within one so-called "super" event that is used for the global fit.\n * In this, a super event refers to the coincidence of events from a variety of emitters. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; typedef JTYPELIST::typelist typelist; typedef JContainer< vector > tripods_container; typedef JContainer< vector > transmitters_container; typedef JContainer< vector > hydrophones_container; JMultipleFileScanner inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); JTriggerParameters parameters; JFileRecorder outputFile; JSoundVelocity V = getSoundVelocity; // default sound velocity string detectorFile; tripods_container tripods; // tripods transmitters_container transmitters; // transmitters hydrophones_container hydrophones; // hydrophones double precision; int debug; try { JParser<> zap("Main program to trigger events in acoustic data."); zap['f'] = make_field(inputFile, "output of JConvertDB -q toashort"); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['@'] = make_field(parameters, "trigger parameters"); zap['o'] = make_field(outputFile, "output file") = "event.root"; zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised(); zap['a'] = make_field(detectorFile, "detector file"); 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['W'] = make_field(getEmitterID, "waveform identification data") = JPARSER::initialised(); zap['p'] = make_field(precision, "precision time-of-arrival") = 1.0e-6; zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } V.set(detector.getUTMZ()); JHashMap receivers; JHashMap emitters; for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) { if ((i->getFloor() == 0 && !i->has(HYDROPHONE_DISABLE)) || (i->getFloor() != 0 && !i->has(PIEZO_DISABLE))) { JPosition3D pos(0.0, 0.0, 0.0); if (i->getFloor() == 0) { // get relative position of hydrophone try { pos = getPosition(hydrophones.begin(), hydrophones.end(), make_predicate(&JHydrophone::getString, i->getString())); } catch(const exception&) { continue; // if no position available, discard hydrophone } } receivers[i->getID()] = JReceiver(i->getID(), i->getPosition() + pos, i->getT0() * 1.0e-9); } } 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&) { continue; // if no string available, discard transmitter } } const JRange unit(0.0, 1.0); outputFile.open(); if (!outputFile.is_open()) { FATAL("Error opening file " << outputFile << endl); } outputFile.put(JMeta(argc, argv)); outputFile.put(parameters); // statistics struct map_type : public map { long long int sum() const { long long int count = 0; for (const_iterator i = this->begin(); i != this->end(); ++i) { count += i->second.getCount(); } return count; } }; map_type number_of_entries; map_type number_of_aliens; map_type number_of_triggers; map_type number_of_events; // input data string oid; // detector identifier typedef vector buffer_type; // data type map > f1; // emitter -> receiver -> data while (inputFile.hasNext()) { if (inputFile.getCounter()%1000 == 0) { STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl); } JToAshort* parameters = inputFile.next(); if (oid == "") { oid = parameters->DETID; } if (oid != parameters->DETID) { // consistency check FATAL("Invalid detector identifier " << parameters->DETID << " != " << oid << endl); } try { const int id = getEmitterID(parameters->EMITTERID); number_of_entries[parameters->EMITTERID].put(1); if (emitters.has(id) && receivers.has(parameters->DOMID)) { const JTransceiver transceiver(emitters[id], receivers[parameters->DOMID]); f1[transceiver.emitter.getID()][transceiver.receiver.getID()].push_back(transceiver.getTransmission(*parameters, V)); } } catch(const exception&) { number_of_aliens[parameters->EMITTERID].put(1); } } STATUS(endl); for (map_type::const_iterator i = number_of_entries.begin(); i != number_of_entries.end(); ++i) { STATUS("number of entries: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl); } for (map_type::const_iterator i = number_of_aliens.begin(); i != number_of_aliens.end(); ++i) { STATUS("number of aliens: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl); } const JEventOverlap match(parameters.TMax_s); for (map >::iterator i = f1.begin(); i != f1.end(); ++i) { buffer_type buffer; for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) { // filter similar hits sort(receiver->second.begin(), receiver->second.end(), JTransmission::compare(precision)); buffer_type::iterator __end = unique(receiver->second.begin(), receiver->second.end(), JTransmission::equals(precision)); // selection based on quality for (buffer_type::const_iterator p = receiver->second.begin(); p != __end; ++p) { if (p->getQ() >= parameters.Q * (unit(parameters.Q) ? p->getW() : 1.0)) { buffer.push_back(*p); } } } sort(buffer.begin(), buffer.end()); // sort according time-of-emisson JQuantile Q; // monitor of trigger JEvent out; // temporary storage of event(s) JTriggerOutput data; // for (buffer_type::const_iterator p = buffer.begin(); p != buffer.end(); ++p) { buffer_type::const_iterator q = p; // basic correlator while (++q != buffer.end() && q->getToE() - p->getToE() <= parameters.TMax_s) {} Q.put(distance(p,q)); if (distance(p,q) >= parameters.numberOfHits) { STATUS("trigger: " << setw(2) << i->first << ' ' << setw(8) << number_of_triggers.sum() << '\r' << flush); DEBUG(endl); JEvent event(oid, number_of_triggers.sum(), i->first, p, q); number_of_triggers[i->first].put(event.size()); DEBUG("trigger: " << endl << event); if (out.empty()) { out = event; } else if (match(out, event)) { out.merge(event); } else { data.push_back(out); out = event; } } } if (!out.empty()) { data.push_back(out); // pending event } STATUS("number of toes: " << setw(3) << i->first << ' ' << setw(8) << buffer.size()); if (Q.getCount() > 0) { STATUS(' ' << FIXED(6,1) << Q.getMean() << ' ' << FIXED(6,1) << Q.getXmin() << ' ' << FIXED(6,1) << Q.getXmax()); } STATUS(endl); for (JTriggerOutput::const_iterator event = data.begin(); event != data.end(); ++event) { number_of_events[event->getID()].put(event->size()); outputFile.put(*event); } } STATUS(endl); for (map_type::const_iterator i = number_of_triggers.begin(); i != number_of_triggers.end(); ++i) { STATUS("number of triggers: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl); } for (map_type::const_iterator i = number_of_events.begin(); i != number_of_events.end(); ++i) { STATUS("number of events: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << ' ' << FIXED(7,3) << i->second.getMean() << endl); } JMultipleFileScanner io(inputFile); io >> outputFile; outputFile.close(); return 0; }