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();
153 simplex.debug =
debug;
154 gandalf.debug =
debug;
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) {
191 if (geometry.hasLocation(module->getLocation())) {
192 module->set(estimator.detector[module->getString()].getPosition(model.
string[module->getString()], module->getFloor()));
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) {
223 if (geometry.hasLocation(module->getLocation())) {
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);
238 t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
239 toa_s + T_s.getUpperLimit());
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);
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.
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)...
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.
do for((RUN=${RANGE%%-*};$RUN<=${RANGE##*-};RUN+=1))
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
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...
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Custom probability density function of time-of-arrival.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
do set_variable DETECTOR_TXT $WORKDIR detector