Jpp  pmt_effective_area_update_2
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 <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 "JROOT/JManager.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/JFitParameters.hh"
#include "JAcoustics/JKatoomba.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.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 53 of file software/JAcoustics/JKatoomba.cc.

54 {
55  using namespace std;
56  using namespace JPP;
57 
58  typedef JContainer< vector<JTripod> > tripods_container;
59  typedef JContainer< vector<JHydrophone> > hydrophones_container;
60 
63  string detectorFile;
64  JLimit_t& numberOfEvents = inputFile.getLimit();
65  JSoundVelocity V = getSoundVelocity; // default sound velocity
66  tripods_container tripods; // tripods
67  hydrophones_container hydrophones; // hydrophones
68  JFitParameters parameters; // fit parameters
69  int fit; // type of fit
70  bool unify; // unify weighing of pings
71  int debug;
72 
73  try {
74 
75  JParser<> zap("Application to fit position calibration model to acoustic data.");
76 
77  zap['f'] = make_field(inputFile, "output of JAcousticEventBuilder[.sh]");
78  zap['o'] = make_field(outputFile);
79  zap['n'] = make_field(numberOfEvents) = JLimit::max();
80  zap['a'] = make_field(detectorFile);
82  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
83  zap['T'] = make_field(tripods, "tripod data");
84  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
85  zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
86  zap['F'] = make_field(fit, "fit type") = linear_t, simplex_t, gandalf_t;
87  zap['u'] = make_field(unify, "unify weighing of pings");
88  zap['d'] = make_field(debug) = 1;
89 
90  zap(argc, argv);
91  }
92  catch(const exception &error) {
93  FATAL(error.what() << endl);
94  }
95 
96 
97  cout.tie(&cerr);
98 
100 
101  try {
102  load(detectorFile, detector);
103  }
104  catch(const JException& error) {
105  FATAL(error);
106  }
107 
108 
109  JHashMap<int, JLocation> receivers;
110  JHashMap<int, JEmitter> emitters;
111 
112  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
113  receivers[i->getID()] = i->getLocation();
114  }
115 
116  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
117  emitters[i->getID()] = JEmitter(i->getID(),
118  i->getUTMPosition() - detector.getUTMPosition());
119  }
120 
121  V.set(detector.getUTMZ()); // sound velocity at detector depth
122 
123  JGeometry geometry(detector, hydrophones);
124 
125  DEBUG(geometry);
126 
127 
128  JKatoomba<JEstimator> estimator(geometry, V);
129  JKatoomba<JAbstractMinimiser> evaluator(geometry, V);
130  JKatoomba<JSimplex> simplex (geometry, V);
131  JKatoomba<JGandalf> gandalf (geometry, V);
132 
133  evaluator.estimator.reset(getMEstimator(parameters.mestimator));
134  simplex .estimator.reset(getMEstimator(parameters.mestimator));
135  gandalf .estimator.reset(getMEstimator(parameters.mestimator));
136 
137  simplex.debug = debug;
138  gandalf.debug = debug;
139 
140  if (parameters.fixStrings) {
141  JKatoomba_t::setOption(false);
142  }
143 
145 
146  typedef JHit<JPDFGauss> hit_type;
147 
148  TH1D h0("chi2/NDF", NULL, 5000, 0.0, 100.0);
149  TH1D h1("h1", NULL, 51, -0.5, 50.5);
150 
151  JManager<int, TH2D> HA(new TH2D("ha[%]", NULL,
153  getNumberOfFloors (detector) + 1, -0.5, getNumberOfFloors (detector) + 0.5));
154 
155  JManager<int, TH2D> HB(new TH2D("hb[%]", NULL,
157  getNumberOfFloors (detector) + 1, -0.5, getNumberOfFloors (detector) + 0.5));
158 
159  for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
160  HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
161  HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
162  }
163 
164  // input data
165 
166  string oid; // detector identifier
167 
168  vector<JEvent> zbuf;
169 
170  while (inputFile.hasNext()) {
171 
172  STATUS("input " << setw(6) << inputFile.getCounter() << '\r'); DEBUG(endl);
173 
174  const JEvent* evt = inputFile.next();
175 
176  if (oid == "") {
177  oid = evt->getOID();
178  }
179 
180  if (oid != evt->getOID()) { // consistency check
181  FATAL("Invalid detector identifier " << evt->getOID() << " != " << oid << endl);
182  }
183 
184  zbuf.push_back(*evt);
185  }
186  STATUS(endl);
187 
188  sort(zbuf.begin(), zbuf.end()); // sort according first time-of-emission
189 
190 
191  outputFile.open();
192 
193  outputFile.put(JMeta(argc, argv));
194 
195 
196  for (vector<JEvent>::const_iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
197 
198  STATUS("output " << setw(6) << distance(zbuf.cbegin(), p) << '\r'); DEBUG(endl);
199 
200  for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}
201 
202  if (getNumberOfEmitters(p,q) >= parameters.Nmin) {
203 
204  h1.Fill((double) getNumberOfEmitters(p,q));
205 
206  map<int, int> numberOfPings;
207 
208  for (vector<JEvent>::const_iterator i = p; i != q; ++i) {
209  numberOfPings[i->getID()] += 1;
210  }
211 
212  for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
213  DEBUG("Number of pings " << setw(2) << i->first << ' ' << setw(3) << i->second << endl);
214  }
215 
216  int minimum_number_of_pings = numeric_limits<int>::max();
217 
218  for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
219  minimum_number_of_pings = min(minimum_number_of_pings, i->second);
220  }
221 
222  vector<hit_type> data;
223 
224  for (vector<JEvent>::const_iterator i = p; i != q; ++i) {
225 
226  const JEmitter& emitter = emitters[i->getID()];
227 
228  const double signal = (unify ? (double) minimum_number_of_pings / (double) numberOfPings[i->getID()] : 1.0);
229 
230  // select transmission with highest quality
231 
232  map<int, vector<JTransmission> > buffer;
233 
234  for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
235  buffer[hit->getID()].push_back(*hit);
236  }
237 
238  for (map<int, vector<JTransmission> >::iterator ps = buffer.begin(); ps != buffer.end(); ++ps) {
239 
240  if (receivers.has(ps->first) && geometry.hasLocation(receivers[ps->first])) {
241 
242  sort(ps->second.begin(), ps->second.end(), make_comparator(&JTransmission::getQ, JComparison::gt()));
243 
244  if (ps->second.begin()->getQ() >= parameters.Qmin) {
245 
246  data.push_back(hit_type(emitter,
247  i->getCounter(),
248  receivers[ps->first],
249  JPDFGauss(ps->second.begin()->getToA(), parameters.sigma_s, signal, parameters.background)));
250  }
251  }
252  }
253  }
254 
255  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
256  DEBUG("hit: " << *hit << endl);
257  }
258 
259 
260 
261  JModel result = estimator(data.begin(), data.end()); // pre-fit
262 
263  DEBUG("prefit:" << endl
264  << result << endl);
265 
266 
267  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
268  HA[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
269  }
270 
271 
272  // remove outliers
273 
274  if (parameters.stdev > 0.0) {
275 
276  double chi2_old = evaluator(result, data.begin(), data.end());
277 
278  for ( ; ; ) {
279 
280  vector<hit_type>::iterator out = data.end();
281 
282  double xmax = 0.0;
283 
284  for (vector<hit_type>::iterator hit = data.begin(); hit != data.end(); ++hit) {
285 
286  const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->sigma;
287 
288  if (x > xmax) {
289  xmax = x;
290  out = hit;
291  }
292  }
293 
294  if (xmax > parameters.stdev) {
295 
296  vector<hit_type>::iterator __end = data.end();
297 
298  iter_swap(out, --__end);
299 
300  result = estimator(data.begin(), __end);
301 
302  double chi2_new = evaluator(result, data.begin(), __end);
303 
304  if (chi2_old - chi2_new > parameters.stdev * parameters.stdev) {
305 
306  DEBUG("Remove outlier " << __end->getLocation() << ' ' << xmax << endl);
307 
308  data.pop_back();
309 
310  chi2_old = chi2_new;
311 
312  continue;
313  }
314  }
315 
316  break;
317  }
318  }
319 
320  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
321  DEBUG("hit: " << *hit << endl);
322  }
323 
324  DEBUG("prefit:" << endl
325  << result << endl);
326 
327  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
328  HB[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
329  }
330 
331 
332  double chi2 = numeric_limits<double>::max();
333 
334  switch (fit) {
335 
336  case linear_t:
337 
338  result = estimator(data.begin(), data.end());
339  chi2 = evaluator(result, data.begin(), data.end());
340  break;
341 
342  case simplex_t:
343 
344  simplex.value = result; // start value
345 
346  chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
347  result = simplex.value;
348  break;
349 
350  case gandalf_t:
351 
352  gandalf.value = result; // start value
353 
354  chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
355  result = gandalf.value;
356  break;
357 
358  default:
359  break;
360  }
361 
362  double W = 0.0;
363 
364  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
365  W += hit->getWeight();
366  }
367 
368  const int ndf = data.size() - result.getN();
369 
370  DEBUG("result:" << endl
371  << FIXED(9,3) << chi2 << '/' << ndf << endl
372  << result << endl);
373 
374  h0.Fill(chi2 / (W - result.getN()));
375 
376  // store results
377 
378  const JEvt evt = getEvt(JHead(oid, data.begin()->getValue(), data.rbegin()->getValue(), ndf, (W - result.getN()), chi2), result);
379 
380  outputFile.put(evt);
381 
382  for (vector<JEvent>::const_iterator i = p; i != q; ++i) {
383 
384  JEvent out(*i); // event with fitted times of emission
385 
386  const double toe = result.emitter[JEKey(i->getID(), i->getCounter())].t1;
387 
388  for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
389  hit->setToE(toe);
390  }
391 
392  outputFile.put(out);
393  }
394  }
395  }
396  STATUS(endl);
397 
398  JMultipleFileScanner<JMeta> io(inputFile);
399 
400  io >> outputFile;
401 
402  outputFile.put(h0);
403  outputFile.put(h1);
404 
405  for (JManager<int, TH2D>::const_iterator i = HA.begin(); i != HA.end(); ++i) { outputFile.put(*(i->second)); }
406  for (JManager<int, TH2D>::const_iterator i = HB.begin(); i != HB.end(); ++i) { outputFile.put(*(i->second)); }
407 
408  outputFile.close();
409 }
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
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
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:81
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
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:446
string outputFile
Acoustics hit.
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
Linear fit.
Definition: JKatoomba.hh:45
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:466
Acoustic transmission.
do for((RUN=${RANGE%%-*};$RUN<=${RANGE##*-};RUN+=1))
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
Detector file.
Definition: JHead.hh:196
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:1961
const std::string & getOID() const
Get detector identifier.
return result
Definition: JPolint.hh:727
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
int debug
debug level
Definition: JSirene.cc:63
double getQ(const double D_m, const double f_kHz, const double d_m)
Get relative quality for given frequency at given distance.
Monte Carlo run header.
Definition: JHead.hh:1113
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
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:226
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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
Simplex fit.
Definition: JKatoomba.hh:46
do if[[!-f $ACOUSTICS_WORKDIR/${KEY}.txt]]
Definition: JAcoustics.sh:38
Gandalf fit.
Definition: JKatoomba.hh:47