82 disable_container disable;
83 waveform_container waveforms;
88 numberOfPings[
WILDCARD] = NUMBER_OF_PINGS;
89 waveforms [
WILDCARD] = EMITTER_WAVEFORM;
93 JParser<> zap(
"Example application to test fit of model to simulated acoustic data.");
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;
114 catch(
const exception &error) {
115 FATAL(error.what() << endl);
119 gRandom->SetSeed(seed);
121 if (numberOfEvents <= 0) {
FATAL(
"Invalid number of events " << numberOfEvents << endl); }
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()));
139 for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
141 emitters.push_back(
JEmitter(i->getID(),
142 i->getPosition() +
detector.getModule(i->getLocation()).getPosition()));
144 catch(
const exception&) {
149 if (
detector.empty()) {
FATAL(
"No modules in detector." << endl); }
150 if (emitters.empty()) {
FATAL(
"No emitters in system." << endl); }
155 const double D_m = -
detector.getUTMZ();
168 TH1D h0(
"cpu", NULL, 100, 0.0, 1.0e3);
169 TH1D h1(
"chi2/NDF", NULL, 100, 0.0, 5.0);
175 for (
int number_of_events = 0, count = 0; number_of_events != numberOfEvents; ++number_of_events) {
177 STATUS(
"event: " << setw(10) << number_of_events <<
'\r');
DEBUG(endl);
181 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
186 gRandom->Uniform(-2.0e-2, +2.0e-2));
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()));
203 int minimum_number_of_pings = numeric_limits<int>::max();
206 minimum_number_of_pings = min(minimum_number_of_pings,
getValue(numberOfPings, emitter->getID(), numberOfPings[
WILDCARD]));
213 const int number_of_pings =
getValue(numberOfPings, emitter->getID(), numberOfPings[
WILDCARD]);
215 const double weight = (unify ? (double) minimum_number_of_pings / (
double) number_of_pings : 1.0);
217 for (
int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
221 const double toe_s = emitter->getID() * 1000.0 + ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
223 model.emission[emitter->getID()][count] = toe_s;
225 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
227 if (disable.count(
JTransmission_t(emitter->getID(), module->getID())) == 0) {
229 if (geometry.hasLocation(module->getLocation())) {
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);
237 if (Q >= parameters.
Qmin) {
241 if (gRandom->Rndm() >= background)
242 t1_s = gRandom->Gaus(toa_s, parameters.
sigma_s);
244 t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
245 toa_s + T_s.getUpperLimit());
247 const JHit hit(*emitter,
249 module->getLocation(),
254 DEBUG(
"hit: " << hit <<
' ' <<
FIXED(7,1) << Q << endl);
264 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
268 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
272 cout <<
"result:" <<
' '
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;
281 h0.Fill(chrono::duration_cast<chrono::milliseconds>(t1 - t0).count());
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;
290 H2[i->first]->Fill(tx, ty);
295 const double t1 = i->second.t1 -
result.value.emission[i->first].t1;
298 H1[i->first.getID()]->Fill(t1);
#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.
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.
JSoundVelocity & set(const double z0)
Set depth.
virtual double getTime(const double D_m, const double z1, const double z2) const override
Get propagation time of sound.
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)...