Jpp  test_elongated_shower_pde
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
software/JAcoustics/JKatoomba.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 
10 #include "TROOT.h"
11 #include "TFile.h"
12 #include "TH1D.h"
13 #include "TH2D.h"
14 
15 #include "JLang/JPredicate.hh"
16 #include "JLang/JComparator.hh"
17 #include "JLang/JComparison.hh"
18 
19 #include "JROOT/JManager.hh"
20 
21 #include "JDetector/JDetector.hh"
23 #include "JDetector/JTripod.hh"
24 #include "JDetector/JModule.hh"
25 #include "JDetector/JHydrophone.hh"
26 
27 #include "JTools/JHashMap.hh"
28 
31 #include "JSupport/JMeta.hh"
32 
34 #include "JAcoustics/JEmitter.hh"
36 #include "JAcoustics/JHit.hh"
38 #include "JAcoustics/JKatoomba.hh"
39 #include "JAcoustics/JEvent.hh"
40 #include "JAcoustics/JEvt.hh"
42 #include "JAcoustics/JSupport.hh"
44 
45 #include "Jeep/JContainer.hh"
46 #include "Jeep/JParser.hh"
47 #include "Jeep/JMessage.hh"
48 
49 
50 /**
51  * \file
52  *
53  * Application to make a global fit of the detector geometry to acoustic data.\n
54  * \author mdejong
55  */
56 int main(int argc, char **argv)
57 {
58  using namespace std;
59  using namespace JPP;
60 
61  typedef JContainer< vector<JTripod> > tripods_container;
62  typedef JContainer< vector<JHydrophone> > hydrophones_container;
63  typedef JContainer< set<JTransmission_t> > disable_container;
64 
67  string detectorFile;
68  JLimit_t& numberOfEvents = inputFile.getLimit();
69  JSoundVelocity V = getSoundVelocity; // default sound velocity
70  tripods_container tripods; // tripods
71  hydrophones_container hydrophones; // hydrophones
72  JFitParameters parameters; // fit parameters
73  int fit; // type of fit
74  bool unify; // unify weighing of pings
75  disable_container disable; // disable tansmissions
76  int debug;
77 
78  try {
79 
80  JParser<> zap("Application to fit position calibration model to acoustic data.");
81 
82  zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
83  zap['o'] = make_field(outputFile);
84  zap['n'] = make_field(numberOfEvents) = JLimit::max();
85  zap['a'] = make_field(detectorFile);
87  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
88  zap['T'] = make_field(tripods, "tripod data");
89  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
90  zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
91  zap['F'] = make_field(fit, "fit type") = linear_t, simplex_t, gandalf_t;
92  zap['u'] = make_field(unify, "unify weighing of pings");
93  zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
94  zap['d'] = make_field(debug) = 1;
95 
96  zap(argc, argv);
97  }
98  catch(const exception &error) {
99  FATAL(error.what() << endl);
100  }
101 
102 
103  cout.tie(&cerr);
104 
106 
107  try {
108  load(detectorFile, detector);
109  }
110  catch(const JException& error) {
111  FATAL(error);
112  }
113 
115 
116  JHashMap<int, JLocation> receivers;
117  JHashMap<int, JEmitter> emitters;
118 
119  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
120  receivers[i->getID()] = i->getLocation();
121  }
122 
123  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
124  emitters[i->getID()] = JEmitter(i->getID(),
125  i->getUTMPosition() - detector.getUTMPosition());
126  }
127 
128  V.set(detector.getUTMZ()); // sound velocity at detector depth
129 
130  JGeometry geometry(detector, hydrophones);
131 
132  DEBUG(geometry);
133 
134 
135  JKatoomba<JEstimator> estimator(geometry, V);
136  JKatoomba<JAbstractMinimiser> evaluator(geometry, V);
137  JKatoomba<JSimplex> simplex (geometry, V);
138  JKatoomba<JGandalf> gandalf (geometry, V);
139 
140  evaluator.estimator.reset(getMEstimator(parameters.mestimator));
141  simplex .estimator.reset(getMEstimator(parameters.mestimator));
142  gandalf .estimator.reset(getMEstimator(parameters.mestimator));
143 
144  simplex.debug = debug;
145  gandalf.debug = debug;
146 
147  if (parameters.fixStrings) {
148  JKatoomba_t::setOption(false);
149  }
150 
151  JSimplex <JModel> ::MAXIMUM_ITERATIONS = 10000;
153 
154  typedef JHit<JPDFGauss> hit_type;
155 
156  TH1D h0("chi2/NDF", NULL, 50000, 0.0, 1000.0);
157  TH1D h1("h1", NULL, 51, -0.5, 50.5);
158  TH1D hn("hn", NULL, 100, 0.0, 6.0);
159 
160  JManager<int, TH2D> HA(new TH2D("ha[%]", NULL,
162  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
163 
164  JManager<int, TH2D> HB(new TH2D("hb[%]", NULL,
166  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
167 
168  for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
169  HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
170  HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
171  }
172 
173  string oid; // detector identifier
174 
175  { // sort input files
176  map<double, string> zmap;
177 
178  for (const string& file_name : inputFile) {
179 
180  STATUS(file_name << '\r'); DEBUG(endl);
181 
182  for (JMultipleFileScanner<JEvent> in(file_name, 1); in.hasNext(); ) {
183 
184  const JEvent* evt = in.next();
185 
186  if (oid == "")
187  oid = evt->getOID();
188  else if (oid != evt->getOID()) // consistency check
189  FATAL("Invalid detector identifier " << evt->getOID() << " != " << oid << endl);
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 
206  outputFile.open();
207 
208  outputFile.put(JMeta(argc, argv));
209 
210  int counter[] = { 0, 0 };
211 
212  typedef deque<JEvent> buffer_type;
213 
214  for (buffer_type zbuf; inputFile.hasNext(); ) {
215 
216  STATUS(inputFile.getFilename() << '\r'); DEBUG(endl);
217 
218  // read one file at a time
219 
220  for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
221 
222  const JEvent* evt = inputFile.next();
223 
224  zbuf.push_back(*evt);
225  }
226 
227  sort(zbuf.begin(), zbuf.end()); // sort according first time-of-emission
228 
229  for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
230 
231  for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}
232 
233  if (q == zbuf.end()) {
234 
235  if (inputFile.hasNext()) {
236 
237  zbuf.erase(zbuf.begin(), p); // remove processed data and continue reading
238 
239  break;
240  }
241  }
242 
243  if (getNumberOfEmitters(p,q) >= parameters.Nmin) {
244 
245  h1.Fill((double) getNumberOfEmitters(p,q));
246 
247  map<int, int> numberOfPings;
248 
249  for (buffer_type::const_iterator i = p; i != q; ++i) {
250  numberOfPings[i->getID()] += 1;
251  }
252 
253  for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
254  DEBUG("Number of pings " << setw(2) << i->first << ' ' << setw(3) << i->second << endl);
255  }
256 
257  int minimum_number_of_pings = numeric_limits<int>::max();
258 
259  for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
260  minimum_number_of_pings = min(minimum_number_of_pings, i->second);
261  }
262 
263  vector<hit_type> data;
264 
265  for (buffer_type::iterator i = p; i != q; ++i) {
266 
267  const JEmitter& emitter = emitters[i->getID()];
268 
269  const double signal = (unify ? (double) minimum_number_of_pings / (double) numberOfPings[i->getID()] : 1.0);
270 
271  // select transmission with highest quality
272 
274 
275  for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
276  if (disable.count(JTransmission_t(i->getID(), hit->getID())) == 0 &&
277  disable.count(JTransmission_t(-1, hit->getID())) == 0) {
278  buffer[hit->getID()].push_back(*hit);
279  }
280  }
281 
282  for (map<int, vector<JTransmission> >::iterator ps = buffer.begin(); ps != buffer.end(); ++ps) {
283 
284  if (receivers.has(ps->first) && geometry.hasLocation(receivers[ps->first])) {
285 
286  sort(ps->second.begin(), ps->second.end(), make_comparator(&JTransmission::getQ, JComparison::gt()));
287 
288  if (ps->second.begin()->getQ() >= parameters.Qmin) {
289 
290  data.push_back(hit_type(emitter,
291  distance(p,i),
292  receivers[ps->first],
293  JPDFGauss(ps->second.begin()->getToA(), parameters.sigma_s, signal, parameters.background)));
294  }
295  }
296  }
297  }
298 
299  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
300  DEBUG("hit: " << *hit << endl);
301  }
302 
303 
304  JModel result = estimator(data.begin(), data.end()); // pre-fit
305 
306  DEBUG("prefit:" << endl
307  << result << endl);
308 
309 
310  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
311  HA[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
312  }
313 
314 
315  // remove outliers
316 
317  if (parameters.stdev > 0.0) {
318 
319  double chi2_old = evaluator(result, data.begin(), data.end());
320 
321  for ( ; ; ) {
322 
323  vector<hit_type>::iterator out = data.end();
324 
325  double xmax = 0.0;
326 
327  for (vector<hit_type>::iterator hit = data.begin(); hit != data.end(); ++hit) {
328 
329  const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->sigma;
330 
331  if (x > xmax) {
332  xmax = x;
333  out = hit;
334  }
335  }
336 
337  if (xmax > parameters.stdev) {
338 
339  vector<hit_type>::iterator __end = data.end();
340 
341  iter_swap(out, --__end);
342 
343  result = estimator(data.begin(), __end);
344 
345  double chi2_new = evaluator(result, data.begin(), __end);
346 
347  if (chi2_old - chi2_new > parameters.stdev * parameters.stdev) {
348 
349  DEBUG("Remove outlier " << __end->getEKey() << ' ' << __end->getLocation() << ' ' << xmax << endl);
350 
351  data.pop_back();
352 
353  chi2_old = chi2_new;
354 
355  continue;
356 
357  } else {
358 
359  result = estimator(data.begin(), ++__end);
360  }
361  }
362 
363  break;
364  }
365  }
366 
367  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
368  DEBUG("hit: " << *hit << endl);
369  }
370 
371  DEBUG("prefit:" << endl
372  << result << endl);
373 
374  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
375  HB[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
376  }
377 
378 
379  double chi2 = numeric_limits<double>::max();
380 
381  switch (fit) {
382 
383  case linear_t:
384 
385  result = estimator(data.begin(), data.end());
386  chi2 = evaluator(result, data.begin(), data.end());
387  break;
388 
389  case simplex_t:
390 
391  simplex.value = result; // start value
392 
393  chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
394  result = simplex.value;
395 
396  hn.Fill(log10(simplex.numberOfIterations));
397  break;
398 
399  case gandalf_t:
400 
401  gandalf.value = result; // start value
402 
403  chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
404  result = gandalf.value;
405 
406  hn.Fill(log10(gandalf.numberOfIterations));
407  break;
408 
409  default:
410  break;
411  }
412 
413  double W = 0.0;
414 
415  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
416  W += hit->getWeight();
417  }
418 
419  const int ndf = data.size() - result.getN();
420 
421  DEBUG("result:" << endl
422  << FIXED(9,3) << chi2 << '/' << ndf << endl
423  << result << endl);
424 
425  h0.Fill(chi2 / (W - result.getN()));
426 
427  // store results
428 
429  if (chi2 / ndf <= parameters.chi2perNDF) {
430 
431  const JEvt evt = getEvt(JHead(oid, data.begin()->getValue(), data.rbegin()->getValue(), ndf, (W - result.getN()), chi2), result);
432 
433  outputFile.put(evt);
434 
435  for (buffer_type::iterator i = p; i != q; ++i) {
436 
437  JEvent out(*i); // event with fitted times of emission
438 
439  const double toe = result.emitter[JEKey(i->getID(), distance(p,i))].t1;
440 
441  for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
442  hit->setToE(toe);
443  }
444 
445  outputFile.put(out);
446  }
447 
448  counter[0] += 1;
449 
450  } else {
451 
452  WARNING(endl << "Event not written: " << chi2 << '/' << ndf << endl);
453 
454  counter[1] += 1;
455  }
456  }
457  }
458  }
459  STATUS(endl);
460 
461  STATUS("Number of events written / rejected: " << counter[0] << " / " << counter[1] << endl);
462 
463  JMultipleFileScanner<JMeta> io(inputFile);
464 
465  io >> outputFile;
466 
467  outputFile.put(h0);
468  outputFile.put(h1);
469  outputFile.put(hn);
470 
471  for (JManager<int, TH2D>::const_iterator i = HA.begin(); i != HA.end(); ++i) { outputFile.put(*(i->second)); }
472  for (JManager<int, TH2D>::const_iterator i = HB.begin(); i != HB.end(); ++i) { outputFile.put(*(i->second)); }
473 
474  outputFile.close();
475 }
JDetector detector
Definition: JRunAnalyzer.hh:23
static int debug
debug level (default is off).
Definition: JMessage.hh:45
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
Acoustic hit.
size_t getN() const
Get number of fit parameters.
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
int main(int argc, char *argv[])
Definition: Main.cc:15
Greater than.
Definition: JComparison.hh:73
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
Sound velocity.
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
General purpose class for hash map of unique elements.
#define STATUS(A)
Definition: JMessage.hh:63
ROOT TTree parameter settings.
Detector data structure.
Definition: JDetector.hh:89
Recording of objects on file according a format that follows from the file name extension.
Template specialisation of fit function of acoustic model based on linear approximation.
Definition: JKatoomba.hh:285
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
Definition: JKatoomba.hh:385
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:151
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Acoustic event.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
Dynamic ROOT object management.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
int getIndex(const T &value) const
Get index of given value.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
string outputFile
Acoustic emitter.
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Data structure for detector geometry and calibration.
Acoustics hit.
Data structure for hydrophone.
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
Acoustic fit parameters.
Linear fit.
Definition: JKatoomba.hh:45
Model for fit to acoustics data.
Acoustic event fit.
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:466
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
Detector file.
Definition: JHead.hh:224
Acoustic event fit.
Acoustic emitter.
Definition: JEmitter.hh:27
Acoustics toolkit.
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
set_variable E_E log10(E_{fit}/E_{#mu})"
const std::string & getOID() const
Get detector identifier.
static int debug
debug level
Definition: JKatoomba.hh:599
return result
Definition: JPolint.hh:743
ROOT I/O of application specific meta data.
int debug
debug level
Definition: JSirene.cc:68
double getQ(const double D_m, const double f_kHz, const double d_m)
Get relative quality for given frequency at given distance.
int numberOfIterations
Definition: JSimplex.hh:242
JModel_t value
Definition: JSimplex.hh:240
General purpose messaging.
Monte Carlo run header.
Definition: JHead.hh:1164
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Implementation for depth dependend velocity of sound.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
JACOUSTICS::JModel::emitter_type emitter
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.
Utility class to parse command line options.
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
Definition: JKatoomba.hh:226
Acoustic transmission identifier.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
double getToA(const JModel &model, const JHit< JPDF_t > &hit) const
Get estimated time-of-arrival for given hit.
Definition: JKatoomba.hh:78
Emitter key.
Definition: JEKey.hh:32
bool hasLocation(const JLocation &location) const
Check if this detector has given location.
Definition: JGeometry.hh:545
Acoustic event.
Custom probability density function of time-of-arrival.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
Fit functions of acoustic model.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:166
Data structure for tripod.
then fatal Not enough tripods
Acoustic event fit.
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:42
Acoustic transmission identifier.
Simplex fit.
Definition: JKatoomba.hh:46
Template definition of fit function of acoustic model.
Definition: JKatoomba.hh:218
Container I/O.
do if[[!-f $ACOUSTICS_WORKDIR/${KEY}.txt]]
Definition: JAcoustics.sh:39
Data structure for optical module.
Gandalf fit.
Definition: JKatoomba.hh:47