Jpp 20.0.0-195-g190c9e876
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 JDetectorMechanics mechanics; // mechanical model data
80 JTimeRange T_s(-1.0e-4, +1.0e-4); // time range of background [s]
81 JFitParameters parameters; // fit parameters
82 bool unify; // unify weighing of pings
83 disable_container disable; // disable tansmissions
84 waveform_container waveforms; // emitter power and frequency
85 double background;
86 ULong_t seed;
87 int debug;
88
89 numberOfPings[WILDCARD] = NUMBER_OF_PINGS;
90 waveforms [WILDCARD] = EMITTER_WAVEFORM;
91
92 try {
93
94 JParser<> zap("Example application to test fit of model to simulated acoustic data.");
95
96 zap['a'] = make_field(detectorFile);
97 zap['o'] = make_field(outputFile);
98 zap['n'] = make_field(numberOfEvents);
99 zap['@'] = make_field(parameters) = JPARSER::initialised();
100 zap['N'] = make_field(numberOfPings, "number of pings per emitter") = JPARSER::initialised();
101 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
102 zap['T'] = make_field(tripods, "tripod data");
103 zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised();
104 zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised();
105 zap['u'] = make_field(unify, "unify weighing of pings");
106 zap['M'] = make_field(mechanics, "mechanics data") = JPARSER::initialised();
107 zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised();
108 zap['W'] = make_field(waveforms, "emitter power and frequency") = JPARSER::initialised();
109 zap['B'] = make_field(background, "background probability") = 0.0;
110 zap['S'] = make_field(seed) = 0;
111 zap['d'] = make_field(debug) = 1;
112
113 zap(argc, argv);
114 }
115 catch(const exception &error) {
116 FATAL(error.what() << endl);
117 }
118
119
120 gRandom->SetSeed(seed);
121
122 if (numberOfEvents <= 0) { FATAL("Invalid number of events " << numberOfEvents << endl); }
123
124 JDetector detector;
125
126 try {
127 load(detectorFile, detector);
128 }
129 catch(const JException& error) {
130 FATAL(error);
131 }
132
133 vector<JEmitter> emitters;
134
135 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
136 emitters.push_back(JEmitter(i->getID(),
137 i->getUTMPosition() - detector.getUTMPosition()));
138 }
139
140 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
141 try {
142 emitters.push_back(JEmitter(i->getID(),
143 i->getPosition() + detector.getModule(i->getLocation()).getPosition()));
144 }
145 catch(const exception&) {
146 continue; // if no string available, discard transmitter
147 }
148 }
149
150 if (detector.empty()) { FATAL("No modules in detector." << endl); }
151 if (emitters.empty()) { FATAL("No emitters in system." << endl); }
152
153
154 V.set(detector.getUTMZ()); // sound velocity at detector depth
155
156 const double D_m = -detector.getUTMZ(); // depth [m]
157
158 JGeometry geometry(detector, mechanics, hydrophones);
159 JGlobalfit katoomba(geometry, V, parameters);
160
161 JKatoomba<JGandalf>::debug = debug;
162 JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;
163
164 DEBUG(geometry);
165
166 typedef vector<JHit> data_type;
167
168
169 TH1D h0("cpu", NULL, 100, 0.0, 1.0e3);
170 TH1D h1("chi2/NDF", NULL, 100, 0.0, 5.0);
171
172 JManager<int, TH1D> H1(new TH1D("emitter[%]", NULL, 100, -5.0e-5, +5.0e-5));
173 JManager<int, TH2D> H2(new TH2D("string[%]", NULL, 100, -2.0, +2.0, 100, -2.0, +2.0));
174
175
176 for (int number_of_events = 0, count = 0; number_of_events != numberOfEvents; ++number_of_events) {
177
178 STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
179
180 JModel model; // randomised model data
181
182 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
183
184 if (model.string[module->getString()] == JMODEL::JString()) {
185
186 model.string[module->getString()] = JMODEL::JString(gRandom->Uniform(-2.0e-2, +2.0e-2),
187 gRandom->Uniform(-2.0e-2, +2.0e-2));
188 }
189 }
190
191 DEBUG("Model" << endl << model << endl);
192
193 // set module positions according actual model data
194
195 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
196 if (module->getFloor() != 0 && geometry.hasLocation(module->getLocation())) {
197 module->set(geometry[module->getString()].getPosition(model.string[module->getString()], module->getFloor()));
198 }
199 }
200
201
202 // generate hits
203
204 int minimum_number_of_pings = numeric_limits<int>::max();
205
206 for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
207 minimum_number_of_pings = min(minimum_number_of_pings, getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]));
208 }
209
210 data_type data;
211
212 for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
213
214 const int number_of_pings = getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]);
215
216 const double weight = (unify ? (double) minimum_number_of_pings / (double) number_of_pings : 1.0);
217
218 for (int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
219
220 ++count;
221
222 const double toe_s = emitter->getID() * 1000.0 + ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
223
224 model.emission[emitter->getID()][count] = toe_s;
225
226 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
227
228 if (disable.count(JTransmission_t(emitter->getID(), module->getID())) == 0) {
229
230 if (geometry.hasLocation(module->getLocation())) {
231
232 const JWaveform& waveform = getValue(waveforms, emitter->getID(), waveforms[WILDCARD]);
233
234 const double d_m = getDistance(module->getPosition(), emitter->getPosition());
235 const double toa_s = toe_s + V.getTime(d_m, emitter->getZ(), module->getZ());
236 const double Q = waveform.getQ(D_m, d_m);
237
238 if (Q >= parameters.Qmin) {
239
240 double t1_s = toa_s;
241
242 if (gRandom->Rndm() >= background)
243 t1_s = gRandom->Gaus(toa_s, parameters.sigma_s);
244 else
245 t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
246 toa_s + T_s.getUpperLimit());
247
248 const JHit hit(*emitter,
249 count,
250 module->getLocation(),
251 t1_s,
252 parameters.sigma_s,
253 weight);
254
255 DEBUG("hit: " << hit << ' ' << FIXED(7,1) << Q << endl);
256
257 data.push_back(hit);
258 }
259 }
260 }
261 }
262 }
263 }
264
265 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
266
267 const auto result = katoomba(data.begin(), data.end());
268
269 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
270
271 if (debug >= debug_t) {
272
273 cout << "result:" << ' '
274 << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl
275 << result.value << endl;
276
277 for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
278 cout << "hit: " << *hit << ' ' << FIXED(9,3) << katoomba.evaluator(result.value, *hit) << endl;
279 }
280 }
281
282 h0.Fill(chrono::duration_cast<chrono::milliseconds>(t1 - t0).count());
283 h1.Fill(result.chi2 / result.ndf);
284
285 for (JHashMap<int, JMODEL::JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
286
287 const double tx = (i->second.tx - result.value.string [i->first].tx) * 1.0e3; // mrad
288 const double ty = (i->second.ty - result.value.string [i->first].ty) * 1.0e3; // mrad
289
290 H2 ->Fill(tx, ty);
291 H2[i->first]->Fill(tx, ty);
292 }
293
294 for (JHashMap<JEKey, JMODEL::JEmission>::const_iterator i = model.emission.begin(); i != model.emission.end(); ++i) {
295
296 const double t1 = i->second.t1 - result.value.emission[i->first].t1;
297
298 H1 ->Fill(t1);
299 H1[i->first.getID()]->Fill(t1);
300 }
301 }
302 STATUS(endl);
303
304
305 TFile out(outputFile.c_str(), "recreate");
306
307 out << h0
308 << h1
309 << H1 << *H1
310 << H2 << *H2;
311
312 out.Write();
313 out.Close();
314}
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:74
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
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JROOT::JManager
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
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.
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:48
JRECONSTRUCTION::data_type
std::vector< event_type > data_type
Definition JPerth.cc:105
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::JDetectorMechanics
Auxiliary data structure for mechanical model parameters with commented data.
Definition JMechanics.hh:38
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
JACOUSTICS::JWaveform::getQ
double getQ(const double D_m, const double d_m) const
Get quality at given distance.
Definition JAcousticsSupportkit.hh:146
JCALIBRATE::data_type
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition JFitK40.hh:103
JEEP::JContainer
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Definition JContainer.hh:42
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.12.0