#include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "JDB/JToAshort.hh" #include "JDB/JSupport.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JTripod.hh" #include "JDetector/JTransmitter.hh" #include "JDetector/JHydrophone.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JTreeScanner.hh" #include "JTools/JHashMap.hh" #include "JTools/JRange.hh" #include "JLang/JComparator.hh" #include "JLang/JComparison.hh" #include "JROOT/JManager.hh" #include "JROOT/JRootToolkit.hh" #include "JAcoustics/JEmitterID.hh" #include "JAcoustics/JSoundVelocity.hh" #include "JAcoustics/JEmitter.hh" #include "JAcoustics/JReceiver.hh" #include "JAcoustics/JTransceiver.hh" #include "JAcoustics/JAcousticsToolkit.hh" #include "JAcoustics/JAcousticsSupportkit.hh" #include "JAcoustics/JHit.hh" #include "JAcoustics/JEvent.hh" #include "JAcoustics/JSupport.hh" #include "Jeep/JContainer.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { /** * Auxiliary data structure for organistation of histograms. */ struct key_type : public std::pair { /** * Constructor. * * \param first first value * \param second second value */ key_type(const int first, const int second) : std::pair(first, second) {} /** * Read key from input. * * \param in input stream * \param key key * \return input stream */ friend inline std::istream& operator>>(std::istream& in, key_type& key) { in >> key.first; in >> key.second; return in; } /** * Write key to output. * * \param out output stream * \param key key * \return output stream */ friend inline std::ostream& operator<<(std::ostream& out, const key_type& key) { out << key.first; out << '.'; out << key.second; return out; } /** * Less=than operator. * * \param first first key * \param second second key * \return true if first key less than second key; else false */ friend inline bool operator<(const key_type& first, const key_type& second) { if (first.first == second.first) return first.second < second.second; else return first.first < second.first; } }; } /** * \file * * Example program to plot hydrophone data. * \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; JMultipleFileScanner inputFile; string detectorFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string outputFile; JSoundVelocity V = getSoundVelocity; // default sound velocity tripods_container tripods; // tripods transmitters_container transmitters; // transmitters hydrophones_container hydrophones; // hydrophones double Q; int debug; try { JParser<> zap("Example program to plot hydrophone data."); zap['f'] = make_field(inputFile); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['o'] = make_field(outputFile) = "hydrophone.root"; zap['a'] = make_field(detectorFile); 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['W'] = make_field(getEmitterID, "waveform identification data") = JPARSER::initialised(); zap['Q'] = make_field(Q, "quality") = 0.0; zap['d'] = make_field(debug) = 2; 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 (hydrophones_container::const_iterator i = hydrophones.begin(); i != hydrophones.end(); ++i) { for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) { if (i->getLocation() == module->getLocation()) { receivers[module->getID()] = JReceiver(module->getID(), module->getPosition() + i->getPosition(), module->getT0() * 1.0e-9); break; } } } 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 T0(-0.05, +0.05); const JRange T1(-0.01, +0.01); TH1D h0(MAKE_CSTRING("Q0 " << T0), NULL, 100, 0.0, 8.0); TH1D h1(MAKE_CSTRING("Q1 " << T1), NULL, 100, 0.0, 8.0); JManager H1(new TH1D("[%]", NULL, 50000, -5.0e-2, +5.0e-2)); while (inputFile.hasNext()) { STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); const JEvent* evt = inputFile.next(); if (!evt->empty() && emitters.has(evt->getID())) { map > buffer; for (JEvent::const_iterator hit = evt->begin(); hit != evt->end(); ++hit) { if (receivers.has(hit->getID())) { buffer[hit->getID()].push_back(*hit); } } for (map >::iterator i = buffer.begin(); i != buffer.end(); ++i) { sort(i->second.begin(), i->second.end(), make_comparator(&JTransmission::getQ, JComparison::gt())); vector::const_iterator hit = i->second.begin(); if (hit->getQ() >= Q) { const JPosition3D& p1 = emitters [evt->getID()].getPosition(); const JPosition3D& p2 = receivers[hit->getID()].getPosition(); const double D = p2.getDistance(p1); const double Vi = V.getInverseVelocity(D, p1.getZ(), p2.getZ()); const double t0 = evt->begin()->getToE() + D * Vi; const double t1 = hit->getToA() - t0; H1[key_type(hit->getID(), evt->getID())]->Fill(t1); const double Q = log10(hit->getQ()); if (T0(t1)) { h0.Fill(Q); } if (T1(t1)) { h1.Fill(Q); } } } } } STATUS(endl); TFile out(outputFile.c_str(), "recreate"); out << h0 << h1 << H1; out.Write(); out.Close(); }