44 int main(
int argc,
char **argv)
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); }
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) {
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) {
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);
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);
static int debug
debug level (default is off).
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)...
Dynamic ROOT object management.
Auxiliary data structure for floating point format specification.
Data structure for detector geometry and calibration.
Utility class to parse parameter values.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
Data structure for hydrophone.
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...
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
Auxiliary wrapper for I/O of container with optional comment (see JComment).
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
static int debug
debug level
Auxiliary class for CPU timing and usage.
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
General purpose messaging.
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++...
Utility class to parse command line options.
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
JPosition3D getPosition() const
Get position.
bool hasLocation(const JLocation &location) const
Check if this detector has given location.
Custom probability density function of time-of-arrival.
const JDetector & detector
Fit functions of acoustic model.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for tripod.
std::vector< double > weight
#define DEBUG(A)
Message macros.
Data structure for a composite optical module.
int main(int argc, char *argv[])