Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
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
68 typedef JContainer< set<JTransmission_t> > disable_container;
69
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 JFitParameters parameters; // fit parameters
79 disable_container disable; // disable tansmissions
80 JROOTClassSelection selection = getROOTClassSelection<JTYPELIST<JEvent>::typelist>();
82 bool squash;
83 Long64_t autoflush;
84 int debug;
85
86 try {
87
88 JParser<> zap("Application to fit position calibration model to acoustic data.");
89
90 zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
91 zap['o'] = make_field(outputFile);
92 zap['n'] = make_field(numberOfEvents) = JLimit::max();
93 zap['a'] = make_field(detectorFile);
94 zap['@'] = make_field(parameters) = JPARSER::initialised();
95 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
96 zap['T'] = make_field(tripods, "tripod data");
97 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
98 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
99 zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
100 zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
101 zap['C'] = make_field(selection,
102 "Precede name of data structure by a '+' or '-' "
103 "to add or remove data types in the output, respectively."
104 "\nROOT wildcards are accepted.") = JPARSER::initialised();
105 zap['u'] = make_field(utc, "select UTC time range") = JRange<double>();
106 zap['q'] = make_field(squash, "squash meta data");
107 zap['R'] = make_field(autoflush) = 5000;
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
118
119 try {
120 load(detectorFile, detector);
121 }
122 catch(const JException& error) {
123 FATAL(error);
124 }
125
127
128 JHashMap<int, JLocation> receivers;
130
131 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
132 receivers[i->getID()] = i->getLocation();
133 }
134
135 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
136 emitters[i->getID()] = JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition());
137 }
138
139 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
140 try {
141 emitters[i->getID()] = JEmitter(i->getID(), i->getPosition() + detector.getModule(i->getLocation()).getPosition());
142 }
143 catch(const exception&) {} // if no module available, discard transmitter
144 }
145
146 V.set(detector.getUTMZ()); // sound velocity at detector depth
147
148 JGeometry geometry(detector, hydrophones);
149 JGlobalfit katoomba(geometry, V, parameters);
150
153
154 DEBUG(geometry);
155
156 typedef vector<JHit> data_type;
157
158 TH1D h0("chi2/NDF", NULL, 50000, 0.0, 1000.0);
159 TH1D h1("h1", NULL, 51, -0.5, 50.5);
160 TH1D hn("hn", NULL, 100, 0.0, 6.0);
161
162 JManager<int, TH2D> HA(new TH2D("ha[%]", NULL,
164 range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
165
166 JManager<int, TH2D> HB(new TH2D("hb[%]", NULL,
168 range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
169
170 for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
171 HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
172 HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
173 }
174
175 if (inputFile.size() > 1u) { // sort input files
176
178
179 for (const string& file_name : inputFile) {
180
181 STATUS(file_name << '\r'); DEBUG(endl);
182
183 for (JSingleFileScanner<JEvent> in(file_name, 1); in.hasNext(); ) {
184
185 const JEvent* evt = in.next();
186
187 if (detector.getID() != evt->getDetectorID()) {
188 FATAL("Invalid detector identifier " << evt->getDetectorID() << " != " << evt->getDetectorID() << endl);
189 }
190
191 if (!evt->empty()) {
192 zmap[evt->begin()->getToE()] = file_name;
193 }
194 }
195 }
196 STATUS(endl);
197
198 inputFile.clear();
199
200 for (map<double, string>::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
201 inputFile.push_back(i->second);
202 }
203 }
204
205 outputFile.open();
206
207 outputFile.put(JMeta(argc, argv));
208 outputFile.put(parameters);
209
210 try { // limit RAM usage
211 dynamic_cast<JTreeWriterObjectOutput<JEvent>&>(*outputFile).getTreeWriter().SetAutoFlush(autoflush);
212 dynamic_cast<JTreeWriterObjectOutput<JEvent>&>(*outputFile).getTreeWriter().SetCacheSize();
213 }
214 catch(const exception&) {}
215
216 int counter[] = { 0, 0 };
217
218 typedef deque<JEvent> buffer_type;
219
220 for (buffer_type zbuf; inputFile.hasNext(); ) {
221
222 STATUS(inputFile.getFilename() << '\r'); DEBUG(endl);
223
224 // read one file at a time
225
226 for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
227
228 const JEvent* evt = inputFile.next();
229
230 if (!evt->empty() && emitters.has(evt->getID())) {
231 if (utc(evt->begin()->getToA()) && utc(evt->rbegin()->getToA())) {
232 zbuf.push_back(*evt);
233 }
234 }
235 }
236
237 sort(zbuf.begin(), zbuf.end()); // sort according first time-of-emission
238
239 for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
240
241 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}
242
243 if (q == zbuf.end()) {
244
245 if (inputFile.hasNext()) {
246
247 zbuf.erase(zbuf.begin(), p); // remove processed data and continue reading
248
249 break;
250 }
251 }
252
253 JEvent::overlap(p, q, parameters.deadTime_s); // empty overlapping events
254
255 if (getNumberOfEmitters(p,q) >= parameters.Nmin) {
256
257 h1.Fill((double) getNumberOfEmitters(p,q));
258
259 const JWeight getWeight(p, q);
260
262
263 for (buffer_type::iterator evt = p; evt != q; ++evt) {
264
265 sort(evt->begin(), evt->end(), JKatoomba<>::compare);
266
267 JEvent::iterator __end = unique(evt->begin(), evt->end(), make_comparator(&JTransmission::getID, JComparison::eq()));
268
269 const JEmitter& emitter = emitters [evt->getID()];
270 const double weight = getWeight(evt->getID());
271
272 for (JEvent::const_iterator i = evt->begin(); i != __end; ++i) {
273
274 if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 &&
275 disable.count(JTransmission_t(-1, i->getID())) == 0) {
276
277 if (receivers.has(i->getID()) && geometry.hasLocation(receivers[i->getID()]) && i->getQ() >= parameters.Qmin * (parameters.Qmin <= 1.0 ? i->getW() : 1.0)) {
278
279 data.push_back(JHit(emitter,
280 distance(p,evt),
281 receivers[i->getID()],
282 i->getToA(),
283 parameters.sigma_s,
284 weight));
285 }
286 }
287 }
288 }
289
290 if (getMinimumNumberOfEmitters(data.begin(), data.end()) >= parameters.Nmin) {
291
292 for (data_type::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
293 HA[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
294 }
295
296 // fit
297
298 const auto result = katoomba(data.begin(), data.end());
299
300 for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
301 HB[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
302 }
303
304 if (debug >= debug_t) {
305
306 cout << "result:" << ' '
307 << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl
308 << result.value << endl;
309
310 for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
311 cout << "hit: " << *hit << ' ' << FIXED(9,3) << katoomba.evaluator(result.value, *hit) << endl;
312 }
313 }
314
315 hn.Fill(log10(katoomba.gandalf.numberOfIterations));
316 h0.Fill(result.chi2 / result.ndf);
317
318 // store results
319
320 if ((parameters.chi2perNDF > 0.0 && result.chi2 / result.ndf <= +parameters.chi2perNDF) ||
321 (parameters.chi2perNDF < 0.0 && result.chi2 / result.ndf >= -parameters.chi2perNDF)) {
322
323 const JEvt evt = getEvt(JHead(detector.getID(),
324 result.getTimeRange(),
325 data .size(),
326 result.size(),
327 result.value.getN(),
328 result.ndf,
329 result.chi2,
330 katoomba.gandalf.numberOfIterations),
331 result.value);
332
333 outputFile.put(evt);
334
335 if (selection.is_valid<JEvent>()) {
336
337 for (buffer_type::iterator i = p; i != q; ++i) {
338
339 JEvent out(*i); // event with fitted times of emission
340
341 const double toe = result.value.emission[JEKey(i->getID(), distance(p,i))].t1;
342
343 for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
344 hit->setToE(toe);
345 }
346
347 if (!out.empty()) {
348 outputFile.put(out);
349 }
350 }
351 }
352
353 counter[0] += 1;
354
355 } else {
356
357 WARNING(endl << "Event not written - chi2/NDF " << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl);
358
359 counter[1] += 1;
360 }
361
362 } else {
363
364 WARNING(endl << "Event not accepted - minimum number of emitters " << getMinimumNumberOfEmitters(data.begin(), data.end()) << endl);
365
366 counter[1] += 1;
367 }
368 }
369 }
370 }
371 STATUS(endl);
372
373 STATUS("Number of events written / rejected: " << counter[0] << " / " << counter[1] << endl);
374
375 if (!squash) {
376
377 JMultipleFileScanner<JMeta> io(inputFile);
378
379 io >> outputFile;
380 }
381
382 outputFile.put(h0);
383 outputFile.put(h1);
384 outputFile.put(hn);
385
386 for (JManager<int, TH2D>::const_iterator i = HA.begin(); i != HA.end(); ++i) { outputFile.put(*(i->second)); }
387 for (JManager<int, TH2D>::const_iterator i = HB.begin(); i != HB.end(); ++i) { outputFile.put(*(i->second)); }
388
389 outputFile.close();
390}
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 WARNING(A)
Definition JMessage.hh:65
#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
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Monte Carlo run header.
Definition JHead.hh:1236
Detector data structure.
Definition JDetector.hh:96
General exception.
Definition JException.hh:24
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
JTreeWriter< T > & getTreeWriter()
Get TreeWriter.
Object writing to file.
General purpose class for object reading from a list of file names.
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.
Range of values.
Definition JRange.hh:42
T getLength() const
Get length (difference between upper and lower limit).
Definition JRange.hh:289
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
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.
JEvt getEvt(const JHead &header, const JModel &model)
Get event.
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
@ debug_t
debug
Definition JMessage.hh:29
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
std::vector< event_type > data_type
Definition JPerth.cc:81
return result
Definition JPolint.hh:862
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Detector file.
Definition JHead.hh:227
Emitter key.
Definition JEKey.hh:36
Acoustic emitter.
Definition JEmitter.hh:30
int getID() const
Get emitter identifier.
const int getDetectorID() const
Get detector identifier.
Acoustic event fit.
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]
double chi2perNDF
maximal chi2/NDF to store event
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...
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition JFitK40.hh:103
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition JContainer.hh:42
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary class for ROOT class selection.
bool is_valid() const
Get status of given data type.
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