Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JAcousticsTriggerProcessor.cc File Reference

Main program to trigger acoustic data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JHydrophone.hh"
#include "JLang/JClonable.hh"
#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"
#include "JTools/JRange.hh"
#include "JTools/JHashMap.hh"
#include "JROOT/JManager.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"
#include "JAcoustics/JToA.hh"
#include "JAcoustics/JTransmission.hh"
#include "JAcoustics/JReceiver.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JConstantSoundVelocity.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JAcousticsSupportkit.hh"
#include "JAcoustics/JTriggerParameters.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JSupport.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVertex3D.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "Jeep/JContainer.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JTrigger/JMatch.hh"
#include "JTrigger/JAlgorithm.hh"

Go to the source code of this file.

Classes

struct  JACOUSTICS::hit_type
 Acoustic hit. More...
 
class  JACOUSTICS::JMatch3D
 3D match criterion for acoustic signals. More...
 
struct  JACOUSTICS::event_type
 Event. More...
 
struct  JACOUSTICS::JEventOverlap
 Match of two events considering overlap in time and position. More...
 
struct  JACOUSTICS::vertex_type
 Vertex. More...
 
struct  JACOUSTICS::JVelo
 Vertex locator. More...
 

Namespaces

namespace  JACOUSTICS
 Auxiliary classes and methods for acoustic position calibration.
 

Functions

double JACOUSTICS::getQuality (const JEvent &evt)
 Get average quality.
 
int main (int argc, char **argv)
 

Detailed Description

Main program to trigger acoustic data.

An acoustic event is based on coincidences between times of arrival.
If the number of coincident times of arrival exceeds a preset minimum, the event is triggered and subsequently stored in the output file.

Author
mdejong

Definition in file JAcousticsTriggerProcessor.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 452 of file JAcousticsTriggerProcessor.cc.

453{
454 using namespace std;
455 using namespace JPP;
456
458
460
462 JLimit_t& numberOfEvents = inputFile.getLimit();
463 JTriggerParameters parameters;
464 int factoryLimit = 10000;
465 double TMaxExtra_s = 500.0e-6;
466 double RMax_m = 0.0;
467 double DMax_m = 0.0;
468 double Xv_m = 0.0;
469 double factor = 2.0;
470 double Z_m = 0.0; // z-position of emitters
471 double fraction = 0.75; // fraction of working modules
473 JSoundVelocity V = getSoundVelocity; // default sound velocity
474 string detectorFile;
475 hydrophones_container hydrophones; // hydrophones
476 double precision;
477 set<int> ID;
478 int debug;
479
480 try {
481
482 JProperties properties;
483
484 properties.insert(gmake_property(parameters.Q));
485 properties.insert(gmake_property(parameters.TMax_s));
486 properties.insert(gmake_property(parameters.numberOfHits));
487 properties.insert(gmake_property(fraction));
488 properties.insert(gmake_property(factoryLimit));
489 properties.insert(gmake_property(TMaxExtra_s));
490 properties.insert(gmake_property(RMax_m));
491 properties.insert(gmake_property(DMax_m));
492 properties.insert(gmake_property(Xv_m));
493 properties.insert(gmake_property(factor));
494 properties.insert(gmake_property(Z_m));
495
496 JParser<> zap("Main program to trigger acoustic data.");
497
498 zap['f'] = make_field(inputFile, "output of JToA");
499 zap['n'] = make_field(numberOfEvents) = JLimit::max();
500 zap['@'] = make_field(properties, "trigger parameters");
501 zap['o'] = make_field(outputFile, "output file") = "event.root";
502 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
503 zap['a'] = make_field(detectorFile, "detector file");
504 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
505 zap['p'] = make_field(precision, "precision time-of-arrival") = 1.0e-6;
506 zap['W'] = make_field(ID) = JPARSER::initialised();
507 zap['d'] = make_field(debug) = 1;
508
509 zap(argc, argv);
510 }
511 catch(const exception &error) {
512 FATAL(error.what() << endl);
513 }
514
515
517
518 try {
519 load(detectorFile, detector);
520 }
521 catch(const JException& error) {
522 FATAL(error);
523 }
524
525 V.set(detector.getUTMZ());
526
527 const JCylinder3D cylinder(detector.begin(), detector.end());
528
529 const JConstantSoundVelocity V0 = V(cylinder.getZmax());
530
531 const JVelo velo(V0, JVector3D(cylinder.getX(), cylinder.getY(), Z_m), RMax_m, Xv_m, factor);
532
533 const JMatch3D match(V0, TMaxExtra_s);
534 const JEventOverlap overlap2D(parameters.TMax_s, DMax_m);
535 const JEventOverlap overlap1D(parameters.TMax_s);
536
537 JHashMap<int, JReceiver> receivers;
538
539 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
540
541 JPosition3D pos(0.0, 0.0, 0.0);
542
543 if (i->getFloor() == 0) { // get relative position of hydrophone
544
545 try {
546 pos = getPosition(hydrophones.begin(),
547 hydrophones.end(),
548 make_predicate(&JHydrophone::getString, i->getString()));
549 }
550 catch(const exception&) {
551 continue; // if no position available, discard hydrophone
552 }
553 }
554
555 receivers[i->getID()] = JReceiver(i->getID(),
556 i->getPosition() + pos,
557 i->getT0() * 1.0e-9);
558 }
559
560 // monitoring
561
562 JManager<int, TH1D> M1(new TH1D("M1[%]", NULL, 2001, -0.5, 2000.5));
563 JManager<int, TH1D> M2(new TH1D("M2[%]", NULL, 2001, -0.5, 2000.5));
564 JManager<int, TH1D> M3(new TH1D("M3[%]", NULL, 2001, -0.5, 2000.5));
565 JManager<int, TH1D> M4(new TH1D("M4[%]", NULL, 2001, -0.5, 2000.5));
566 JManager<int, TH1D> M5(new TH1D("M5[%]", NULL, 2001, -0.5, 2000.5));
567
568
569 outputFile.open();
570
571 if (!outputFile.is_open()) {
572 FATAL("Error opening file " << outputFile << endl);
573 }
574
575 outputFile.put(JMeta(argc, argv));
576 outputFile.put(parameters);
577
578
579 // input data
580
581 typedef vector<hit_type> buffer_type; // data type
582
583 map<int, map<int, buffer_type> > f1; // emitter -> receiver -> data
584
585 set<int> ids; // working receivers
586
587 static double TOA_s = numeric_limits<double>::max();
588
589 while (inputFile.hasNext()) {
590
591 if (inputFile.getCounter()%100000 == 0) {
592 STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
593 }
594
595 JToA* parameters = inputFile.next();
596
597 if (detector.getID() != parameters->DETID) { // consistency check
598 FATAL("Invalid detector identifier " << parameters->DETID << " != " << detector.getID() << endl);
599 }
600
601 if (receivers.has(parameters->DOMID)) {
602
603 ids.insert(parameters->DOMID);
604
605 const JReceiver& receiver = receivers[parameters->DOMID];
606 const double toa_s = parameters->TOA_S();
607
608 if (toa_s < TOA_s) {
609 TOA_s = toa_s;
610 }
611
612 const JTransmission transmission(parameters->RUN,
613 parameters->DOMID,
614 parameters->QUALITYFACTOR,
615 parameters->QUALITYNORMALISATION,
616 receiver.getT(toa_s),
617 receiver.getT(toa_s));
618
619 if (ID.empty() || ID.count(getWaveformID(parameters->WAVEFORMID)) != 0) {
620 f1[getWaveformID(parameters->WAVEFORMID)][receiver.getID()].push_back(hit_type(receiver.getPosition(), transmission));
621 }
622 }
623 }
624 STATUS(endl);
625
626 if (parameters.numberOfHits == 0) {
627
628 parameters.numberOfHits = (int) (ids.size() * fraction);
629
630 STATUS("Number of hits " << parameters.numberOfHits << endl);
631 }
632
633 vector<event_type> queue;
634
635 for (map<int, map<int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
636
638
639 for (map<int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
640
641 // filter similar hits
642
643 sort(receiver->second.begin(), receiver->second.end(), JTransmission::compare(precision));
644
645 buffer_type::iterator __end = unique(receiver->second.begin(), receiver->second.end(), JTransmission::equals(precision));
646
647 // selection based on quality
648
649 for (buffer_type::const_iterator p = receiver->second.begin(); p != __end; ++p) {
650 if (p->getQ() >= parameters.Q * (parameters.Q <= 1.0 ? p->getW() : 1.0)) {
651 data.push_back(*p);
652 }
653 }
654 }
655
656 if (!data.empty()) {
657
658 sort(data.begin(), data.end(), make_comparator(&hit_type::getToA, JComparison::lt())); // sort according time-of-arrival
659
660 double t0 = TOA_s;
661
662 buffer_type buffer(factoryLimit); // local buffer of hits
663 event_type out[2]; // FIFO of events
664
665 TH1D* h1 = M1[i->first];
666 TH1D* h2 = M2[i->first];
667 TH1D* h3 = M3[i->first];
668 TH1D* h4 = M4[i->first];
669 TH1D* h5 = M5[i->first];
670
671 for (buffer_type::const_iterator p = data.begin(); p != data.end(); ++p) {
672
673 if (distance(data.cbegin(), p)%100 == 0) {
674 STATUS("processed: " << setw(3) << i->first << ' ' << FIXED(9,3) << p->getToA() - data.begin()->getToA() << " [s]" << '\r' << flush); DEBUG(endl);
675 }
676
677 buffer_type::const_iterator q = p;
678
679 while (++q != data.end() && q->getToA() - p->getToA() <= parameters.TMax_s) {}
680
681 h1->Fill(distance(p,q));
682
683 if (distance(p,q) >= parameters.numberOfHits) {
684
685 if (distance(p,q) < factoryLimit) {
686
687 if ((int) out[1].size() >= parameters.numberOfHits && velo(out[1].getVertex(), p, q) >= parameters.numberOfHits) {
688
689 // ?
690
691 } else {
692
693 buffer_type::iterator root = buffer.begin();
694 buffer_type::iterator __p = buffer.begin();
695 buffer_type::iterator __q = buffer.begin();
696
697 *root = *p;
698
699 ++__p;
700 ++__q;
701
702 for (buffer_type::const_iterator i = p; ++i != q; ) {
703 if (match(*p,*i)) {
704 *__q = *i;
705 ++__q;
706 }
707 }
708
709 h2->Fill(distance(root,__q));
710
711 if (distance(root,__q) >= parameters.numberOfHits) {
712
713 __q = clusterize(__p, __q, match);
714
715 h3->Fill(distance(root,__q));
716
717 if (distance(root,__q) >= parameters.numberOfHits) {
718
719 const vertex_type vx = velo(*root, __p, __q, parameters.numberOfHits);
720
721 h4->Fill(vx.N);
722
723 if (vx.N >= parameters.numberOfHits) {
724 out[1] = event_type(JEvent(detector.getID(), out[0].getCounter() + 1, i->first, root, __q), vx.getVertex());
725 }
726 }
727 }
728 }
729
730 } else {
731
732 out[1] = event_type(JEvent(detector.getID(), out[0].getCounter() + 1, i->first, p, q));
733 }
734
735 if (!out[1].empty()) {
736
737 if (out[0].empty()) {
738
739 out[0] = out[1]; // shift
740
741 } else if (overlap2D(out[0],out[1])) {
742
743 out[0].merge(out[1]); // merge
744
745 } else {
746
747 const double t1 = out[0].begin()->getToA();
748
749 STATUS(endl);
750 STATUS("trigger: "
751 << setw(3) << i->first << ' '
752 << FIXED(9,3) << t1 - t0 << " [s]" << ' '
753 << setw(4) << out[0].size() << ' '
754 << FIXED(9,2) << getQuality(out[0]) << ' '
755 << out[0].getPosition() << endl);
756
757 //t0 = t1;
758
759 h5->Fill(out[0].size());
760
761 queue.push_back(out[0]); // store
762
763 out[0] = out[1]; // shift
764 }
765
766 out[1].clear();
767 }
768 }
769 }
770
771 if (!out[0].empty()) {
772
773 queue.push_back(out[0]); // store
774 }
775 STATUS(endl);
776 }
777 }
778
779 if (!queue.empty()) {
780
781 struct time_t {
782
783 void operator=(const double t0) { this->t0 = t0; }
784
785 operator double() const { return this->t0; }
786
787 double t0 = TOA_s;
788 };
789
791
792 sort(queue.begin(), queue.end());
793
794 for (vector<event_type>::iterator p = queue.begin(); p != queue.end(); ) {
795
798
799 // select event with highest quality
800
801 while (++q != queue.end() && overlap1D(*p,*q)) {
802 if (getQuality(*q) > getQuality(*r)) {
803 r = q;
804 }
805 }
806
807 const double t0 = ts[r->getID()];
808 const double t1 = r->begin()->getToA();
809
810 STATUS("event: "
811 << setw(3) << r->getID() << ' '
812 << FIXED(9,3) << t1 - t0 << " [s]" << ' '
813 << setw(4) << r->size() << ' '
814 << FIXED(9,2) << getQuality(*r) << ' '
815 << r->getPosition() << endl);
816
817 //ts[r->getID()] = t1;
818
819 outputFile.put(*r);
820
821 p = q;
822 }
823 }
824
825 JMultipleFileScanner<JMeta> io(inputFile);
826
827 io >> outputFile;
828
829 for (const auto M : { &M1, &M2, &M3, &M4, &M5 }){
830 for (const auto& i : *M) {
831 outputFile.put(*i.second);
832 }
833 }
834
835 outputFile.close();
836}
string outputFile
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
3D match criterion for acoustic signals.
Definition JBillabong.cc:96
Detector data structure.
Definition JDetector.hh:96
int getString() const
Get string number.
Definition JLocation.hh:135
Utility class to parse parameter values.
Data structure for position in three dimensions.
const JPosition3D & getPosition() const
Get position.
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
General exception.
Definition JException.hh:24
int getID() const
Get identifier.
Definition JObjectID.hh:50
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
Object writing to file.
General purpose class for object reading from a list of file names.
virtual bool hasNext() override
Check availability of next element.
counter_type getCounter() const
Get counter.
virtual const pointer_type & next() override
Get next element.
bool has(const T &value) const
Test whether given value is present.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
JVertex3D getVertex(const Trk &track)
Get vertex.
JPosition3D getPosition(const Vec &pos)
Get position.
int getWaveformID(int id)
Get waveform identifier.
double getQuality(const JEvent &evt)
Get average quality.
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition JSydney.cc:80
static const JSoundVelocity getSoundVelocity(1541.0, -17.0e-3, -2000.0)
Function object for velocity of sound.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::vector< JHitW0 > buffer_type
hits
Definition JPerth.cc:70
JFIT::JEvent JEvent
Definition JHistory.hh:404
JHitIterator_t clusterize(JHitIterator_t __begin, JHitIterator_t __end, const JMatch_t &match, const int Nmin=1)
Partition data according given binary match operator.
Definition JAlgorithm.hh:38
Definition root.py:1
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Detector file.
Definition JHead.hh:227
Implementation for depth independend velocity of sound.
Match of two events considering overlap in time and position.
void merge(const JEvent &event)
Merge event.
Acoustic receiver.
Definition JReceiver.hh:30
double getT(const double t_s) const
Get corrected time.
Definition JReceiver.hh:72
Implementation for depth dependend velocity of sound.
JSoundVelocity & set(const double z0)
Set depth.
Time-of-arrival data from acoustic piezo sensor or hydrophone.
Definition JToA.hh:28
uint32_t DOMID
DAQ run number.
Definition JToA.hh:34
uint32_t QUALITYFACTOR
The ticks (16ns) part of the DAQ frame timestamp.
Definition JToA.hh:39
uint32_t QUALITYNORMALISATION
A measure of how good the waveform match was to the signal.
Definition JToA.hh:40
int32_t WAVEFORMID
DOM unique identifeir.
Definition JToA.hh:35
int32_t DETID
Definition JToA.hh:32
int32_t RUN
detector identifier
Definition JToA.hh:33
double TOA_S() const
Time of Arrival, expressed in seconds relative to Unix epoch (1 January 1970 00:00:00 UTC)
Definition JToAImp.cc:40
Auxiliary class to compare transmissions.
Auxiliary class to compare transmissions.
Acoustic transmission.
double Q
minimal quality if larger than one; else minimal normalised quality
double TMax_s
maximal difference between times of emission [s]
int numberOfHits
minimal number of hits to trigger event
Acoustic hit.
Definition JBillabong.cc:70
const JVertex3D & getVertex() const
Get vertex.
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition JContainer.hh:42
Type list.
Definition JTypeList.hh:23
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
General purpose class for hash map of unique keys.
Definition JHashMap.hh:75