43 static const int WILDCARD = -1;
45 static const int NUMBER_OF_PINGS = 11;
46 static const JWaveform EMITTER_WAVEFORM(2.5e6,
57 int main(
int argc,
char **argv)
71 tripods_container tripods;
72 hydrophones_container hydrophones;
77 waveform_container waveforms;
81 numberOfPings[WILDCARD] = NUMBER_OF_PINGS;
82 waveforms [WILDCARD] = EMITTER_WAVEFORM;
86 JParser<> zap(
"Example application to test fit of model to simulated acoustic data.");
97 zap[
'u'] =
make_field(unify,
"unify weighing of pings");
104 catch(
const exception &error) {
105 FATAL(error.what() << endl);
111 gRandom->SetSeed(seed);
113 if (numberOfEvents <= 0) {
FATAL(
"Invalid number of events " << numberOfEvents << endl); }
127 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
128 emitters.push_back(
JEmitter(i->getID(), i->getUTMPosition() -
detector.getUTMPosition()));
131 if (
detector.empty()) {
FATAL(
"No modules in detector." << endl); }
132 if (emitters.empty()) {
FATAL(
"No emitters in system." << endl); }
137 const double D_m = -
detector.getUTMZ();
162 TH1D h0(
"cpu", NULL, 100, 1.0, 7.0);
163 TH1D
h1(
"chi2/NDF", NULL, 100, 0.0, 5.0);
171 for (
int number_of_events = 0,
count = 0; number_of_events != numberOfEvents; ++number_of_events) {
173 STATUS(
"event: " << setw(10) << number_of_events <<
'\r');
DEBUG(endl);
177 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
182 gRandom->Uniform(-2.0e-2, +2.0e-2));
186 DEBUG(
"Model" << endl << model << endl);
190 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
199 int minimum_number_of_pings = numeric_limits<int>::max();
202 minimum_number_of_pings = min(minimum_number_of_pings,
getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]));
209 const int number_of_pings =
getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]);
211 const double signal = (unify ? (double) minimum_number_of_pings / (
double) number_of_pings : 1.0);
213 for (
int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
217 const double toe_s = ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
221 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
225 const JWaveform& waveform =
getValue(waveforms, emitter->getID(), waveforms[WILDCARD]);
227 const double d_m =
getDistance(module->getPosition(), emitter->getPosition());
228 const double toa_s = toe_s + V.getTime(d_m, emitter->getZ(), module->getZ());
229 const double Q = waveform.
getQ(D_m, d_m);
236 t1_s = gRandom->Gaus(toa_s,
parameters.sigma_s);
243 module->getLocation(),
246 DEBUG(
"hit: " << hit <<
' ' <<
FIXED(7,1) << Q << endl);
260 double chi2 = numeric_limits<double>::max();
266 chi2 = evaluator(result, data.begin(), data.end());
273 chi2 = simplex(data.begin(), data.end()) / simplex.
estimator->getRho(1.0);
274 result = simplex.
value;
281 chi2 = gandalf(data.begin(), data.end()) / gandalf.
estimator->getRho(1.0);
282 result = gandalf.
value;
294 W += hit->getWeight();
297 const int ndf = data.size() - result.
getN();
299 DEBUG(
"Final values" << endl
300 <<
FIXED(9,3) << chi2 <<
'/' << ndf << endl
304 h0.Fill(log10((
double) timer.usec_wall));
305 h1.Fill(chi2 / (
double) (W - result.
getN()));
308 H2[i->first]->Fill((i->second.tx - result.
string [i->first].tx) * 1.0e3,
309 (i->second.ty - result.
string [i->first].ty) * 1.0e3);
313 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.
Q(UTCMax_s-UTCMin_s)-livetime_s
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
int main(int argc, char *argv[])
Template specialisation of fit function of acoustic model based on linear approximation.
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
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.
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
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...
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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.
Data structure for optical module.