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
examples/JAcoustics/JKatoomba.cc File Reference

Example application to test fit of model to simulated acoustic data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <algorithm>
#include "TRandom3.h"
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "JLang/JPredicate.hh"
#include "JLang/JSTDToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JHydrophone.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JAcousticsSupportkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JKatoomba.hh"
#include "Jeep/JContainer.hh"
#include "Jeep/JTimer.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

Example application to test fit of model to simulated acoustic data.

Author
mdejong

Definition in file examples/JAcoustics/JKatoomba.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 57 of file examples/JAcoustics/JKatoomba.cc.

58 {
59  using namespace std;
60  using namespace JPP;
61 
62  typedef JContainer< vector<JTripod> > tripods_container;
63  typedef JContainer< vector<JHydrophone> > hydrophones_container;
64  typedef JContainer< map<int, JWaveform> > waveform_container;
65 
66  string detectorFile;
67  string outputFile;
68  int numberOfEvents;
69  map<int, int> numberOfPings;
70  JSoundVelocity V = getSoundVelocity; // default sound velocity
71  tripods_container tripods; // tripods
72  hydrophones_container hydrophones; // hydrophones
73  JTimeRange T_s(-1.0e-4, +1.0e-4); // time range of background [s]
74  JFitParameters parameters; // fit parameters
75  int fit; // type of fit
76  bool unify; // unify weighing of pings
77  waveform_container waveforms; // emitter power and frequency
78  UInt_t seed;
79  int debug;
80 
81  numberOfPings[WILDCARD] = NUMBER_OF_PINGS;
82  waveforms [WILDCARD] = EMITTER_WAVEFORM;
83 
84  try {
85 
86  JParser<> zap("Example application to test fit of model to simulated acoustic data.");
87 
88  zap['a'] = make_field(detectorFile);
89  zap['o'] = make_field(outputFile);
90  zap['n'] = make_field(numberOfEvents);
92  zap['N'] = make_field(numberOfPings, "number of pings per emitter") = JPARSER::initialised();
93  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
94  zap['T'] = make_field(tripods, "tripod data");
95  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
96  zap['F'] = make_field(fit, "fit type") = linear_t, simplex_t, gandalf_t;
97  zap['u'] = make_field(unify, "unify weighing of pings");
98  zap['W'] = make_field(waveforms, "emitter power and frequency") = JPARSER::initialised();
99  zap['S'] = make_field(seed) = 0;
100  zap['d'] = make_field(debug) = 1;
101 
102  zap(argc, argv);
103  }
104  catch(const exception &error) {
105  FATAL(error.what() << endl);
106  }
107 
108 
109  cout.tie(&cerr);
110 
111  gRandom->SetSeed(seed);
112 
113  if (numberOfEvents <= 0) { FATAL("Invalid number of events " << numberOfEvents << endl); }
114 
115 
117 
118  try {
119  load(detectorFile, detector);
120  }
121  catch(const JException& error) {
122  FATAL(error);
123  }
124 
125  vector<JEmitter> emitters;
126 
127  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
128  emitters.push_back(JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition()));
129  }
130 
131  if (detector.empty()) { FATAL("No modules in detector." << endl); }
132  if (emitters.empty()) { FATAL("No emitters in system." << endl); }
133 
134 
135  V.set(detector.getUTMZ()); // sound velocity at detector depth
136 
137  const double D_m = -detector.getUTMZ(); // depth [m]
138 
139 
140  const JGeometry geometry(detector, hydrophones);
141 
142  DEBUG(geometry);
143 
144 
145  JKatoomba<JEstimator> estimator(geometry, V);
146  JKatoomba<JAbstractMinimiser> evaluator(geometry, V);
147  JKatoomba<JSimplex> simplex (geometry, V);
148  JKatoomba<JGandalf> gandalf (geometry, V);
149 
150  simplex.estimator.reset(getMEstimator(parameters.mestimator));
151  gandalf.estimator.reset(getMEstimator(parameters.mestimator));
152 
153  simplex.debug = debug;
154  gandalf.debug = debug;
155 
157 
158 
159  typedef JHit<JPDFGauss> hit_type;
160 
161 
162  TH1D h0("cpu", NULL, 100, 1.0, 7.0);
163  TH1D h1("chi2/NDF", NULL, 100, 0.0, 5.0);
164 
165  JManager<int, TH1D> H1(new TH1D("emitter[%]", NULL, 100, -5.0e-5, +5.0e-5));
166  JManager<int, TH2D> H2(new TH2D("string[%]", NULL, 100, -1.0, +1.0, 100, -1.0, +1.0));
167 
168 
169  JTimer timer;
170 
171  for (int number_of_events = 0, count = 0; number_of_events != numberOfEvents; ++number_of_events) {
172 
173  STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
174 
175  JModel model; // randomised model data
176 
177  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
178 
179  if (model.string[module->getString()] == JMODEL::JString()) {
180 
181  model.string[module->getString()] = JMODEL::JString(gRandom->Uniform(-2.0e-2, +2.0e-2),
182  gRandom->Uniform(-2.0e-2, +2.0e-2));
183  }
184  }
185 
186  DEBUG("Model" << endl << model << endl);
187 
188  // set module positions according actual model data
189 
190  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
191  if (geometry.hasLocation(module->getLocation())) {
192  module->set(estimator.detector[module->getString()].getPosition(model.string[module->getString()], module->getFloor()));
193  }
194  }
195 
196 
197  // generate hits
198 
199  int minimum_number_of_pings = numeric_limits<int>::max();
200 
201  for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
202  minimum_number_of_pings = min(minimum_number_of_pings, getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]));
203  }
204 
205  vector<hit_type> data;
206 
207  for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
208 
209  const int number_of_pings = getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]);
210 
211  const double signal = (unify ? (double) minimum_number_of_pings / (double) number_of_pings : 1.0);
212 
213  for (int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
214 
215  ++count;
216 
217  const double toe_s = ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
218 
219  model.emitter[emitter->getID()][count] = toe_s;
220 
221  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
222 
223  if (geometry.hasLocation(module->getLocation())) {
224 
225  const JWaveform& waveform = getValue(waveforms, emitter->getID(), waveforms[WILDCARD]);
226 
227  const double d_m = getDistance(module->getPosition(), emitter->getPosition());
228  const double toa_s = toe_s + V.getTime(d_m, emitter->getZ(), module->getZ());
229  const double Q = waveform.getQ(D_m, d_m);
230 
231  if (Q >= parameters.Qmin) {
232 
233  double t1_s = toa_s;
234 
235  if (gRandom->Rndm() >= parameters.background)
236  t1_s = gRandom->Gaus(toa_s, parameters.sigma_s);
237  else
238  t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
239  toa_s + T_s.getUpperLimit());
240 
241  const hit_type hit(*emitter,
242  count,
243  module->getLocation(),
244  JPDFGauss(t1_s, parameters.sigma_s, signal, parameters.background));
245 
246  DEBUG("hit: " << hit << ' ' << FIXED(7,1) << Q << endl);
247 
248  data.push_back(hit);
249  }
250  }
251  }
252  }
253  }
254 
255 
256  timer.reset();
257  timer.start();
258 
259  JModel result = estimator(data.begin(), data.end()); // pre-fit
260  double chi2 = numeric_limits<double>::max();
261 
262  switch (fit) {
263 
264  case linear_t:
265 
266  chi2 = evaluator(result, data.begin(), data.end());
267  break;
268 
269  case simplex_t:
270 
271  simplex.value = result; // start value
272 
273  chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
274  result = simplex.value;
275  break;
276 
277  case gandalf_t:
278 
279  gandalf.value = result; // start value
280 
281  chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
282  result = gandalf.value;
283  break;
284 
285  default:
286  break;
287  }
288 
289  timer.stop();
290 
291  double W = 0.0;
292 
293  for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
294  W += hit->getWeight();
295  }
296 
297  const int ndf = data.size() - result.getN();
298 
299  DEBUG("Final values" << endl
300  << FIXED(9,3) << chi2 << '/' << ndf << endl
301  << result << endl);
302 
303 
304  h0.Fill(log10((double) timer.usec_wall));
305  h1.Fill(chi2 / (double) (W - result.getN()));
306 
307  for (JHashMap<int, JMODEL::JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
308  H2[i->first]->Fill((i->second.tx - result.string [i->first].tx) * 1.0e3, // mrad
309  (i->second.ty - result.string [i->first].ty) * 1.0e3); // mrad
310  }
311 
312  for (JHashMap<JEKey, JMODEL::JEmitter>::const_iterator i = model.emitter.begin(); i != model.emitter.end(); ++i) {
313  H1[i->first.getID()]->Fill(i->second.t1 - result.emitter[i->first].t1);
314  }
315  }
316  STATUS(endl);
317 
318 
319  TFile out(outputFile.c_str(), "recreate");
320 
321  out << h0
322  << h1
323  << H1
324  << H2;
325 
326  out.Write();
327  out.Close();
328 }
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
Q(UTCMax_s-UTCMin_s)-livetime_s
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
do rm f tmp H1
#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
*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
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.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Linear fit.
Definition: JKatoomba.hh:45
Model for fit to acoustics data.
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
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 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
return result
Definition: JPolint.hh:727
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:32
double getQ(const double D_m, const double d_m) const
Get quality at given distance.
int debug
debug level
Definition: JSirene.cc:63
Implementation for velocity of sound.
#define FATAL(A)
Definition: JMessage.hh:67
JACOUSTICS::JModel::string_type string
JACOUSTICS::JModel::emitter_type emitter
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
std::vector< int > count
Definition: JAlgorithm.hh:180
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
Custom probability density function of time-of-arrival.
container_type::const_iterator const_iterator
Definition: JHashMap.hh:86
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:166
do set_variable DETECTOR_TXT $WORKDIR detector
Utility class for emitter power and frequency.
Simplex fit.
Definition: JKatoomba.hh:46
Gandalf fit.
Definition: JKatoomba.hh:47