Jpp test-rotations-old
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 "TGraph.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 "JROOT/JGraph.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 454 of file JAcousticsTriggerProcessor.cc.

455{
456 using namespace std;
457 using namespace JPP;
458
460
462
464 JLimit_t& numberOfEvents = inputFile.getLimit();
465 JTriggerParameters parameters;
466 int factoryLimit = 10000;
467 double TMaxExtra_s = 500.0e-6;
468 double RMax_m = 0.0;
469 double DMax_m = 0.0;
470 double Xv_m = 0.0;
471 double factor = 2.0;
472 double Z_m = 0.0; // z-position of emitters
473 double fraction = 0.75; // fraction of working modules
475 JSoundVelocity V = getSoundVelocity; // default sound velocity
476 string detectorFile;
477 hydrophones_container hydrophones; // hydrophones
478 double precision;
479 set<int> ID;
480 int debug;
481
482 try {
483
484 JProperties properties;
485
486 properties.insert(gmake_property(parameters.Q));
487 properties.insert(gmake_property(parameters.TMax_s));
488 properties.insert(gmake_property(parameters.numberOfHits));
489 properties.insert(gmake_property(fraction));
490 properties.insert(gmake_property(factoryLimit));
491 properties.insert(gmake_property(TMaxExtra_s));
492 properties.insert(gmake_property(RMax_m));
493 properties.insert(gmake_property(DMax_m));
494 properties.insert(gmake_property(Xv_m));
495 properties.insert(gmake_property(factor));
496 properties.insert(gmake_property(Z_m));
497
498 JParser<> zap("Main program to trigger acoustic data.");
499
500 zap['f'] = make_field(inputFile, "output of JToA");
501 zap['n'] = make_field(numberOfEvents) = JLimit::max();
502 zap['@'] = make_field(properties, "trigger parameters");
503 zap['o'] = make_field(outputFile, "output file") = "event.root";
504 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
505 zap['a'] = make_field(detectorFile, "detector file");
506 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
507 zap['p'] = make_field(precision, "precision time-of-arrival") = 1.0e-6;
508 zap['W'] = make_field(ID) = JPARSER::initialised();
509 zap['d'] = make_field(debug) = 1;
510
511 zap(argc, argv);
512 }
513 catch(const exception &error) {
514 FATAL(error.what() << endl);
515 }
516
517
519
520 try {
521 load(detectorFile, detector);
522 }
523 catch(const JException& error) {
524 FATAL(error);
525 }
526
527 V.set(detector.getUTMZ());
528
529 const JCylinder3D cylinder(detector.begin(), detector.end());
530
531 const JConstantSoundVelocity V0 = V(cylinder.getZmax());
532
533 const JVelo velo(V0, JVector3D(cylinder.getX(), cylinder.getY(), Z_m), RMax_m, Xv_m, factor);
534
535 const JMatch3D match(V0, TMaxExtra_s);
536 const JEventOverlap overlap2D(parameters.TMax_s, DMax_m);
537 const JEventOverlap overlap1D(parameters.TMax_s);
538
539 JHashMap<int, JReceiver> receivers;
540
541 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
542
543 JPosition3D pos(0.0, 0.0, 0.0);
544
545 if (i->getFloor() == 0) { // get relative position of hydrophone
546
547 try {
548 pos = getPosition(hydrophones.begin(),
549 hydrophones.end(),
550 make_predicate(&JHydrophone::getString, i->getString()));
551 }
552 catch(const exception&) {
553 continue; // if no position available, discard hydrophone
554 }
555 }
556
557 receivers[i->getID()] = JReceiver(i->getID(),
558 i->getPosition() + pos,
559 i->getT0() * 1.0e-9);
560 }
561
562 // monitoring
563
564 JManager<int, TH1D> M1(new TH1D("M1[%]", NULL, 2001, -0.5, 2000.5));
565 JManager<int, TH1D> M2(new TH1D("M2[%]", NULL, 2001, -0.5, 2000.5));
566 JManager<int, TH1D> M3(new TH1D("M3[%]", NULL, 2001, -0.5, 2000.5));
567 JManager<int, TH1D> M4(new TH1D("M4[%]", NULL, 2001, -0.5, 2000.5));
568 JManager<int, TH1D> M5(new TH1D("M5[%]", NULL, 2001, -0.5, 2000.5));
569
571
572 outputFile.open();
573
574 if (!outputFile.is_open()) {
575 FATAL("Error opening file " << outputFile << endl);
576 }
577
578 outputFile.put(JMeta(argc, argv));
579 outputFile.put(parameters);
580
581
582 // input data
583
584 typedef vector<hit_type> buffer_type; // data type
585
586 map<int, map<int, buffer_type> > f1; // emitter -> receiver -> data
587
588 set<int> ids; // working receivers
589
590 static double TOA_s = numeric_limits<double>::max();
591
592 while (inputFile.hasNext()) {
593
594 if (inputFile.getCounter()%100000 == 0) {
595 STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
596 }
597
598 JToA* parameters = inputFile.next();
599
600 if (detector.getID() != parameters->DETID) { // consistency check
601 FATAL("Invalid detector identifier " << parameters->DETID << " != " << detector.getID() << endl);
602 }
603
604 if (receivers.has(parameters->DOMID)) {
605
606 ids.insert(parameters->DOMID);
607
608 const JReceiver& receiver = receivers[parameters->DOMID];
609 const double toa_s = parameters->TOA_S();
610
611 if (toa_s < TOA_s) {
612 TOA_s = toa_s;
613 }
614
615 const JTransmission transmission(parameters->RUN,
616 parameters->DOMID,
617 parameters->QUALITYFACTOR,
618 parameters->QUALITYNORMALISATION,
619 receiver.getT(toa_s),
620 receiver.getT(toa_s));
621
622 if (ID.empty() || ID.count(getWaveformID(parameters->WAVEFORMID)) != 0) {
623 f1[getWaveformID(parameters->WAVEFORMID)][receiver.getID()].push_back(hit_type(receiver.getPosition(), transmission));
624 }
625 }
626 }
627 STATUS(endl);
628
629 if (parameters.numberOfHits == 0) {
630
631 parameters.numberOfHits = (int) (ids.size() * fraction);
632
633 STATUS("Number of hits " << parameters.numberOfHits << endl);
634 }
635
636 vector<event_type> queue;
637
638 for (map<int, map<int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
639
641
642 for (map<int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {
643
644 // filter similar hits
645
646 sort(receiver->second.begin(), receiver->second.end(), JTransmission::compare(precision));
647
648 buffer_type::iterator __end = unique(receiver->second.begin(), receiver->second.end(), JTransmission::equals(precision));
649
650 // selection based on quality
651
652 for (buffer_type::const_iterator p = receiver->second.begin(); p != __end; ++p) {
653 if (p->getQ() >= parameters.Q * (parameters.Q <= 1.0 ? p->getW() : 1.0)) {
654 data.push_back(*p);
655 }
656 }
657 }
658
659 if (!data.empty()) {
660
661 sort(data.begin(), data.end(), make_comparator(&hit_type::getToA, JComparison::lt())); // sort according time-of-arrival
662
663 double t0 = TOA_s;
664
665 buffer_type buffer(factoryLimit); // local buffer of hits
666 event_type out[2]; // FIFO of events
667
668 TH1D* h1 = M1[i->first];
669 TH1D* h2 = M2[i->first];
670 TH1D* h3 = M3[i->first];
671 TH1D* h4 = M4[i->first];
672 TH1D* h5 = M5[i->first];
673
674 for (buffer_type::const_iterator p = data.begin(); p != data.end(); ++p) {
675
676 if (distance(data.cbegin(), p)%100 == 0) {
677 STATUS("processed: " << setw(3) << i->first << ' ' << FIXED(9,3) << p->getToA() - data.begin()->getToA() << " [s]" << '\r' << flush); DEBUG(endl);
678 }
679
680 buffer_type::const_iterator q = p;
681
682 while (++q != data.end() && q->getToA() - p->getToA() <= parameters.TMax_s) {}
683
684 h1->Fill(distance(p,q));
685
686 if (distance(p,q) >= parameters.numberOfHits) {
687
688 if (distance(p,q) < factoryLimit) {
689
690 if ((int) out[1].size() >= parameters.numberOfHits && velo(out[1].getVertex(), p, q) >= parameters.numberOfHits) {
691
692 // ?
693
694 } else {
695
696 buffer_type::iterator root = buffer.begin();
697 buffer_type::iterator __p = buffer.begin();
698 buffer_type::iterator __q = buffer.begin();
699
700 *root = *p;
701
702 ++__p;
703 ++__q;
704
705 for (buffer_type::const_iterator i = p; ++i != q; ) {
706 if (match(*p,*i)) {
707 *__q = *i;
708 ++__q;
709 }
710 }
711
712 h2->Fill(distance(root,__q));
713
714 if (distance(root,__q) >= parameters.numberOfHits) {
715
716 __q = clusterize(__p, __q, match);
717
718 h3->Fill(distance(root,__q));
719
720 if (distance(root,__q) >= parameters.numberOfHits) {
721
722 const vertex_type vx = velo(*root, __p, __q, parameters.numberOfHits);
723
724 h4->Fill(vx.N);
725
726 if (vx.N >= parameters.numberOfHits) {
727 out[1] = event_type(JEvent(detector.getID(), out[0].getCounter() + 1, i->first, root, __q), vx.getVertex());
728 }
729 }
730 }
731 }
732
733 } else {
734
735 out[1] = event_type(JEvent(detector.getID(), out[0].getCounter() + 1, i->first, p, q));
736 }
737
738 if (!out[1].empty()) {
739
740 if (out[0].empty()) {
741
742 out[0] = out[1]; // shift
743
744 } else if (overlap2D(out[0],out[1])) {
745
746 out[0].merge(out[1]); // merge
747
748 } else {
749
750 const double t1 = out[0].begin()->getToA();
751
752 STATUS(endl);
753 STATUS("trigger: "
754 << setw(3) << i->first << ' '
755 << FIXED(9,3) << t1 - t0 << " [s]" << ' '
756 << setw(4) << out[0].size() << ' '
757 << FIXED(9,2) << getQuality(out[0]) << ' '
758 << out[0].getPosition() << endl);
759
760 //t0 = t1;
761
762 h5->Fill(out[0].size());
763
764 queue.push_back(out[0]); // store
765
766 out[0] = out[1]; // shift
767 }
768
769 out[1].clear();
770 }
771 }
772 }
773
774 if (!out[0].empty()) {
775
776 queue.push_back(out[0]); // store
777 }
778 STATUS(endl);
779 }
780 }
781
782 if (!queue.empty()) {
783
784 struct time_t {
785
786 void operator=(const double t0) { this->t0 = t0; }
787
788 operator double() const { return this->t0; }
789
790 double t0 = TOA_s;
791 };
792
794
795 sort(queue.begin(), queue.end());
796
797 for (vector<event_type>::iterator p = queue.begin(); p != queue.end(); ) {
798
801
802 // select event with highest quality
803
804 while (++q != queue.end() && overlap1D(*p,*q)) {
805 if (getQuality(*q) > getQuality(*r)) {
806 r = q;
807 }
808 }
809
810 const double t0 = ts[r->getID()];
811 const double t1 = r->begin()->getToA();
812
813 STATUS("event: "
814 << setw(3) << r->getID() << ' '
815 << FIXED(9,3) << t1 - t0 << " [s]" << ' '
816 << setw(4) << r->size() << ' '
817 << FIXED(9,2) << getQuality(*r) << ' '
818 << r->getPosition() << endl);
819
820 //ts[r->getID()] = t1;
821
822 outputFile.put(*r);
823
824 p = q;
825
826 G2[r->getID()].put(r->getX(), r->getY());
827 }
828 }
829
830 JMultipleFileScanner<JMeta> io(inputFile);
831
832 io >> outputFile;
833
834 for (const auto M : { &M1, &M2, &M3, &M4, &M5 }){
835 for (const auto& i : *M) {
836 outputFile.put(*i.second);
837 }
838 }
839
840 for (map<int, JGraph_t>::const_iterator i = G2.begin(); i != G2.end(); ++i) {
841 outputFile.put(JGraph(i->second, MAKE_CSTRING("G[" << FILL(2,'0') << i->first << "]")));
842 }
843
844 outputFile.close();
845}
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 MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
#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 sequence of same character.
Definition JManip.hh:330
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 data structure to build TGraph.
Definition JGraph.hh:44
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