Jpp  17.2.1-pre0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
software/JAcoustics/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/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/JKatoomba.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

int main ( int  argc,
char **  argv 
)

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

60 {
61  using namespace std;
62  using namespace JPP;
63 
64  typedef JContainer< vector<JTripod> > tripods_container;
65  typedef JContainer< vector<JTransmitter> > transmitters_container;
66  typedef JContainer< vector<JHydrophone> > hydrophones_container;
67  typedef JContainer< set<JTransmission_t> > disable_container;
68 
71  string detectorFile;
72  JLimit_t& numberOfEvents = inputFile.getLimit();
73  JSoundVelocity V = getSoundVelocity; // default sound velocity
74  tripods_container tripods; // tripods
75  transmitters_container transmitters; // transmitters
76  hydrophones_container hydrophones; // hydrophones
77  JFitParameters parameters; // fit parameters
78  bool unify; // unify weighing of pings
79  disable_container disable; // disable tansmissions
80  JROOTClassSelection selection = getROOTClassSelection<JTYPELIST<JEvent>::typelist>();
81  double Tmax_s; // deadtime [s]
82  int debug;
83 
84  try {
85 
86  JParser<> zap("Application to fit position calibration model to acoustic data.");
87 
88  zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
89  zap['o'] = make_field(outputFile);
90  zap['n'] = make_field(numberOfEvents) = JLimit::max();
91  zap['a'] = make_field(detectorFile);
93  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
94  zap['T'] = make_field(tripods, "tripod data");
95  zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
96  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
97  zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
98  zap['u'] = make_field(unify, "unify weighing of pings");
99  zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
100  zap['C'] = make_field(selection,
101  "Precede name of data structure by a '+' or '-' "
102  "to add or remove data types in the output, respectively."
103  "\nROOT wildcards are accepted.") = JPARSER::initialised();
104  zap['D'] = make_field(Tmax_s, "deadtime [s]") = 100.0e-3;
105  zap['d'] = make_field(debug) = 1;
106 
107  zap(argc, argv);
108  }
109  catch(const exception &error) {
110  FATAL(error.what() << endl);
111  }
112 
113 
114  cout.tie(&cerr);
115 
117 
118  try {
119  load(detectorFile, detector);
120  }
121  catch(const JException& error) {
122  FATAL(error);
123  }
124 
126 
127  JHashMap<int, JLocation> receivers;
128  JHashMap<int, JEmitter> emitters;
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(),
136  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(),
142  i->getPosition() + detector.getModule(i->getLocation()).getPosition());
143  }
144  catch(const exception&) {
145  continue; // if no module available, discard transmitter
146  }
147  }
148 
149  V.set(detector.getUTMZ()); // sound velocity at detector depth
150 
151  JGeometry geometry(detector, hydrophones);
152  JKatoomba_t katoomba(geometry, V, parameters);
153 
156 
159 
160  DEBUG(geometry);
161 
162  typedef JHit<JPDFGauss> hit_type;
163  typedef vector<hit_type> data_type;
164 
165  TH1D h0("chi2/NDF", NULL, 50000, 0.0, 1000.0);
166  TH1D h1("h1", NULL, 51, -0.5, 50.5);
167  TH1D hn("hn", NULL, 100, 0.0, 6.0);
168 
169  JManager<int, TH2D> HA(new TH2D("ha[%]", NULL,
171  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
172 
173  JManager<int, TH2D> HB(new TH2D("hb[%]", NULL,
175  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
176 
177  for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
178  HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
179  HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
180  }
181 
182  string oid; // detector identifier
183 
184  { // sort input files
185  map<double, string> zmap;
186 
187  for (const string& file_name : inputFile) {
188 
189  STATUS(file_name << '\r'); DEBUG(endl);
190 
191  for (JMultipleFileScanner<JEvent> in(file_name, 1); in.hasNext(); ) {
192 
193  const JEvent* evt = in.next();
194 
195  if (oid == "")
196  oid = evt->getOID();
197  else if (oid != evt->getOID()) // consistency check
198  FATAL("Invalid detector identifier " << evt->getOID() << " != " << oid << endl);
199 
200  if (!evt->empty()) {
201  zmap[evt->begin()->getToE()] = file_name;
202  }
203  }
204  }
205  STATUS(endl);
206 
207  inputFile.clear();
208 
209  for (map<double, string>::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
210  inputFile.push_back(i->second);
211  }
212  }
213 
214 
215  outputFile.open();
216 
217  outputFile.put(JMeta(argc, argv));
218  outputFile.put(parameters);
219 
220  try { // limit RAM usage
221  dynamic_cast<JTreeWriterObjectOutput<JEvent>&>(*outputFile).getTreeWriter().SetAutoFlush(5000);
222  }
223  catch(const exception&) {}
224 
225  int counter[] = { 0, 0 };
226 
227  typedef deque<JEvent> buffer_type;
228 
229  for (buffer_type zbuf; inputFile.hasNext(); ) {
230 
231  STATUS(inputFile.getFilename() << '\r'); DEBUG(endl);
232 
233  // read one file at a time
234 
235  for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
236 
237  const JEvent* evt = inputFile.next();
238 
239  if (!evt->empty() && emitters.has(evt->getID())) {
240  zbuf.push_back(*evt);
241  }
242  }
243 
244  sort(zbuf.begin(), zbuf.end()); // sort according first time-of-emission
245 
246  for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
247 
248  for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}
249 
250  if (q == zbuf.end()) {
251 
252  if (inputFile.hasNext()) {
253 
254  zbuf.erase(zbuf.begin(), p); // remove processed data and continue reading
255 
256  break;
257  }
258  }
259 
260  // empty overlapping events
261 
262  for (buffer_type::iterator __q = p, __p = __q++; __p != q && __q != q; __p = __q++) {
263 
264  if (__q->begin()->getToE() < __p->rbegin()->getToE() + Tmax_s) {
265 
266  __p->clear(); // clear first
267 
268  for (__p = __q++; __p != q && __q != q; __p = __q++) {
269 
270  if (__q->begin()->getToE() < __p->rbegin()->getToE() + Tmax_s)
271  __p->clear();
272  else
273  break;
274  }
275 
276  __p->clear(); // clear last
277  }
278  }
279 
280  if (getNumberOfEmitters(p,q) >= parameters.Nmin) {
281 
282  h1.Fill((double) getNumberOfEmitters(p,q));
283 
284  map<int, int> numberOfPings;
285 
286  for (buffer_type::const_iterator i = p; i != q; ++i) {
287  numberOfPings[i->getID()] += 1;
288  }
289 
290  for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
291  DEBUG("Number of pings " << setw(2) << i->first << ' ' << setw(3) << i->second << endl);
292  }
293 
294  int minimum_number_of_pings = numeric_limits<int>::max();
295 
296  for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
297  minimum_number_of_pings = min(minimum_number_of_pings, i->second);
298  }
299 
300  set<int> buffer;
301  data_type data;
302 
303  for (buffer_type::iterator evt = p; evt != q; ++evt) {
304 
305  sort(evt->begin(), evt->end(), compare);
306 
307  JEvent::iterator __end = unique(evt->begin(), evt->end(), make_comparator(&JTransmission::getID, JComparison::eq()));
308 
309  const JEmitter& emitter = emitters[evt->getID()];
310  const double signal = (unify ? (double) minimum_number_of_pings / (double) numberOfPings[evt->getID()] : 1.0);
311 
312  for (JEvent::const_iterator i = evt->begin(); i != __end; ++i) {
313 
314  if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 &&
315  disable.count(JTransmission_t(-1, i->getID())) == 0) {
316 
317  if (receivers.has(i->getID()) && geometry.hasLocation(receivers[i->getID()]) && i->getQ() >= parameters.Qmin) {
318 
319  data.push_back(hit_type(emitter,
320  distance(p,evt),
321  receivers[i->getID()],
322  JPDFGauss(i->getToA(), parameters.sigma_s, signal, parameters.background)));
323 
324  buffer.insert(evt->getID());
325  }
326  }
327  }
328  }
329 
330  if (buffer.size() >= parameters.Nmin) {
331 
332  for (data_type::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
333  HA[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
334  }
335 
336  // fit
337 
338  const auto result = katoomba(data.begin(), data.end());
339 
340  for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
341  HB[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
342  }
343 
344  if (debug >= debug_t) {
345 
346  cout << "result:" << ' '
347  << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl
348  << result.value << endl;
349 
350  for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
351  cout << "hit: " << *hit << ' ' << FIXED(9,3) << katoomba.evaluator(result.value, *hit) << endl;
352  }
353  }
354 
355  hn.Fill(log10(katoomba.gandalf.numberOfIterations));
356  h0.Fill(result.chi2 / result.ndf);
357 
358  // store results
359 
360  if ((parameters.chi2perNDF > 0.0 && result.chi2 / result.ndf <= +parameters.chi2perNDF) ||
361  (parameters.chi2perNDF < 0.0 && result.chi2 / result.ndf >= -parameters.chi2perNDF)) {
362 
364 
365  const JEvt evt = getEvt(JHead(oid,
366  range.getLowerLimit(),
367  range.getUpperLimit(),
368  data.size(),
369  result.value.getN(),
370  result.ndf,
371  result.chi2),
372  result.value);
373 
374  outputFile.put(evt);
375 
376  if (selection.is_valid<JEvent>()) {
377 
378  for (buffer_type::iterator i = p; i != q; ++i) {
379 
380  JEvent out(*i); // event with fitted times of emission
381 
382  const double toe = result.value.emitter[JEKey(i->getID(), distance(p,i))].t1;
383 
384  for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
385  hit->setToE(toe);
386  }
387 
388  outputFile.put(out);
389  }
390 
391  }
392 
393  counter[0] += 1;
394 
395  } else {
396 
397  WARNING(endl << "Event not written: " << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl);
398 
399  counter[1] += 1;
400  }
401  }
402  }
403  }
404  }
405  STATUS(endl);
406 
407  STATUS("Number of events written / rejected: " << counter[0] << " / " << counter[1] << endl);
408 
409  JMultipleFileScanner<JMeta> io(inputFile);
410 
411  io >> outputFile;
412 
413  outputFile.put(h0);
414  outputFile.put(h1);
415  outputFile.put(hn);
416 
417  for (JManager<int, TH2D>::const_iterator i = HA.begin(); i != HA.end(); ++i) { outputFile.put(*(i->second)); }
418  for (JManager<int, TH2D>::const_iterator i = HB.begin(); i != HB.end(); ++i) { outputFile.put(*(i->second)); }
419 
420  outputFile.close();
421 }
Worker class for fit function of acoustic model.
Definition: JKatoomba.hh:869
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1517
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
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.
Definition: JComparator.hh:185
static JDetectorMechanics getMechanics
Function object to get string mechanics.
Definition: JMechanics.hh:243
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:72
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
Auxiliary class for ROOT class selection.
JEvt getEvt(const JHead &header, const JModel &model)
Get event.
*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
Definition: diff-Tuna.sh:38
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:136
then fatal Number of tripods
Definition: JFootprint.sh:45
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
string outputFile
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Acoustics hit.
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
static const JSoundVelocity getSoundVelocity(1541.0,-17.0e-3,-2000.0)
Function object for velocity of sound.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
Detector file.
Definition: JHead.hh:226
Acoustic event fit.
Acoustic emitter.
Definition: JEmitter.hh:27
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition: JContainer.hh:39
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1993
set_variable E_E log10(E_{fit}/E_{#mu})"
const std::string & getOID() const
Get detector identifier.
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition: JVectorize.hh:54
return result
Definition: JPolint.hh:764
JPosition3D getPosition(const Vec &pos)
Get position.
static struct JACOUSTICS::@4 compare
Auxiliary data structure to sort transmissions.
Monte Carlo run header.
Definition: JHead.hh:1167
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Implementation for depth dependend velocity of sound.
#define FATAL(A)
Definition: JMessage.hh:67
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
JTreeWriter< T > & getTreeWriter()
Get TreeWriter.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
Emitter key.
Definition: JEKey.hh:32
Acoustic event.
Custom probability density function of time-of-arrival.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
do echo n Creating graphics for string $STRING for((FLOOR=$FIRST_FLOOR;$FLOOR<=$LAST_FLOOR;FLOOR+=1))
do set_variable DETECTOR_TXT $WORKDIR detector
int getID() const
Get identifier.
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
Definition: JCanberra.sh:46
Acoustic transmission identifier.
Template definition of fit function of acoustic model.
Definition: JKatoomba.hh:115
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62