Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
software/JAcoustics/JKatoomba.cc File Reference

Application to fit position calibration model to acoustic data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#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 "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JHydrophone.hh"
#include "JTools/JHashMap.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/JKatoomba.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"
#include "Jeep/JProperties.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 fit position calibration model to acoustic data.

Author
mdejong

Definition in file software/JAcoustics/JKatoomba.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

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

52 {
53  using namespace std;
54  using namespace JPP;
55 
56  typedef JContainer< vector<JTripod> > tripods_container;
57  typedef JContainer< vector<JHydrophone> > hydrophones_container;
58 
61  string detectorFile;
62  JLimit_t& numberOfEvents = inputFile.getLimit();
63  JSoundVelocity V = getSoundVelocity; // default sound velocity
64  tripods_container tripods; // tripods
65  hydrophones_container hydrophones; // hydrophones
66  map<int, JMechanics> mechanics;
67  int mestimator = EM_LINEAR; // M-estimator
68  double sigma_s = 50.0e-6; // time-of-arrival resolution [s]
69  double background = 1.0e-4; // probability background transmission
70  double Tmax_s = 300.0; // time window [s]
71  size_t Nmin = 3; // minimum number of emitters
72  double stdev = 5.0; // standard deviation for outlier removal
73  int fit; // type of fit
74  bool unique; // one ping per cycle
75  int debug;
76 
77  try {
78 
79  JProperties properties;
80 
81  properties.insert(gmake_property(mestimator));
82  properties.insert(gmake_property(sigma_s));
83  properties.insert(gmake_property(background));
84  properties.insert(gmake_property(Tmax_s));
85  properties.insert(gmake_property(Nmin));
86  properties.insert(gmake_property(stdev));
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(properties) = JPARSER::initialised();
95  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
96  zap['T'] = make_field(tripods, "tripod data");
97  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
98  zap['M'] = make_field(mechanics, "mechanics data") = JPARSER::initialised();
99  zap['F'] = make_field(fit, "fit type") = linear_t, simplex_t, gandalf_t;
100  zap['u'] = make_field(unique, "one ping per cycle");
101  zap['d'] = make_field(debug) = 1;
102 
103  zap(argc, argv);
104  }
105  catch(const exception &error) {
106  FATAL(error.what() << endl);
107  }
108 
109 
110  cout.tie(&cerr);
111 
113 
114  try {
115  load(detectorFile, detector);
116  }
117  catch(const JException& error) {
118  FATAL(error);
119  }
120 
121 
122  JHashMap<int, JLocation> receivers;
123  JHashMap<int, JEmitter> emitters;
124 
125  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
126  receivers[i->getID()] = i->getLocation();
127  }
128 
129  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
130  emitters[i->getID()] = JEmitter(i->getID(),
131  i->getUTMPosition() - detector.getUTMPosition());
132  }
133 
134 
135  JGeometry geometry(detector, hydrophones);
136 
137  for (JGeometry::iterator i = geometry.begin(); i != geometry.end(); ++i) {
138  i->second.mechanics = mechanics[i->first];
139  }
140 
141  DEBUG(geometry);
142 
143 
144  JKatoomba<JEstimator> estimator(geometry, V[detector.getUTMZ()]);
145  JKatoomba<JAbstractMinimiser> evaluator(geometry, V[detector.getUTMZ()]);
146  JKatoomba<JSimplex> simplex (geometry, V[detector.getUTMZ()]);
147  JKatoomba<JGandalf> gandalf (geometry, V[detector.getUTMZ()]);
148 
149  simplex.estimator.reset(getMEstimator(mestimator));
150  gandalf.estimator.reset(getMEstimator(mestimator));
151 
152  simplex.debug = debug;
153  gandalf.debug = debug;
154 
156 
157  typedef JHit<JPDFGauss> hit_type;
158 
159  TH1D h0("chi2/NDF", NULL, 5000, 0.0, 100.0);
160 
161  TH2D ha("ha", NULL,
164 
165  TH2D hb("hb", NULL,
168 
169  for (Int_t i = 1; i <= ha.GetXaxis()->GetNbins(); ++i) {
170  ha.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
171  hb.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
172  }
173 
174  // input data
175 
176  string oid; // detector identifier
177 
178  vector<JEvent> zbuf;
179 
180  while (inputFile.hasNext()) {
181 
182  STATUS("input " << setw(6) << inputFile.getCounter() << '\r'); DEBUG(endl);
183 
184  const JEvent* evt = inputFile.next();
185 
186  if (oid == "") {
187  oid = evt->getOID();
188  }
189 
190  if (oid != evt->getOID()) { // consistency check
191  FATAL("Invalid detector identifier " << evt->getOID() << " != " << oid << endl);
192  }
193 
194  zbuf.push_back(*evt);
195  }
196  STATUS(endl);
197 
198  sort(zbuf.begin(), zbuf.end());
199 
200 
201  outputFile.open();
202 
203  outputFile.put(JMeta(argc, argv));
204 
205 
206  for (vector<JEvent>::const_iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
207 
208  STATUS("output " << setw(6) << distance(zbuf.cbegin(), p) << '\r'); DEBUG(endl);
209 
210  for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + Tmax_s; ) {}
211 
212  if (getNumberOfEmitters(p,q) >= Nmin) {
213 
214  map<int, int> numberOfPings;
215 
216  for (vector<JEvent>::const_iterator i = p; i != q; ++i) {
217  numberOfPings[i->getID()] += 1;
218  }
219 
220  int weight = numeric_limits<int>::max();
221 
222  for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
223  DEBUG("Number of pings " << setw(2) << i->first << ' ' << setw(3) << i->second << endl);
224  }
225 
226  for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
227  if (i->second < weight) {
228  weight = i->second;
229  }
230  }
231 
232  vector<hit_type> data;
233 
234  for (vector<JEvent>::const_iterator i = p; i != q; ++i) {
235 
236  {
237  const JEmitter& emitter = emitters[i->getID()];
238 
239  const double signal = (unique ? (double) weight / (double) numberOfPings[i->getID()] : 1.0);
240 
241  // select first transmission
242 
243  map<int, vector<JTransmission> > buffer;
244 
245  for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
246  buffer[hit->getID()].push_back(*hit);
247  }
248 
249  for (map<int, vector<JTransmission> >::iterator ps = buffer.begin(); ps != buffer.end(); ++ps) {
250 
251  if (receivers.has(ps->first) && geometry.hasLocation(receivers[ps->first])) {
252 
253  sort(ps->second.begin(), ps->second.end(), make_comparator(&JTransmission::getToA));
254 
255  data.push_back(hit_type(emitter,
256  i->getCounter(),
257  receivers[ps->first],
258  JPDFGauss(ps->second.begin()->getToA(), sigma_s, signal, background)));
259  }
260  }
261  }
262  }
263 
264  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
265  DEBUG("hit: " << *hit << endl);
266  }
267 
268 
269 
270  JModel result = estimator(data.begin(), data.end()); // pre-fit
271 
272  DEBUG("prefit:" << endl
273  << result << endl);
274 
275 
276  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
277  ha.Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
278  }
279 
280 
281  // remove outliers
282 
283  if (stdev > 0.0) {
284 
285  for ( ; ; ) {
286 
287  vector<hit_type>::iterator out = data.end();
288 
289  double xmax = 0.0;
290 
291  for (vector<hit_type>::iterator hit = data.begin(); hit != data.end(); ++hit) {
292 
293  const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->sigma;
294 
295  if (x > xmax) {
296  xmax = x;
297  out = hit;
298  }
299  }
300 
301  if (xmax > stdev) {
302 
303  DEBUG("Remove outlier " << out->getLocation() << ' ' << xmax << endl);
304 
305  data.erase(out);
306 
307  result = estimator(data.begin(), data.end());
308 
309  } else {
310 
311  break;
312  }
313  }
314  }
315 
316  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
317  hb.Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
318  }
319 
320 
321  double chi2 = numeric_limits<double>::max();
322 
323  switch (fit) {
324 
325  case linear_t:
326 
327  result = estimator(data.begin(), data.end());
328  chi2 = evaluator(result, data.begin(), data.end());
329  break;
330 
331  case simplex_t:
332 
333  simplex.value = result; // start value
334 
335  chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
336  result = simplex.value;
337  break;
338 
339  case gandalf_t:
340 
341  gandalf.value = result; // start value
342 
343  chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
344  result = gandalf.value;
345  break;
346 
347  default:
348  break;
349  }
350 
351  double W = 0.0;
352 
353  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
354  W += hit->getWeight();
355  }
356 
357  const int ndf = data.size() - result.getN();
358 
359  DEBUG("result:" << endl
360  << FIXED(9,3) << chi2 << '/' << ndf << endl
361  << result << endl);
362 
363  h0.Fill(chi2 / (W - result.getN()));
364 
365  // store results
366 
367  const JEvt evt = getEvt(JHead(oid, data.begin()->getValue(), data.rbegin()->getValue(), ndf, chi2), result);
368 
369  outputFile.put(evt);
370 
371  for (vector<JEvent>::const_iterator i = p; i != q; ++i) {
372 
373  JEvent out(*i);
374 
375  const double toe = result.emitter[JEKey(i->getID(), i->getCounter())].t1;
376 
377  for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
378  hit->setToE(toe);
379  }
380 
381  outputFile.put(out);
382  }
383  }
384  }
385  STATUS(endl);
386 
387  JMultipleFileScanner<JMeta> io(inputFile);
388 
389  io >> outputFile;
390 
391  outputFile.put(h0);
392  outputFile.put(ha);
393  outputFile.put(hb);
394 
395  outputFile.close();
396 }
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:1500
size_t getN() const
Get number of fit parameters.
General exception.
Definition: JException.hh:23
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
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:71
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:80
Template specialisation of fit function of acoustic model based on linear approximation.
Definition: JKatoomba.hh:250
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
Definition: JKatoomba.hh:345
Utility class to parse parameter values.
Definition: JProperties.hh:496
JEvt getEvt(const JHead &header, const JModel &model)
Get event.
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:151
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:445
string outputFile
Acoustics hit.
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
Model for fit to acoustics data.
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.
Template specialisation of fit function of acoustic model based on JGandalf minimiser.
Definition: JKatoomba.hh:426
Acoustic transmission.
Detector file.
Definition: JHead.hh:196
Acoustic event fit.
Acoustic emitter.
Definition: JEmitter.hh:27
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
Definition: JKatoomba.hh:94
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:1961
const std::string & getOID() const
Get detector identifier.
return result
Definition: JPolint.hh:727
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
do if[[!-f $ACOUSTICS_WORKDIR/${KEY}.txt &&-f $JPP_LIB/${KEY}_ ${(l:8::0::0:) ID}.txt]]
Definition: JAcoustics.sh:35
int debug
debug level
Definition: JSirene.cc:63
Monte Carlo run header.
Definition: JHead.hh:1091
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Implementation for velocity of sound.
#define FATAL(A)
Definition: JMessage.hh:67
JACOUSTICS::JModel::emitter_type emitter
then $DIR JKatoomba a $DETECTOR o $WORKDIR katoomba root T $TRIPOD n sigma_s
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.
Simple fit method based on Powell&#39;s algorithm, see reference: Numerical Recipes in C++...
Definition: JSimplex.hh:42
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
Definition: JKatoomba.hh:191
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
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:166
do set_variable DETECTOR_TXT $WORKDIR detector
do for((RUN=${RANGE%%-*};$RUN<=${RANGE##*-};RUN+=1))
std::vector< double > weight
Definition: JAlgorithm.hh:428
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
container_type::iterator iterator
Definition: JHashMap.hh:87