57 tripods_container tripods;
58 hydrophones_container hydrophones;
61 double background = 1.0e-4;
68 numberOfPings[-1] = 11;
79 JParser<> zap(
"Example application to test fit of model to simukated acoustic data.");
90 zap[
'u'] =
make_field(unique,
"one ping per cycle");
96 catch(
const exception &error) {
97 FATAL(error.what() << endl);
103 gRandom->SetSeed(seed);
105 if (numberOfEvents <= 0) {
FATAL(
"Invalid number of events " << numberOfEvents << endl); }
118 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
119 emitters.push_back(
JEmitter(i->getID(), i->getUTMPosition() -
detector.getUTMPosition()));
122 if (
detector.empty()) {
FATAL(
"No modules in detector." << endl); }
123 if (emitters.empty()) {
FATAL(
"No emitters in system." << endl); }
139 simplex.debug =
debug;
140 gandalf.debug =
debug;
148 TH1D h0(
"cpu", NULL, 100, 1.0, 7.0);
149 TH1D
h1(
"chi2/NDF", NULL, 100, 0.0, 5.0);
157 for (
int number_of_events = 0,
count = 0; number_of_events != numberOfEvents; ++number_of_events) {
159 STATUS(
"event: " << setw(10) << number_of_events <<
'\r');
DEBUG(endl);
163 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
168 gRandom->Uniform(-2.0e-2, +2.0e-2));
176 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
177 if (geometry.hasLocation(module->getLocation())) {
178 module->set(estimator.detector[module->getString()].getPosition(model.
string[module->getString()], module->getFloor()));
189 int number_of_pings = numberOfPings[-1];
191 if (numberOfPings.find(emitter->getID()) != numberOfPings.end()) {
192 number_of_pings = numberOfPings[emitter->getID()];
195 int weight = numeric_limits<int>::max();
198 if (i->second < weight) {
203 const double signal = (unique ? (double) weight / (
double) number_of_pings : 1.0);
205 for (
int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
209 const double toe_s = ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
213 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
215 if (geometry.hasLocation(module->getLocation())) {
217 const double toa_s = toe_s + V.getTime(
getDistance(module->getPosition(), emitter->getPosition()), emitter->getZ(), module->getZ());
221 if (gRandom->Rndm() >= background)
222 t1_s = gRandom->Gaus(toa_s,
sigma_s);
224 t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
225 toa_s + T_s.getUpperLimit());
227 data.push_back(hit_type(*emitter,
229 module->getLocation(),
236 DEBUG(
"Model" << endl << model << endl);
239 DEBUG(
"hit: " << *hit << endl);
247 double chi2 = numeric_limits<double>::max();
253 chi2 = evaluator(result, data.begin(), data.end());
260 chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
261 result = simplex.value;
268 chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
269 result = gandalf.value;
281 W += hit->getWeight();
284 const int ndf = data.size() - model.
getN();
286 DEBUG(
"Final values" << endl
287 <<
FIXED(9,3) << chi2 <<
'/' << ndf << endl
291 h0.Fill(log10((
double) timer.usec_wall));
292 h1.Fill(chi2 / (
double) (W - model.
getN()));
295 H2[i->first]->Fill((i->second.tx - result.
string [i->first].tx) * 1.0e3,
296 (i->second.ty - result.
string [i->first].ty) * 1.0e3);
300 H1[i->first.getID()]->Fill(i->second.t1 - result.
emitter[i->first].t1);
Utility class to parse command line options.
size_t getN() const
Get number of fit parameters.
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
Template specialisation of fit function of acoustic model based on linear approximation.
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
Utility class to parse parameter values.
then for HISTOGRAM in h0 h1
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary data structure for floating point format specification.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Model for fit to acoustics data.
static const JSoundVelocity getSoundVelocity(1541.0,-17.0e-3,-2000.0)
Function object for velocity of sound.
Template specialisation of fit function of acoustic model based on JGandalf minimiser.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Auxiliary wrapper for I/O of container with optional comment (see JComment).
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Auxiliary class for CPU timing and usage.
Implementation for velocity of sound.
JACOUSTICS::JModel::string_type string
JACOUSTICS::JModel::emitter_type emitter
then $DIR JKatoomba a $DETECTOR o $WORKDIR katoomba root T $TRIPOD n sigma_s
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++...
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
Custom probability density function of time-of-arrival.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
do set_variable DETECTOR_TXT $WORKDIR detector
do for((RUN=${RANGE%%-*};$RUN<=${RANGE##*-};RUN+=1))
std::vector< double > weight
#define DEBUG(A)
Message macros.