10 #include <type_traits>
68 namespace JACOUSTICS {
74 using namespace JDETECTOR;
119 bool hasString(
const int id)
const
135 static constexpr
double RADIUS_M = 1.0;
156 std::
set<int>(buffer.begin(), buffer.end())
170 std::set_difference(A.begin(), A.end(), B.begin(), B.end(), std::inserter(*
this, this->begin()));
186 std::set_difference(A.begin(), A.end(), B.begin(), B.end(), std::inserter(*
this, this->begin()));
199 for (
int id; in >> id; ) {
220 for (
const int id :
object) {
283 virtual void apply(
const double step)
override
328 virtual void apply(
const double step)
override
330 for (
const auto i : index) {
372 virtual void apply(
const double step)
override
374 for (
const auto i : index) {
421 virtual void apply(
const double step)
override
423 tripods[index].add(direction * step);
446 hydrophones (setup.hydrophones),
447 transmitters(setup.transmitters)
462 virtual void apply(
const double step)
override
468 if (index[0] != hydrophones .size()) { hydrophones [index[0]].rotate(R); }
469 if (index[1] != transmitters.size()) { transmitters[index[1]].rotate(R); }
489 Nmax (std::numeric_limits<size_t>::max()),
515 for (
double value; in >> value; ) {
516 object.steps.push_back(value);
539 out << setw(2) <<
object.option <<
' '
540 << setw(2) <<
object.mestimator <<
' '
543 << setw(3) <<
object.Nextra;
545 for (
const double value :
object.steps) {
546 out <<
' ' <<
FIXED(9,5) << value;
570 const size_t threads,
572 filenames(filenames),
579 setup.tripods.load(filenames.
tripod.c_str());
584 for (JDetector::const_iterator
i = setup.detector.begin();
i != setup.detector.end(); ++
i) {
585 receivers[
i->getID()] =
i->getLocation();
592 for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) {
598 fits.initialise(setup);
600 this->V.
set(setup.detector.getUTMZ());
607 ROOT::EnableThreadSafety();
647 (*this)[module.
getString()].push_back(module);
668 for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) {
670 if (fits.strings .count(module->getString()) == 0)
671 A.push_back(*module);
672 else if (fits.transmitters.count(module->getString()) != 0)
673 B.push_back(*module);
675 C.push_back(*module);
678 setup.detector.swap(
A);
680 for (
const auto& element : B) {
681 setup.detector.push_back(element.second.base);
688 for (
const int i : fits.tripods) {
694 for (
const int i : fits.transmitters) {
700 const double chi2 = fit(*
this);
702 for (
auto& element : B) {
704 const JModule& base = setup.detector.getModule(router->getAddress(
JLocation(element.second.base.getString(),0)));
707 for (string_type::iterator module = element.second.begin(); module != element.second.end(); ++module) {
711 setup.detector.push_back(*module);
715 copy(
C.begin(),
C.end(), back_inserter(setup.detector));
739 for (
const int i : fits.strings) {
745 for (
const int i : fits.tripods) {
751 for (
const int i :
ids_t(fits.hydrophones, fits.transmitters)) {
755 for (
const int i : fits.transmitters) {
757 JModule& module = setup.detector.getModule(router->getAddress(
JLocation(i,0)));
763 const double chi2 = fit(*
this);
785 for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) {
787 buffer[module->getString()].push_back(*module);
789 if (module->getZ() > z0) {
796 for (transmitters_container::iterator
i = setup.transmitters.begin();
i != setup.transmitters.end(); ++
i) {
798 tx.push_back(router->getModule(
i->getLocation()));
800 catch(
const exception&) {}
803 for (
const int i : fits.strings) {
805 setup.detector.swap(buffer[i]);
807 copy(tx.begin(), tx.end(), back_inserter(setup.detector));
816 const double chi2 = fit(*
this);
818 STATUS(
"string: " << setw(4) << i <<
' ' <<
FIXED(9,4) << chi2 << endl);
825 setup.detector.clear();
827 for (
const auto& element : buffer) {
828 copy(element.second.begin(), element.second.end(), back_inserter(setup.detector));
849 for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) {
850 if (fits.strings.count(module->getString()) != 0 && module->getFloor() != 0) {
855 const double chi2 = fit(*
this);
872 const JGeometry geometry(setup.detector, setup.hydrophones);
876 for (tripods_container::const_iterator
i = setup.tripods.begin();
i != setup.tripods.end(); ++
i) {
878 emitters[
i->getID()] =
JEmitter(
i->getID(),
i->getUTMPosition() - setup.detector.getUTMPosition());
882 for (transmitters_container::const_iterator
i = setup.transmitters.begin();
i != setup.transmitters.end(); ++
i) {
884 emitters[
i->getID()] =
JEmitter(
i->getID(),
i->getPosition() + router->getModule(
i->getLocation()).
getPosition());
886 catch(
const exception&) {}
892 this->output.clear();
908 for (input_type::const_iterator evt = superevt.begin(); evt != superevt.end(); ++evt) {
910 if (emitters.has(evt->getID())) {
912 const JEmitter& emitter = emitters [evt->getID()];
913 const double weight =
getWeight(evt->getID());
915 for (JEvent::const_iterator
i = evt->begin();
i != evt->end(); ++
i) {
919 data.push_back(
JHit(emitter,
921 receivers[
i->getID()],
926 counter.insert(evt->getID());
940 }
else if (option == 2) {
952 return numeric_limits<float>::max();
967 ifstream
in(script.c_str());
975 if (buffer.empty() || buffer[0] ==
skip_t) {
979 istringstream
is(buffer);
985 fits.initialise(setup);
987 }
else if (key ==
fix_t) {
992 if (is >> type >>
id) {
994 fits.strings .
fix(
id);
995 fits.hydrophones .fix(
id);
996 fits.transmitters.fix(
id);
998 fits.tripods .fix(
id);
1009 if (is >> stage >> input) {
1011 STATUS(
"stage " << setw(3) << stage <<
" {" << input <<
"}" << endl);
1017 ofstream out(
MAKE_CSTRING(
"stage-" << stage <<
".log"));
1022 switch (stage[stage.size() - 1]) {
1076 using namespace JPP;
1078 if (
getFileStatus(dir.c_str()) || (mkdir(dir.c_str(), S_IRWXU | S_IRWXG) != -1)) {
1082 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
1083 module->set(router->getModule(module->getLocation()).
getPosition());
1088 setup.tripods.store(
getFilename(dir, filenames.tripod).c_str());
1090 if (filenames.hydrophone !=
"") { setup.hydrophones .store(
getFilename(dir, filenames.hydrophone) .c_str()); }
1091 if (filenames.transmitter !=
"") { setup.transmitters.store(
getFilename(dir, filenames.transmitter).c_str()); }
1109 for (JDetector::const_iterator
i = setup.detector.begin();
i != setup.detector.end(); ++
i) {
1111 buffer.insert(
i->getID());
1126 using namespace std;
1130 for (tripods_container::const_iterator
i = setup.tripods.begin();
i != setup.tripods.end(); ++
i) {
1131 buffer.insert(
i->getID());
1134 for (transmitters_container::const_iterator
i = setup.transmitters.begin();
i != setup.transmitters.end(); ++
i) {
1138 const JModule& module = router->getModule(
i->getLocation());
1141 buffer.insert(
i->getID());
1144 catch(
const exception&) {}
1181 using namespace std;
1182 using namespace JPP;
1192 disable_container disable;
1198 JParser<> zap(
"Application to perform acoustic pre-calibration.");
1200 zap[
'f'] =
make_field(inputFile,
"output of JAcousticEventBuilder[.sh]");
1201 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
1206 zap[
'T'] =
make_field(filenames.tripod,
"tripod file");
1211 zap[
'N'] =
make_field(threads,
"number of threads") = 1;
1216 catch(
const exception &error) {
1217 FATAL(error.what() << endl);
1221 FATAL(
"Invalid number of threads " << threads << endl);
1233 while (inputFile.hasNext()) {
1235 const JEvent* evt = inputFile.next();
1237 if (emitters.count(evt->
getID())) {
1238 zbuf.push_back(*evt);
1244 sort(zbuf.begin(), zbuf.end());
1246 for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
1248 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() +
parameters.Tmax_s; ) {}
1254 for (buffer_type::iterator evt = p; evt != q; ++evt) {
1260 for (JEvent::iterator
i = evt->begin();
i != __end; ) {
1265 if (receivers.count(
i->getID()) &&
i->getQ() >=
parameters.Qmin * (unit(
parameters.Qmin) ?
i->getW() : 1.0)) {
1270 iter_swap(
i, --__end);
1273 buffer.push_back(
JEvent(evt->getOID(), buffer.size(), evt->getID(), evt->begin(), __end));
1277 sydney.
input.push_back(buffer);
void initialise(const JSetup &setup)
Initialise.
ids_t(const std::vector< int > &buffer)
Copy constructor.
Auxiliary class to edit length of Dyneema ropes.
Utility class to parse command line options.
std::vector< JHydrophone > & hydrophones
JDetector detector
detector
friend std::ostream & operator<<(std::ostream &out, const JParameters_t &object)
Write parameters to output stream.
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
int main(int argc, char *argv[])
void stage_b(const JParameters_t ¶meters)
Fit procedure to determine the stretching and z-positions of individual strings.
Auxiliary data structure for decomposed string.
void stage_a(const JParameters_t ¶meters)
Fit procedure to determine the positions of the strings and tripods.
friend std::istream & operator>>(std::istream &in, ids_t &object)
Read identifiers from input stream.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
double getWeight(T __begin, T __end)
Get total weight of data points.
static const std::string fix_t
fix objects
int getFloor() const
Get floor number.
Data structure for a composite optical module.
JContainer< std::vector< JTransmitter > > transmitters_container
double operator()(const int option) const
Get chi2.
static JDetectorMechanics getMechanics
Function object to get string mechanics.
std::string detector
detector
virtual void apply(const double step) override
Apply step.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
JVector3D getPosition(T __begin, T __end, const JPredicate< JTypename_t, JComparator_t > &predicate)
Get position from element in data which corresponds to given predicate.
Auxiliary data structure to unify weights of acoustics data according to the number of pings per emit...
General purpose class for hash map of unique elements.
ROOT TTree parameter settings.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
virtual void apply(const double step) override
Apply step.
void push_back(const JModule &module)
Add module.
Auxiliary data structure for setup of complete system.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
#define MAKE_CSTRING(A)
Make C-string.
void run(const std::string &script)
Run.
tripods_container tripods
tripods
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
void enqueue(input_type &data)
Queue data.
then usage $script< detector specific pre-calibration script >< option > nAuxiliary script to make scan of pre stretching of detector strings(see JEditDetector)." "\nPossible options
Auxiliary data structure for floating point format specification.
Auxiliary data structure for handling of file names.
Extended data structure for parameters of stage.
transmitters_container transmitters
transmitters
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
ids_t(const ids_t &A, const ids_t &B)
Difference constructor.
Data structure for detector geometry and calibration.
JMODEL::JString getString(const JFit &fit)
Get model parameters of string.
std::vector< JHitW0 > buffer_type
hits
static const JVector3D JVector3X_t(1, 0, 0)
unit x-vector
Data structure for hydrophone.
static const std::string stage_t
fit stage
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
Auxiliary class to edit (x,y,z) position of string.
#define MAKE_STRING(A)
Make string.
ids_t strings
identifiers of strings
Direct access to location in detector data structure.
static const std::string tripod_t
tripod
Implementation of object output from STD container.
static const double C
Physics constants.
List of object identifiers.
std::vector< double > steps
Auxiliary class for defining the range of iterations of objects.
static const JSoundVelocity getSoundVelocity(1541.0,-17.0e-3,-2000.0)
Function object for velocity of sound.
std::string hydrophone
hydrophone
void store(const std::string &dir=".")
Store data in given directory.
JParameters_t()
Default constuctor.
JModuleEditor(JModule &module)
Constructor.
Auxiliary class to edit (z) position of module.
This class can be used to temporarily redirect one output (input) stream to another output (input) st...
JContainer< std::vector< JHydrophone > > hydrophones_container
void fix(const ids_t &B)
Fix.
Data structure for vector in three dimensions.
Router for direct addressing of location data in detector data structure.
Data structure for transmitter.
Logical location of module.
bool has(const int bit) const
Test PMT status.
const JLocation & getLocation() const
Get location.
Auxiliary wrapper for I/O of container with optional comment (see JComment).
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
ids_t()
Default constructor.
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
std::vector< input_type > input
int getID() const
Get identifier.
static JQuantile Q
chi2/NDF
ids_t transmitters
identifiers of strings with transmitter
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Auxiliary class for CPU timing and usage.
std::vector< JEvent > input_type
std::istream & getline(std::istream &in, JString &object)
Read string from input stream until end of line.
Auxiliary data structure for detector with decomposed strings.
static const int TRANSMITTER_DISABLE
Enable (disable) use of transmitter if this status bit is 0 (1);.
Auxiliary data structure for fit parameter.
double getY() const
Get y position.
std::vector< JSuperEvt > output
Auxiliary class to edit orientation of anchor.
const JPosition3D & getPosition() const
Get position.
JContainer< std::vector< JTripod > > tripods_container
General purpose messaging.
virtual void apply(const double step) override
Apply step.
JDyneemaEditor(JSetup &setup, const int id, const double z0=0.0)
Constructor.
Implementation for depth dependend velocity of sound.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
static JStat getFileStatus
Function object for file status.
ids_t getEmitters() const
Get list of identifiers of emitters.
ids_t getReceivers() const
Get list of identifiers of receivers.
Auxiliary data structure for editable parameter.
fits_t()
Default constructor.
void push_back(const JModule &module)
Add module.
virtual void apply(const double step) override
Apply step.
static JQuantile Q
chi2/NDF
Thread pool for global fits using super events.
then JCookie sh JDataQuality D $DETECTOR_ID R
static void overlap(T p, T q, const double Tmax_s)
Empty overlapping events.
Implementation of object iteration from STD container.
int getString() const
Get string number.
std::vector< JTransmitter > & transmitters
friend std::ostream & operator<<(std::ostream &out, const ids_t &object)
Write identifiers to output stream.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
detector()
Default constructor.
void stage_c(const JParameters_t ¶meters)
Fit procedure to determine the z-positions of the modules.
Auxiliary class to define a range between two values.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
Main class for pre-calibration using acoustics data.
static output_type * output
optional output
JFitParameters parameters
Acoustic transmission identifier.
JTOOLS::JHashMap< int, JLocation > receivers
Fit functions of acoustic model.
static const int PIEZO_DISABLE
Enable (disable) use of piezo if this status bit is 0 (1);.
ids_t hydrophones
identifiers of strings with hydrophone
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
JAnchorEditor(JSetup &setup, const int id)
Constructor.
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Thread pool for global fits.
int getID() const
Get identifier.
ids_t tripods
identifiers of tripods
double getX() const
Get x position.
void copy(const Head &from, JHead &to)
Copy header from from to to.
bool hasLocation(const JLocation &location) const
Check if this detector has given location.
static const char skip_t
Script commands.
std::string getFilename(const std::string &file_name)
Get file name part, i.e. part after last JEEP::PATHNAME_SEPARATOR if any.
std::vector< size_t > index
JTripodEditor(JSetup &setup, const int id, const JVector3D &direction)
Constructor.
static const int HYDROPHONE_DISABLE
Enable (disable) use of hydrophone if this status bit is 0 (1);.
std::vector< size_t > index
JStringEditor(JSetup &setup, const int id, const JVector3D &direction, const bool option)
Constructor.
const JLimit & getLimit() const
Get limit.
static const std::string initialise_t
initialise
Exception for accessing a value in a collection that is outside of its range.
long double getMean() const
Get mean value.
Auxiliary data structure for group of lists of identifiers of to-be-fitted objects.
virtual void apply(const double step) override
Apply step.
int getID() const
Get emitter identifier.
static const JVector3D JVector3Y_t(0, 1, 0)
unit y-vector
std::string transmitter
transmitter
Data structure for tripod.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
Acoustic transmission identifier.
JSoundVelocity & set(const double z0)
Set depth.
Auxiliary data structure for floating point format specification.
JModule & set(const JVector3D &pos)
Set position.
friend std::istream & operator>>(std::istream &in, JParameters_t &object)
Read parameters from input stream.
double getZ() const
Get z position.
std::unique_ptr< JLocationRouter > router
JVector3D & add(const JVector3D &vector)
Add vector.
void stage_0(const JParameters_t ¶meters)
Fit procedure to determine the positions of tripods and transmitters using strings that are fixed...
Template definition of fit function of acoustic model.
std::vector< JTripod > & tripods
unsigned long long usec_wall
JSydney(const JFilenames &filenames, const JSoundVelocity &V, const size_t threads, const int debug)
Constructor.
Data structure for optical module.
static const std::string string_t
string
Auxiliary class to edit (x,y,z) position of tripod.