Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JFremantle.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <vector>
5#include <set>
6#include <map>
7#include <deque>
8#include <algorithm>
9#include <limits>
10
11#include <type_traits>
12#include <functional>
13#include <future>
14#include <mutex>
15#include <thread>
16#include <vector>
17#include <queue>
18
19#include "TROOT.h"
20#include "TFile.h"
21
22#include "JLang/JPredicate.hh"
23#include "JLang/JComparator.hh"
24#include "JLang/JComparison.hh"
25#include "JLang/JFileStream.hh"
26
29#include "JDetector/JTripod.hh"
31#include "JDetector/JModule.hh"
33
34#include "JTools/JHashMap.hh"
35#include "JTools/JRange.hh"
36#include "JTools/JQuantile.hh"
37
41#include "JSupport/JMeta.hh"
42
46#include "JAcoustics/JHit.hh"
49#include "JAcoustics/JEvent.hh"
55
56#include "Jeep/JContainer.hh"
57#include "Jeep/JParser.hh"
58#include "Jeep/JMessage.hh"
59
60
61/**
62 * \file
63 *
64 * Application to make a global fit of the detector geometry to acoustic data.\n
65 * \author mdejong
66 */
67int main(int argc, char **argv)
68{
69 using namespace std;
70 using namespace JPP;
71
75 typedef JContainer< set<JTransmission_t> > disable_container;
76
79 string detectorFile;
80 JLimit_t& numberOfEvents = inputFile.getLimit();
81 JSoundVelocity V = getSoundVelocity; // default sound velocity
82 tripods_container tripods; // tripods
83 transmitters_container transmitters; // transmitters
84 hydrophones_container hydrophones; // hydrophones
85 JFitParameters parameters; // fit parameters
86 disable_container disable; // disable tansmissions
87 size_t threads; // number of parallel threads
88 bool& squash = JFremantle::squash;
89 int debug;
90
91 try {
92
93 JParser<> zap("Application to fit position calibration model to acoustic data.");
94
95 zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
96 zap['o'] = make_field(outputFile) = "";
97 zap['n'] = make_field(numberOfEvents) = JLimit::max();
98 zap['a'] = make_field(detectorFile);
99 zap['@'] = make_field(parameters) = JPARSER::initialised();
100 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
101 zap['T'] = make_field(tripods, "tripod data");
102 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
103 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
104 zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
105 zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
106 zap['N'] = make_field(threads, "number of threads") = 1;
107 zap['s'] = make_field(squash, "squash transmissions in output");
108 zap['d'] = make_field(debug) = 1;
109
110 zap(argc, argv);
111 }
112 catch(const exception &error) {
113 FATAL(error.what() << endl);
114 }
115
116 ROOT::EnableThreadSafety();
117
119
120 try {
121 load(detectorFile, detector);
122 }
123 catch(const JException& error) {
124 FATAL(error);
125 }
126
127 JHashMap<int, JLocation> receivers;
129
130 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
131 receivers[i->getID()] = i->getLocation();
132 }
133
134 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
135 emitters[i->getID()] = JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition());
136 }
137
138 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
139 try {
140 emitters[i->getID()] = JEmitter(i->getID(), i->getPosition() + detector.getModule(i->getLocation()).getPosition());
141 }
142 catch(const exception&) {} // if no module available, discard transmitter
143 }
144
145 V.set(detector.getUTMZ()); // sound velocity at detector depth
146
147 JGeometry geometry(detector, hydrophones);
148 JGlobalfit katoomba(geometry, V, parameters);
149
151
154
155 if (inputFile.size() > 1u) { // sort input files
156
158
159 for (const string& file_name : inputFile) {
160
161 STATUS(file_name << '\r'); DEBUG(endl);
162
163 for (JSingleFileScanner<JEvent> in(file_name, 1); in.hasNext(); ) {
164
165 const JEvent* evt = in.next();
166
167 if (JFremantle::detid != evt->getDetectorID()) {
168 FATAL("Invalid detector identifier " << evt->getDetectorID() << " != " << JFremantle::detid << endl);
169 }
170
171 if (!evt->empty()) {
172 zmap[evt->begin()->getToE()] = file_name;
173 }
174 }
175 }
176 STATUS(endl);
177
178 inputFile.clear();
179
180 for (map<double, string>::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
181 inputFile.push_back(i->second);
182 }
183 }
184
185 if (outputFile.getFilename() != "" &&
186 outputFile.getFilename() != "--") {
187
188 outputFile.open();
189
191 }
192 /*
193 outputFile.put(JMeta(argc, argv));
194 outputFile.put(parameters);
195 */
196
197 try {
198
199 JFremantle fremantle(geometry, V, parameters, threads, 2 * threads);
200
201 typedef deque<JEvent> buffer_type;
202
203 for (buffer_type zbuf; inputFile.hasNext(); ) {
204
205 STATUS(inputFile.getFilename() << '\r'); DEBUG(endl);
206
207 // read one file at a time
208
209 for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
210
211 const JEvent* evt = inputFile.next();
212
213 if (!evt->empty() && emitters.has(evt->getID())) {
214 zbuf.push_back(*evt);
215 }
216 }
217
218 sort(zbuf.begin(), zbuf.end()); // sort according first time-of-emission
219
220 for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
221
222 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}
223
224 if (q == zbuf.end()) {
225
226 if (inputFile.hasNext()) {
227
228 zbuf.erase(zbuf.begin(), p); // remove processed data and continue reading
229
230 break;
231 }
232 }
233
234 JEvent::overlap(p, q, parameters.deadTime_s); // empty overlapping events
235
236 if (getNumberOfEmitters(p,q) >= parameters.Nmin) {
237
238 const JWeight getWeight(p, q);
239
241
242 for (buffer_type::iterator evt = p; evt != q; ++evt) {
243
244 sort(evt->begin(), evt->end(), JKatoomba<>::compare);
245
246 JEvent::iterator __end = unique(evt->begin(), evt->end(), make_comparator(&JTransmission::getID, JComparison::eq()));
247
248 const JEmitter& emitter = emitters [evt->getID()];
249 const double weight = getWeight(evt->getID());
250
251 for (JEvent::const_iterator i = evt->begin(); i != __end; ++i) {
252
253 if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 &&
254 disable.count(JTransmission_t(-1, i->getID())) == 0) {
255
256 if (receivers.has(i->getID()) && geometry.hasLocation(receivers[i->getID()]) && i->getQ() >= parameters.Qmin * (parameters.Qmin <= 1.0 ? i->getW() : 1.0)) {
257
258 data.push_back(JHit(emitter,
259 distance(p,evt),
260 receivers[i->getID()],
261 i->getToA(),
262 parameters.sigma_s,
263 weight));
264 }
265 }
266 }
267 }
268
269 if (getMinimumNumberOfEmitters(data.begin(), data.end()) >= parameters.Nmin) {
270 fremantle.enqueue(data);
271 }
272 }
273 }
274 }
275 STATUS(endl);
276 }
277 catch(const exception& error) {
278 FATAL("main " << error.what());
279 }
280 /*
281 JMultipleFileScanner<JMeta> io(inputFile);
282
283 io >> outputFile;
284 */
285 outputFile.close();
286
287 JFileOutputStream(3) << SCIENTIFIC(1,10) << JFremantle::Q.getMean(numeric_limits<float>::max()) << endl;
288}
Acoustics toolkit.
Acoustic event.
Acoustic hit.
ROOT TTree parameter settings.
Container I/O.
string outputFile
Data structure for detector geometry and calibration.
Acoustic emitter.
Recording of objects on file according a format that follows from the file name extension.
Acoustic fit parameters.
int main(int argc, char **argv)
Definition JFremantle.cc:67
Fit function of acoustic model.
General purpose class for hash map of unique elements.
Data structure for hydrophone.
General purpose messaging.
#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
ROOT I/O of application specific meta data.
Data structure for optical module.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Auxiliary class to define a range between two values.
Scanning of objects from a single file according a format that follows from the extension of each fil...
Sound velocity.
Acoustic super event fit toolkit.
Acoustic event fit.
Acoustic transmission identifier.
Data structure for transmitter.
Data structure for tripod.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Thread pool for global fits.
static output_type * output
optional output
void enqueue(input_type &data)
Queue data.
static JMATH::JQuantile_t Q
chi2/NDF
static int detid
detector identifier
static bool squash
squash transmissions in output
Detector data structure.
Definition JDetector.hh:96
General exception.
Definition JException.hh:24
Streaming of output.
Utility class to parse command line options.
Definition JParser.hh:1698
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.
virtual const pointer_type & next() override
Get next element.
const std::string & getFilename() const
Get current file name.
Object reading from a list of files.
virtual bool hasNext() override
Check availability of next element.
bool has(const T &value) const
Test whether given value is present.
size_t getMinimumNumberOfEmitters(T __begin, T __end)
Get minimum number of emitters for any string in data.
double getWeight(T __begin, T __end)
Get total weight of data points.
JContainer< std::vector< JTripod > > tripods_container
Definition JSydney.cc:79
JContainer< std::vector< JTransmitter > > transmitters_container
Definition JSydney.cc:81
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition JSydney.cc:80
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
static const JSoundVelocity getSoundVelocity(1541.0, -17.0e-3, -2000.0)
Function object for velocity of sound.
static JDetectorMechanics getMechanics
Function object to get string mechanics.
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
Detector file.
Definition JHead.hh:227
Acoustic emitter.
Definition JEmitter.hh:30
int getID() const
Get emitter identifier.
const int getDetectorID() const
Get detector identifier.
double Qmin
minimal quality transmission
double deadTime_s
dead time between events [s]
size_t Nmin
minimum number of emitters
double sigma_s
time-of-arrival resolution [s]
double Tmax_s
time window to combine events [s]
bool hasLocation(const JLocation &location) const
Check if this detector has given location.
Definition JGeometry.hh:588
Global fit of prameterised detector geometry to acoustics data.
Definition JGlobalfit.hh:29
Acoustics hit.
Template definition of fit function of acoustic model.
Implementation for depth dependend velocity of sound.
JSoundVelocity & set(const double z0)
Set depth.
Acoustic transmission identifier.
int getID() const
Get identifier.
Auxiliary data structure to unify weights of acoustics data according to the number of pings per emit...
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition JContainer.hh:42
double getMean() const
Get mean value.
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
General purpose class for hash map of unique keys.
Definition JHashMap.hh:75
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488