Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
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 <chrono>
#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/JGlobalfit.hh"
#include "JAcoustics/JTransmission_t.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

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

Author
mdejong

Definition in file examples/JAcoustics/JKatoomba.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

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

61{
62 using namespace std;
63 using namespace JPP;
64
65 typedef JContainer< vector<JTripod> > tripods_container;
66 typedef JContainer< vector<JTransmitter> > transmitters_container;
67 typedef JContainer< vector<JHydrophone> > hydrophones_container;
68 typedef JContainer< map<int, JWaveform> > waveform_container;
69 typedef JContainer< set<JTransmission_t> > disable_container;
70
71 string detectorFile;
72 string outputFile;
73 int numberOfEvents;
74 map<int, int> numberOfPings;
75 JSoundVelocity V = getSoundVelocity; // default sound velocity
76 tripods_container tripods; // tripods
77 transmitters_container transmitters; // transmitters
78 hydrophones_container hydrophones; // hydrophones
79 JTimeRange T_s(-1.0e-4, +1.0e-4); // time range of background [s]
80 JFitParameters parameters; // fit parameters
81 bool unify; // unify weighing of pings
82 disable_container disable; // disable tansmissions
83 waveform_container waveforms; // emitter power and frequency
84 double background;
85 ULong_t seed;
86 int debug;
87
88 numberOfPings[WILDCARD] = NUMBER_OF_PINGS;
89 waveforms [WILDCARD] = EMITTER_WAVEFORM;
90
91 try {
92
93 JParser<> zap("Example application to test fit of model to simulated acoustic data.");
94
95 zap['a'] = make_field(detectorFile);
96 zap['o'] = make_field(outputFile);
97 zap['n'] = make_field(numberOfEvents);
98 zap['@'] = make_field(parameters) = JPARSER::initialised();
99 zap['N'] = make_field(numberOfPings, "number of pings per emitter") = JPARSER::initialised();
100 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
101 zap['T'] = make_field(tripods, "tripod data");
102 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
103 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
104 zap['u'] = make_field(unify, "unify weighing of pings");
105 zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised();
106 zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
107 zap['W'] = make_field(waveforms, "emitter power and frequency") = JPARSER::initialised();
108 zap['B'] = make_field(background, "background probability") = 0.0;
109 zap['S'] = make_field(seed) = 0;
110 zap['d'] = make_field(debug) = 1;
111
112 zap(argc, argv);
113 }
114 catch(const exception &error) {
115 FATAL(error.what() << endl);
116 }
117
118
119 gRandom->SetSeed(seed);
120
121 if (numberOfEvents <= 0) { FATAL("Invalid number of events " << numberOfEvents << endl); }
122
123 JDetector detector;
124
125 try {
126 load(detectorFile, detector);
127 }
128 catch(const JException& error) {
129 FATAL(error);
130 }
131
132 vector<JEmitter> emitters;
133
134 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
135 emitters.push_back(JEmitter(i->getID(),
136 i->getUTMPosition() - detector.getUTMPosition()));
137 }
138
139 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
140 try {
141 emitters.push_back(JEmitter(i->getID(),
142 i->getPosition() + detector.getModule(i->getLocation()).getPosition()));
143 }
144 catch(const exception&) {
145 continue; // if no string available, discard transmitter
146 }
147 }
148
149 if (detector.empty()) { FATAL("No modules in detector." << endl); }
150 if (emitters.empty()) { FATAL("No emitters in system." << endl); }
151
152
153 V.set(detector.getUTMZ()); // sound velocity at detector depth
154
155 const double D_m = -detector.getUTMZ(); // depth [m]
156
157 JGeometry geometry(detector, hydrophones);
158 JGlobalfit katoomba(geometry, V, parameters);
159
160 JKatoomba<JGandalf>::debug = debug;
161 JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;
162
163 DEBUG(geometry);
164
165 typedef vector<JHit> data_type;
166
167
168 TH1D h0("cpu", NULL, 100, 0.0, 1.0e3);
169 TH1D h1("chi2/NDF", NULL, 100, 0.0, 5.0);
170
171 JManager<int, TH1D> H1(new TH1D("emitter[%]", NULL, 100, -5.0e-5, +5.0e-5));
172 JManager<int, TH2D> H2(new TH2D("string[%]", NULL, 100, -2.0, +2.0, 100, -2.0, +2.0));
173
174
175 for (int number_of_events = 0, count = 0; number_of_events != numberOfEvents; ++number_of_events) {
176
177 STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
178
179 JModel model; // randomised model data
180
181 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
182
183 if (model.string[module->getString()] == JMODEL::JString()) {
184
185 model.string[module->getString()] = JMODEL::JString(gRandom->Uniform(-2.0e-2, +2.0e-2),
186 gRandom->Uniform(-2.0e-2, +2.0e-2));
187 }
188 }
189
190 DEBUG("Model" << endl << model << endl);
191
192 // set module positions according actual model data
193
194 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
195 if (module->getFloor() != 0 && geometry.hasLocation(module->getLocation())) {
196 module->set(geometry[module->getString()].getPosition(model.string[module->getString()], module->getFloor()));
197 }
198 }
199
200
201 // generate hits
202
203 int minimum_number_of_pings = numeric_limits<int>::max();
204
205 for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
206 minimum_number_of_pings = min(minimum_number_of_pings, getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]));
207 }
208
209 data_type data;
210
211 for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
212
213 const int number_of_pings = getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]);
214
215 const double weight = (unify ? (double) minimum_number_of_pings / (double) number_of_pings : 1.0);
216
217 for (int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
218
219 ++count;
220
221 const double toe_s = emitter->getID() * 1000.0 + ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
222
223 model.emission[emitter->getID()][count] = toe_s;
224
225 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
226
227 if (disable.count(JTransmission_t(emitter->getID(), module->getID())) == 0) {
228
229 if (geometry.hasLocation(module->getLocation())) {
230
231 const JWaveform& waveform = getValue(waveforms, emitter->getID(), waveforms[WILDCARD]);
232
233 const double d_m = getDistance(module->getPosition(), emitter->getPosition());
234 const double toa_s = toe_s + V.getTime(d_m, emitter->getZ(), module->getZ());
235 const double Q = waveform.getQ(D_m, d_m);
236
237 if (Q >= parameters.Qmin) {
238
239 double t1_s = toa_s;
240
241 if (gRandom->Rndm() >= background)
242 t1_s = gRandom->Gaus(toa_s, parameters.sigma_s);
243 else
244 t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
245 toa_s + T_s.getUpperLimit());
246
247 const JHit hit(*emitter,
248 count,
249 module->getLocation(),
250 t1_s,
251 parameters.sigma_s,
252 weight);
253
254 DEBUG("hit: " << hit << ' ' << FIXED(7,1) << Q << endl);
255
256 data.push_back(hit);
257 }
258 }
259 }
260 }
261 }
262 }
263
264 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
265
266 const auto result = katoomba(data.begin(), data.end());
267
268 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
269
270 if (debug >= debug_t) {
271
272 cout << "result:" << ' '
273 << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl
274 << result.value << endl;
275
276 for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
277 cout << "hit: " << *hit << ' ' << FIXED(9,3) << katoomba.evaluator(result.value, *hit) << endl;
278 }
279 }
280
281 h0.Fill(chrono::duration_cast<chrono::milliseconds>(t1 - t0).count());
282 h1.Fill(result.chi2 / result.ndf);
283
284 for (JHashMap<int, JMODEL::JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
285
286 const double tx = (i->second.tx - result.value.string [i->first].tx) * 1.0e3; // mrad
287 const double ty = (i->second.ty - result.value.string [i->first].ty) * 1.0e3; // mrad
288
289 H2 ->Fill(tx, ty);
290 H2[i->first]->Fill(tx, ty);
291 }
292
293 for (JHashMap<JEKey, JMODEL::JEmission>::const_iterator i = model.emission.begin(); i != model.emission.end(); ++i) {
294
295 const double t1 = i->second.t1 - result.value.emission[i->first].t1;
296
297 H1 ->Fill(t1);
298 H1[i->first.getID()]->Fill(t1);
299 }
300 }
301 STATUS(endl);
302
303
304 TFile out(outputFile.c_str(), "recreate");
305
306 out << h0
307 << h1
308 << H1 << *H1
309 << H2 << *H2;
310
311 out.Write();
312 out.Close();
313}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JLANG::JException
General exception.
Definition JException.hh:24
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JACOUSTICS::tripods_container
JContainer< std::vector< JTripod > > tripods_container
Definition JSydney.cc:79
JACOUSTICS::transmitters_container
JContainer< std::vector< JTransmitter > > transmitters_container
Definition JSydney.cc:81
JACOUSTICS::hydrophones_container
JContainer< std::vector< JHydrophone > > hydrophones_container
Definition JSydney.cc:80
JACOUSTICS::getSoundVelocity
static const JSoundVelocity getSoundVelocity(1541.0, -17.0e-3, -2000.0)
Function object for velocity of sound.
JACOUSTICS::getMechanics
static JDetectorMechanics getMechanics
Function object to get string mechanics.
Definition JMechanics.hh:242
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:520
JEEP::getValue
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition JScale.hh:47
JEEP::debug_t
@ debug_t
debug
Definition JMessage.hh:29
JFIT::model
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition JGandalf.hh:57
JMATH::getDistance
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Definition JMathToolkit.hh:73
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JRECONSTRUCTION::data_type
std::vector< event_type > data_type
Definition JPerth.cc:82
JTOOLS::result
return result
Definition JPolint.hh:862
KM3NETDAQ::WILDCARD
static const char WILDCARD
Definition JDAQTags.hh:56
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::detector
Detector file.
Definition JHead.hh:227
JACOUSTICS::JEmitter
Acoustic emitter.
Definition JEmitter.hh:30
JACOUSTICS::JFitParameters::Qmin
double Qmin
minimal quality transmission
Definition JFitParameters.hh:62
JACOUSTICS::JFitParameters::sigma_s
double sigma_s
time-of-arrival resolution [s]
Definition JFitParameters.hh:64
JACOUSTICS::JGlobalfit
Global fit of prameterised detector geometry to acoustics data.
Definition JGlobalfit.hh:29
JACOUSTICS::JHit
Acoustics hit.
Definition JAcoustics/JHit.hh:35
JACOUSTICS::JKatoomba
Template definition of fit function of acoustic model.
Definition JKatoomba_t.hh:76
JACOUSTICS::JModel
Model for fit to acoustics data.
Definition JAcoustics/JModel.hh:480
JACOUSTICS::JSoundVelocity
Implementation for depth dependend velocity of sound.
Definition JSoundVelocity.hh:37
JACOUSTICS::JSoundVelocity::getTime
virtual double getTime(const double D_m, const double z1, const double z2) const override
Get propagation time of sound.
Definition JSoundVelocity.hh:147
JACOUSTICS::JSoundVelocity::set
JSoundVelocity & set(const double z0)
Set depth.
Definition JSoundVelocity.hh:62
JACOUSTICS::JTransmission_t
Acoustic transmission identifier.
Definition JTransmission_t.hh:28
JACOUSTICS::JWaveform
Utility class for emitter power and frequency.
Definition JAcousticsSupportkit.hh:103
JCALIBRATE::data_type
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition JFitK40.hh:103
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
JTOOLS::JHashMap::const_iterator
container_type::const_iterator const_iterator
Definition JHashMap.hh:86
Generated by doxygen 1.10.0