61{
64
70
71 string detectorFile;
73 int numberOfEvents;
81 bool unify;
82 disable_container disable;
83 waveform_container waveforms;
84 double background;
85 ULong_t seed;
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
101 zap[
'T'] =
make_field(tripods,
"tripod data");
104 zap[
'u'] =
make_field(unify,
"unify weighing of pings");
108 zap[
'B'] =
make_field(background,
"background probability") = 0.0;
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
124
125 try {
127 }
130 }
131
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;
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
154
155 const double D_m = -
detector.getUTMZ();
156
159
162
164
166
167
168 TH1D h0("cpu", NULL, 100, 0.0, 1.0e3);
169 TH1D h1("chi2/NDF", NULL, 100, 0.0, 5.0);
170
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
180
181 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
182
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
191
192
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
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
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
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,
252 weight);
253
254 DEBUG(
"hit: " << hit <<
' ' <<
FIXED(7,1) << Q << endl);
255
257 }
258 }
259 }
260 }
261 }
262 }
263
264 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
265
267
268 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
269
271
272 cout << "result:" << ' '
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());
283
285
286 const double tx = (i->second.tx -
result.value.string [i->first].tx) * 1.0e3;
287 const double ty = (i->second.ty -
result.value.string [i->first].ty) * 1.0e3;
288
289 H2 ->Fill(tx, ty);
290 H2[i->first]->Fill(tx, ty);
291 }
292
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 }
302
303
305
306 out << h0
307 << h1
308 << H1 << *H1
309 << H2 << *H2;
310
311 out.Write();
312 out.Close();
313}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Utility class to parse command line options.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
JContainer< std::vector< JTripod > > tripods_container
JContainer< std::vector< JTransmitter > > transmitters_container
JContainer< std::vector< JHydrophone > > hydrophones_container
static const JSoundVelocity getSoundVelocity(1541.0, -17.0e-3, -2000.0)
Function object for velocity of sound.
static JDetectorMechanics getMechanics
Function object to get string mechanics.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::vector< event_type > data_type
static const char WILDCARD
Auxiliary data structure for floating point format specification.
double Qmin
minimal quality transmission
double sigma_s
time-of-arrival resolution [s]
Global fit of prameterised detector geometry to acoustics data.
Template definition of fit function of acoustic model.
Model for fit to acoustics data.
Implementation for depth dependend velocity of sound.
virtual double getTime(const double D_m, const double z1, const double z2) const override
Get propagation time of sound.
JSoundVelocity & set(const double z0)
Set depth.
Acoustic transmission identifier.
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Auxiliary wrapper for I/O of container with optional comment (see JComment).
Empty structure for specification of parser element that is initialised (i.e. does not require input)...