Jpp 21.0.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JKatoomba.cc File Reference

Application to make a global fit of the detector geometry to acoustic data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <deque>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"
#include "JLang/JComparison.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JTransmitter.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JHydrophone.hh"
#include "JTools/JHashMap.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JGlobalfit.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JTransmission_t.hh"
#include "Jeep/JContainer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Application to make a global fit of the detector geometry to acoustic data.


Author
mdejong

Definition in file software/JAcoustics/JKatoomba.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 60 of file software/JAcoustics/JKatoomba.cc.

61{
62 using namespace std;
63 using namespace JPP;
64
65 typedef JContainer< vector<JTripod> > tripods_container;
66 typedef JContainer< vector<JTransmitter> > transmitters_container;
67 typedef JContainer< vector<JHydrophone> > hydrophones_container;
68 typedef JContainer< set<JTransmission_t> > disable_container;
69
70 JMultipleFileScanner<JEvent> inputFile;
71 JFileRecorder<JTYPELIST<JEvt, JEvent, JFitParameters, JDetectorMechanics, JMeta, TH1D, TH2D>::typelist> outputFile;
72 string detectorFile;
73 JLimit_t& numberOfEvents = inputFile.getLimit();
74 JSoundVelocity V = getSoundVelocity; // default sound velocity
75 tripods_container tripods; // tripods
76 transmitters_container transmitters; // transmitters
77 hydrophones_container hydrophones; // hydrophones
78 JDetectorMechanics mechanics; // mechanical model data
79 JFitParameters parameters; // fit parameters
80 disable_container disable; // disable tansmissions
81 JROOTClassSelection selection = getROOTClassSelection<JTYPELIST<JEvent>::typelist>();
82 JRange<double> utc;
83 bool squash;
84 Long64_t autoflush;
85 int debug;
86
87 try {
88
89 JParser<> zap("Application to fit position calibration model to acoustic data.");
90
91 zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
92 zap['o'] = make_field(outputFile);
93 zap['n'] = make_field(numberOfEvents) = JLimit::max();
94 zap['a'] = make_field(detectorFile);
95 zap['@'] = make_field(parameters) = JPARSER::initialised();
96 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
97 zap['T'] = make_field(tripods, "tripod data");
98 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
99 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
100 zap['M'] = make_field(mechanics, "mechanics data") = JPARSER::initialised();
101 zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
102 zap['C'] = make_field(selection,
103 "Precede name of data structure by a '+' or '-' "
104 "to add or remove data types in the output, respectively."
105 "\nROOT wildcards are accepted.") = JPARSER::initialised();
106 zap['u'] = make_field(utc, "select UTC time range [s]") = JRange<double>();
107 zap['q'] = make_field(squash, "squash meta data");
108 zap['R'] = make_field(autoflush) = 5000;
109 zap['d'] = make_field(debug) = 1;
110
111 zap(argc, argv);
112 }
113 catch(const exception &error) {
114 FATAL(error.what() << endl);
115 }
116
117
118 JDetector detector;
119
120 try {
121 load(detectorFile, detector);
122 }
123 catch(const JException& error) {
124 FATAL(error);
125 }
126
127 const floor_range range = getRangeOfFloors(detector);
128
129 JHashMap<int, JLocation> receivers;
130 JHashMap<int, JEmitter> emitters;
131
132 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
133 receivers[i->getID()] = i->getLocation();
134 }
135
136 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
137 emitters[i->getID()] = JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition());
138 }
139
140 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
141 try {
142 emitters[i->getID()] = JEmitter(i->getID(), i->getPosition() + detector.getModule(i->getLocation()).getPosition());
143 }
144 catch(const exception&) {} // if no module available, discard transmitter
145 }
146
147 V.set(detector.getUTMZ()); // sound velocity at detector depth
148
149 JGeometry geometry(detector, mechanics, hydrophones);
150 JGlobalfit katoomba(geometry, V, parameters);
151
152 JKatoomba<JGandalf>::debug = debug;
153 JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;
154
155 DEBUG(geometry);
156
157 typedef vector<JHit> data_type;
158
159 TH1D h0("chi2/NDF", NULL, 50000, 0.0, 1000.0);
160 TH1D h1("h1", NULL, 51, -0.5, 50.5);
161 TH1D hn("hn", NULL, 100, 0.0, 6.0);
162
163 JManager<int, TH2D> HA(new TH2D("ha[%]", NULL,
164 getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
165 range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
166
167 JManager<int, TH2D> HB(new TH2D("hb[%]", NULL,
168 getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
169 range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
170
171 for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
172 HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
173 HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
174 }
175
176 if (inputFile.size() > 1u) { // sort input files
177
178 map<double, string> zmap;
179
180 for (const string& file_name : inputFile) {
181
182 STATUS(file_name << '\r'); DEBUG(endl);
183
184 for (JSingleFileScanner<JEvent> in(file_name, 1); in.hasNext(); ) {
185
186 const JEvent* evt = in.next();
187
188 if (detector.getID() != evt->getDetectorID()) {
189 FATAL("Invalid detector identifier " << evt->getDetectorID() << " != " << evt->getDetectorID() << endl);
190 }
191
192 if (!evt->empty()) {
193 zmap[evt->begin()->getToE()] = file_name;
194 }
195 }
196 }
197 STATUS(endl);
198
199 inputFile.clear();
200
201 for (map<double, string>::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
202 inputFile.push_back(i->second);
203 }
204 }
205
206 outputFile.open();
207
208 outputFile.put(JMeta(argc, argv));
209 outputFile.put(parameters);
210 outputFile.put(mechanics);
211
212 try { // limit RAM usage
213 dynamic_cast<JTreeWriterObjectOutput<JEvent>&>(*outputFile).getTreeWriter().SetAutoFlush(autoflush);
214 dynamic_cast<JTreeWriterObjectOutput<JEvent>&>(*outputFile).getTreeWriter().SetCacheSize();
215 }
216 catch(const exception&) {}
217
218 int counter[] = { 0, 0 };
219
220 typedef deque<JEvent> buffer_type;
221
222 for (buffer_type zbuf; inputFile.hasNext(); ) {
223
224 STATUS(inputFile.getFilename() << '\r'); DEBUG(endl);
225
226 // read one file at a time
227
228 for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
229
230 const JEvent* evt = inputFile.next();
231
232 if (!evt->empty() && emitters.has(evt->getID())) {
233 if (utc(evt->begin()->getToA()) && utc(evt->rbegin()->getToA())) {
234 zbuf.push_back(*evt);
235 }
236 }
237 }
238
239 sort(zbuf.begin(), zbuf.end()); // sort according first time-of-emission
240
241 for (buffer_type::const_iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
242
243 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}
244
245 if (q == zbuf.end()) {
246
247 if (inputFile.hasNext()) {
248
249 zbuf.erase(zbuf.begin(), p); // remove processed data and continue reading
250
251 break;
252 }
253 }
254
255 // empty overlapping events and select unique transmissions
256
257 vector<JEvent> buffer(p, q); // super event containing copy of data
258
259 JEvent::overlap(buffer.begin(), buffer.end(), parameters.deadTime_s);
260
261 if (getNumberOfEmitters(buffer.begin(), buffer.end()) >= parameters.Nmin) {
262
263 h1.Fill((double) getNumberOfEmitters(buffer.begin(), buffer.end()));
264
265 const JWeight getWeight(buffer.begin(), buffer.end());
266
267 data_type data;
268
269 for (vector<JEvent>::iterator evt = buffer.begin(); evt != buffer.end(); ++evt) {
270
271 sort(evt->begin(), evt->end(), JKatoomba<>::compare);
272
273 JEvent::const_iterator __end = unique(evt->begin(), evt->end(), make_comparator(&JTransmission::getID, JComparison::eq()));
274
275 const JEmitter& emitter = emitters [evt->getID()];
276 const double weight = getWeight(evt->getID());
277
278 for (JEvent::const_iterator i = evt->begin(); i != __end; ++i) {
279
280 if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 &&
281 disable.count(JTransmission_t(-1, i->getID())) == 0) {
282
283 if (receivers.has(i->getID()) && geometry.hasLocation(receivers[i->getID()]) && i->getQ() >= parameters.Qmin * (parameters.Qmin <= 1.0 ? i->getW() : 1.0)) {
284
285 data.push_back(JHit(emitter,
286 distance(buffer.begin(),evt),
287 receivers[i->getID()],
288 i->getToA(),
289 parameters.sigma_s,
290 weight));
291 }
292 }
293 }
294 }
295
296 if (getMinimumNumberOfEmitters(data.begin(), data.end()) >= parameters.Nmin) {
297
298 for (data_type::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
299 HA[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
300 }
301
302 // fit
303
304 const auto result = katoomba(data.begin(), data.end());
305
306 for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
307 HB[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
308 }
309
310 if (debug >= debug_t) {
311
312 cout << "result:" << ' '
313 << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl
314 << result.value << endl;
315
316 for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
317 cout << "hit: " << *hit << ' ' << FIXED(9,3) << katoomba.evaluator(result.value, *hit) << endl;
318 }
319 }
320
321 hn.Fill(log10(katoomba.gandalf.numberOfIterations));
322 h0.Fill(result.chi2 / result.ndf);
323
324 // store results
325
326 if ((parameters.chi2perNDF > 0.0 && result.chi2 / result.ndf <= +parameters.chi2perNDF) ||
327 (parameters.chi2perNDF < 0.0 && result.chi2 / result.ndf >= -parameters.chi2perNDF)) {
328
329 const JEvt evt = getEvt(JHead(detector.getID(),
330 result.getTimeRange(),
331 data .size(),
332 result.size(),
333 result.value.getN(),
334 result.ndf,
335 result.chi2,
336 katoomba.gandalf.numberOfIterations),
337 result.value);
338
339 outputFile.put(evt);
340
341 if (selection.is_valid<JEvent>()) {
342
343 for (vector<JEvent>::iterator i = buffer.begin(); i != buffer.end(); ++i) {
344
345 JEvent out(*i); // event with fitted times of emission
346
347 const double toe = result.value.emission[JEKey(i->getID(), distance(buffer.begin(),i))].t1;
348
349 for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
350 hit->setToE(toe);
351 }
352
353 if (!out.empty()) {
354 outputFile.put(out);
355 }
356 }
357 }
358
359 counter[0] += 1;
360
361 } else {
362
363 WARNING(endl << "Event not written - chi2/NDF " << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl);
364
365 counter[1] += 1;
366 }
367
368 } else {
369
370 WARNING(endl << "Event not accepted - minimum number of emitters " << getMinimumNumberOfEmitters(data.begin(), data.end()) << endl);
371
372 counter[1] += 1;
373 }
374 }
375 }
376 }
377 STATUS(endl);
378
379 STATUS("Number of events written / rejected: " << counter[0] << " / " << counter[1] << endl);
380
381 if (!squash) {
382
383 JMultipleFileScanner<JMeta> io(inputFile);
384
385 io >> outputFile;
386 }
387
388 outputFile.put(h0);
389 outputFile.put(h1);
390 outputFile.put(hn);
391
392 for (JManager<int, TH2D>::const_iterator i = HA.begin(); i != HA.end(); ++i) { outputFile.put(*(i->second)); }
393 for (JManager<int, TH2D>::const_iterator i = HB.begin(); i != HB.end(); ++i) { outputFile.put(*(i->second)); }
394
395 outputFile.close();
396}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:74
WARNING
#define WARNING(A)
Definition JMessage.hh:65
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2140
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:57
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition build/antares-dataformat/PhysicsEvent.hh:434
JAANET::JHead
Monte Carlo run header.
Definition JHead.hh:1236
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JLANG::JException
General exception.
Definition JException.hh:25
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1697
JROOT::JManager
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
JROOT::JTreeWriterObjectOutput::getTreeWriter
JTreeWriter< T > & getTreeWriter()
Get TreeWriter.
Definition JTreeWriterObjectOutput.hh:78
JSUPPORT::JFileRecorder
Object writing to file.
Definition JFileRecorder.hh:44
JSUPPORT::JMultipleFileScanner
General purpose class for object reading from a list of file names.
Definition JMultipleFileScanner.hh:306
JSUPPORT::JSingleFileScanner
Object reading from a list of files.
Definition JSingleFileScanner.hh:172
JSUPPORT::JSingleFileScanner::hasNext
virtual bool hasNext() override
Check availability of next element.
Definition JSingleFileScanner.hh:257
JTOOLS::JHashCollection::has
bool has(const T &value) const
Test whether given value is present.
Definition JHashCollection.hh:309
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
JTOOLS::JRange::getLength
T getLength() const
Get length (difference between upper and lower limit).
Definition JRange.hh:289
JTOOLS::JRange::getLowerLimit
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
JACOUSTICS::getMinimumNumberOfEmitters
size_t getMinimumNumberOfEmitters(T __begin, T __end)
Get minimum number of emitters for any string in data.
Definition JAcoustics/JEvtToolkit.hh:69
JACOUSTICS::tripods_container
JContainer< std::vector< JTripod > > tripods_container
Definition JSydney.cc:79
JACOUSTICS::transmitters_container
JContainer< std::vector< JTransmitter > > transmitters_container
Definition JSydney.cc:81
JACOUSTICS::hydrophones_container
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition JSydney.cc:80
JACOUSTICS::getNumberOfEmitters
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
Definition JAcoustics/JEvtToolkit.hh:45
JACOUSTICS::getSoundVelocity
static const JSoundVelocity getSoundVelocity(1541.0, -17.0e-3, -2000.0)
Function object for velocity of sound.
JCOMPASS::getEvt
JEvt getEvt(const JHead &header, const JModel &model)
Get event.
Definition JCompass/JEvtToolkit.hh:71
JDETECTOR::getRangeOfFloors
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Definition JDetectorToolkit.hh:295
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:520
JDETECTOR::getNumberOfStrings
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Definition JStringCounter.hh:95
JEEP::debug_t
@ debug_t
debug
Definition JMessage.hh:29
JLANG::make_comparator
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Definition JLang/JComparator.hh:185
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:48
JRECONSTRUCTION::buffer_type
std::vector< JHitW0 > buffer_type
hits
Definition JPerth.cc:74
JRECONSTRUCTION::data_type
std::vector< event_type > data_type
Definition JPerth.cc:105
JRECONSTRUCTION::getWeight
JRECONSTRUCTION::JWeight getWeight
JTOOLS::result
return result
Definition JPolint.hh:862
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::detector
Detector file.
Definition JHead.hh:227
JACOUSTICS::JDetectorMechanics
Auxiliary data structure for mechanical model parameters with commented data.
Definition JMechanics.hh:38
JACOUSTICS::JEKey
Emitter key.
Definition JEKey.hh:36
JACOUSTICS::JEmitter
Acoustic emitter.
Definition JEmitter.hh:30
JACOUSTICS::JEvent::getID
int getID() const
Get emitter identifier.
Definition JAcoustics/JEvent.hh:152
JACOUSTICS::JEvent::getDetectorID
const int getDetectorID() const
Get detector identifier.
Definition JAcoustics/JEvent.hh:130
JACOUSTICS::JEvt
Acoustic event fit.
Definition JAcoustics/JEvt.hh:310
JACOUSTICS::JFitParameters::Qmin
double Qmin
minimal quality transmission
Definition JFitParameters.hh:62
JACOUSTICS::JFitParameters::deadTime_s
double deadTime_s
dead time between events [s]
Definition JFitParameters.hh:70
JACOUSTICS::JFitParameters::Nmin
size_t Nmin
minimum number of emitters
Definition JFitParameters.hh:66
JACOUSTICS::JFitParameters::sigma_s
double sigma_s
time-of-arrival resolution [s]
Definition JFitParameters.hh:64
JACOUSTICS::JFitParameters::Tmax_s
double Tmax_s
time window to combine events [s]
Definition JFitParameters.hh:65
JACOUSTICS::JFitParameters::chi2perNDF
double chi2perNDF
maximal chi2/NDF to store event
Definition JFitParameters.hh:69
JACOUSTICS::JGlobalfit
Global fit of prameterised detector geometry to acoustics data.
Definition JGlobalfit.hh:29
JACOUSTICS::JHit
Acoustics hit.
Definition JAcoustics/JHit.hh:35
JACOUSTICS::JKatoomba
Template definition of fit function of acoustic model.
Definition JKatoomba_t.hh:76
JACOUSTICS::JSoundVelocity
Implementation for depth dependend velocity of sound.
Definition JSoundVelocity.hh:37
JACOUSTICS::JSoundVelocity::set
JSoundVelocity & set(const double z0)
Set depth.
Definition JSoundVelocity.hh:62
JACOUSTICS::JTransmission_t
Acoustic transmission identifier.
Definition JTransmission_t.hh:28
JACOUSTICS::JTransmission::getID
int getID() const
Get identifier.
Definition JTransmission.hh:93
JACOUSTICS::JWeight
Auxiliary data structure to unify weights of acoustics data according to the number of pings per emit...
Definition JAcousticsToolkit.hh:55
JCALIBRATE::data_type
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition JFitK40.hh:107
JEEP::JContainer
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition JContainer.hh:42
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:67
JROOT::JROOTClassSelection
Auxiliary class for ROOT class selection.
Definition JROOTClassSelector.hh:132
JROOT::JROOTClassSelection::is_valid
bool is_valid() const
Get status of given data type.
Definition JROOTClassSelector.hh:255
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
JTOOLS::JHashMap
General purpose class for hash map of unique keys.
Definition JHashMap.hh:75
Generated by doxygen 1.12.0