60int main(
int argc,
char **argv)
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) {
185 model.string[
module->getString()] = JMODEL::JString(gRandom->Uniform(-2.0e-2, +2.0e-2),
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();
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]));
211 for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
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) {
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);
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();
266 const auto result = katoomba(data.begin(), data.end());
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) {
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);