Jpp  18.2.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Enumerations | 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 <set>
#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/JTransmitter.hh"
#include "JDetector/JHydrophone.hh"
#include "JDetector/JModule.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_t.hh"
#include "JAcoustics/JTransmission_t.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.

Enumerations

enum  JFit_t
 Enumeration for fit algorithms. More...
 

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.

Enumeration Type Documentation

enum JFit_t

Enumeration for fit algorithms.

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

55  {
56  linear_t = 0, //!< linear fit
57  simplex_t, //!< simplex fit
58  gandalf_t //!< gandalf fit
59  };

Function Documentation

int main ( int  argc,
char **  argv 
)

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

70 {
71  using namespace std;
72  using namespace JPP;
73 
77  typedef JContainer< map<int, JWaveform> > waveform_container;
78  typedef JContainer< set<JTransmission_t> > disable_container;
79 
80  string detectorFile;
81  string outputFile;
82  int numberOfEvents;
83  map<int, int> numberOfPings;
84  JSoundVelocity V = getSoundVelocity; // default sound velocity
85  tripods_container tripods; // tripods
86  transmitters_container transmitters; // transmitters
87  hydrophones_container hydrophones; // hydrophones
88  JTimeRange T_s(-1.0e-4, +1.0e-4); // time range of background [s]
89  JFitParameters parameters; // fit parameters
90  int fit; // type of fit
91  bool unify; // unify weighing of pings
92  disable_container disable; // disable tansmissions
93  waveform_container waveforms; // emitter power and frequency
94  double background;
95  UInt_t seed;
96  int debug;
97 
98  numberOfPings[WILDCARD] = NUMBER_OF_PINGS;
99  waveforms [WILDCARD] = EMITTER_WAVEFORM;
100 
101  try {
102 
103  JParser<> zap("Example application to test fit of model to simulated acoustic data.");
104 
105  zap['a'] = make_field(detectorFile);
106  zap['o'] = make_field(outputFile);
107  zap['n'] = make_field(numberOfEvents);
109  zap['N'] = make_field(numberOfPings, "number of pings per emitter") = JPARSER::initialised();
110  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
111  zap['T'] = make_field(tripods, "tripod data");
112  zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
113  zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
114  zap['F'] = make_field(fit, "fit type") = linear_t, simplex_t, gandalf_t;
115  zap['u'] = make_field(unify, "unify weighing of pings");
116  zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
117  zap['W'] = make_field(waveforms, "emitter power and frequency") = JPARSER::initialised();
118  zap['B'] = make_field(background, "background probability") = 0.0;
119  zap['S'] = make_field(seed) = 0;
120  zap['d'] = make_field(debug) = 1;
121 
122  zap(argc, argv);
123  }
124  catch(const exception &error) {
125  FATAL(error.what() << endl);
126  }
127 
128 
129  gRandom->SetSeed(seed);
130 
131  if (numberOfEvents <= 0) { FATAL("Invalid number of events " << numberOfEvents << endl); }
132 
134 
135  try {
136  load(detectorFile, detector);
137  }
138  catch(const JException& error) {
139  FATAL(error);
140  }
141 
142  vector<JEmitter> emitters;
143 
144  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
145  emitters.push_back(JEmitter(i->getID(),
146  i->getUTMPosition() - detector.getUTMPosition()));
147  }
148 
149  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
150  try {
151  emitters.push_back(JEmitter(i->getID(),
152  i->getPosition() + detector.getModule(i->getLocation()).getPosition()));
153  }
154  catch(const exception&) {
155  continue; // if no string available, discard transmitter
156  }
157  }
158 
159  if (detector.empty()) { FATAL("No modules in detector." << endl); }
160  if (emitters.empty()) { FATAL("No emitters in system." << endl); }
161 
162 
163  V.set(detector.getUTMZ()); // sound velocity at detector depth
164 
165  const double D_m = -detector.getUTMZ(); // depth [m]
166 
167 
168  const JGeometry geometry(detector, hydrophones);
169 
170  DEBUG(geometry);
171 
172  JKatoomba<JEstimator> estimator(geometry, V, parameters.option);
173  JKatoomba<JAbstractMinimiser> evaluator(geometry, V, parameters.option);
174  JKatoomba<JSimplex> simplex (geometry, V, parameters.option);
175  JKatoomba<JGandalf> gandalf (geometry, V, parameters.option);
176 
177  simplex.estimator.reset(getMEstimator(parameters.mestimator));
178  gandalf.estimator.reset(getMEstimator(parameters.mestimator));
179 
180  simplex.debug = debug;
181  gandalf.debug = debug;
182 
184 
185 
186  TH1D h0("cpu", NULL, 100, 1.0, 7.0);
187  TH1D h1("chi2/NDF", NULL, 100, 0.0, 5.0);
188 
189  JManager<int, TH1D> H1(new TH1D("emitter[%]", NULL, 100, -5.0e-5, +5.0e-5));
190  JManager<int, TH2D> H2(new TH2D("string[%]", NULL, 100, -1.0, +1.0, 100, -1.0, +1.0));
191 
192 
193  JTimer timer;
194 
195  for (int number_of_events = 0, count = 0; number_of_events != numberOfEvents; ++number_of_events) {
196 
197  STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
198 
199  JModel model; // randomised model data
200 
201  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
202 
203  if (model.string[module->getString()] == JMODEL::JString()) {
204 
205  model.string[module->getString()] = JMODEL::JString(gRandom->Uniform(-2.0e-2, +2.0e-2),
206  gRandom->Uniform(-2.0e-2, +2.0e-2));
207  }
208  }
209 
210  DEBUG("Model" << endl << model << endl);
211 
212  // set module positions according actual model data
213 
214  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
215  if (geometry.hasLocation(module->getLocation())) {
216  module->set(estimator.geometry[module->getString()].getPosition(model.string[module->getString()], module->getFloor()));
217  }
218  }
219 
220 
221  // generate hits
222 
223  int minimum_number_of_pings = numeric_limits<int>::max();
224 
225  for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
226  minimum_number_of_pings = min(minimum_number_of_pings, getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]));
227  }
228 
230 
231  for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
232 
233  const int number_of_pings = getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]);
234 
235  const double weight = (unify ? (double) minimum_number_of_pings / (double) number_of_pings : 1.0);
236 
237  for (int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
238 
239  ++count;
240 
241  const double toe_s = ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
242 
243  model.emission[emitter->getID()][count] = toe_s;
244 
245  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
246 
247  if (disable.count(JTransmission_t(emitter->getID(), module->getID())) == 0) {
248 
249  if (geometry.hasLocation(module->getLocation())) {
250 
251  const JWaveform& waveform = getValue(waveforms, emitter->getID(), waveforms[WILDCARD]);
252 
253  const double d_m = getDistance(module->getPosition(), emitter->getPosition());
254  const double toa_s = toe_s + V.getTime(d_m, emitter->getZ(), module->getZ());
255  const double Q = waveform.getQ(D_m, d_m);
256 
257  if (Q >= parameters.Qmin) {
258 
259  double t1_s = toa_s;
260 
261  if (gRandom->Rndm() >= background)
262  t1_s = gRandom->Gaus(toa_s, parameters.sigma_s);
263  else
264  t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
265  toa_s + T_s.getUpperLimit());
266 
267  const JHit hit(*emitter,
268  count,
269  module->getLocation(),
270  t1_s,
271  parameters.sigma_s,
272  weight);
273 
274  DEBUG("hit: " << hit << ' ' << FIXED(7,1) << Q << endl);
275 
276  data.push_back(hit);
277  }
278  }
279  }
280  }
281  }
282  }
283 
284 
285  timer.reset();
286  timer.start();
287 
288  JModel result = estimator(data.begin(), data.end()); // pre-fit
289  double chi2 = numeric_limits<double>::max();
290 
291  switch (fit) {
292 
293  case linear_t:
294 
295  chi2 = evaluator(result, data.begin(), data.end());
296  break;
297 
298  case simplex_t:
299 
300  simplex.value = result; // start value
301 
302  chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
303  result = simplex.value;
304  break;
305 
306  case gandalf_t:
307 
308  gandalf.value = result; // start value
309 
310  chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
311  result = gandalf.value;
312  break;
313 
314  default:
315  break;
316  }
317 
318  timer.stop();
319 
320  double W = 0.0;
321 
322  for (vector<JHit>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
323  W += hit->getWeight();
324  }
325 
326  const int ndf = data.size() - result.getN();
327 
328  DEBUG("Final values" << endl
329  << FIXED(9,3) << chi2 << '/' << ndf << endl
330  << result << endl);
331 
332 
333  h0.Fill(log10((double) timer.usec_wall));
334  h1.Fill(chi2 / (double) (W - result.getN()));
335 
336  for (JHashMap<int, JMODEL::JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
337  H2[i->first]->Fill((i->second.tx - result.string [i->first].tx) * 1.0e3, // mrad
338  (i->second.ty - result.string [i->first].ty) * 1.0e3); // mrad
339  }
340 
341  for (JHashMap<JEKey, JMODEL::JEmission>::const_iterator i = model.emission.begin(); i != model.emission.end(); ++i) {
342  H1[i->first.getID()]->Fill(i->second.t1 - result.emission[i->first].t1);
343  }
344  }
345  STATUS(endl);
346 
347 
348  TFile out(outputFile.c_str(), "recreate");
349 
350  out << h0
351  << h1
352  << H1
353  << H2;
354 
355  out.Write();
356  out.Close();
357 }
Utility class to parse command line options.
Definition: JParser.hh:1514
size_t getN() const
Get number of fit parameters.
General exception.
Definition: JException.hh:24
Q(UTCMax_s-UTCMin_s)-livetime_s
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
JContainer< std::vector< JTransmitter > > transmitters_container
Definition: JSydney.cc:79
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
Template specialisation of fit function of acoustic model based on linear approximation.
Definition: JKatoomba_t.hh:495
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
Definition: JKatoomba_t.hh:642
*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 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
Acoustics hit.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
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_t.hh:735
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
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition: JSydney.cc:78
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:1989
set_variable E_E log10(E_{fit}/E_{#mu})"
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:32
JPosition3D getPosition(const Vec &pos)
Get position.
JACOUSTICS::JModel::emission_type emission
double getQ(const double D_m, const double d_m) const
Get quality at given distance.
JContainer< std::vector< JTripod > > tripods_container
Definition: JSydney.cc:77
Implementation for depth dependend velocity of sound.
#define FATAL(A)
Definition: JMessage.hh:67
JACOUSTICS::JModel::string_type string
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
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_t.hh:329
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
container_type::const_iterator const_iterator
Definition: JHashMap.hh:86
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203
do echo n Creating graphics for string $STRING for((FLOOR=$FIRST_FLOOR;$FLOOR<=$LAST_FLOOR;FLOOR+=1))
do set_variable DETECTOR_TXT $WORKDIR detector
Acoustic transmission identifier.
Utility class for emitter power and frequency.
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62