45 static const int WILDCARD = -1;
47 static const int NUMBER_OF_PINGS = 11;
48 static const JWaveform EMITTER_WAVEFORM(2.5e6,
59 int main(
int argc,
char **argv)
75 hydrophones_container hydrophones;
80 disable_container disable;
81 waveform_container waveforms;
85 numberOfPings[WILDCARD] = NUMBER_OF_PINGS;
86 waveforms [WILDCARD] = EMITTER_WAVEFORM;
90 JParser<> zap(
"Example application to test fit of model to simulated acoustic data.");
101 zap[
'u'] =
make_field(unify,
"unify weighing of pings");
109 catch(
const exception &error) {
110 FATAL(error.what() << endl);
116 gRandom->SetSeed(seed);
118 if (numberOfEvents <= 0) {
FATAL(
"Invalid number of events " << numberOfEvents << endl); }
132 for (tripods_container::const_iterator i =
tripods.begin(); i !=
tripods.end(); ++i) {
133 emitters.push_back(
JEmitter(i->getID(), i->getUTMPosition() -
detector.getUTMPosition()));
136 if (
detector.empty()) {
FATAL(
"No modules in detector." << endl); }
137 if (emitters.empty()) {
FATAL(
"No emitters in system." << endl); }
142 const double D_m = -
detector.getUTMZ();
167 TH1D h0(
"cpu", NULL, 100, 1.0, 7.0);
168 TH1D
h1(
"chi2/NDF", NULL, 100, 0.0, 5.0);
176 for (
int number_of_events = 0,
count = 0; number_of_events != numberOfEvents; ++number_of_events) {
178 STATUS(
"event: " << setw(10) << number_of_events <<
'\r');
DEBUG(endl);
182 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
187 gRandom->Uniform(-2.0e-2, +2.0e-2));
191 DEBUG(
"Model" << endl << model << endl);
195 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
204 int minimum_number_of_pings = numeric_limits<int>::max();
207 minimum_number_of_pings = min(minimum_number_of_pings,
getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]));
214 const int number_of_pings =
getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]);
216 const double signal = (unify ? (double) minimum_number_of_pings / (
double) number_of_pings : 1.0);
218 for (
int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
222 const double toe_s = ping_counter * 10.0 + gRandom->Uniform(-1.0, +1.0);
226 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
228 if (disable.count(
JTransmission_t(emitter->getID(), module->getID())) == 0) {
232 const JWaveform& waveform =
getValue(waveforms, emitter->getID(), waveforms[WILDCARD]);
234 const double d_m =
getDistance(module->getPosition(), emitter->getPosition());
235 const double toa_s = toe_s +
V.getTime(d_m, emitter->getZ(), module->getZ());
236 const double Q = waveform.
getQ(D_m, d_m);
243 t1_s = gRandom->Gaus(toa_s,
parameters.sigma_s);
250 module->getLocation(),
253 DEBUG(
"hit: " << hit <<
' ' <<
FIXED(7,1) << Q << endl);
268 double chi2 = numeric_limits<double>::max();
274 chi2 = evaluator(result, data.begin(), data.end());
281 chi2 = simplex(data.begin(), data.end()) / simplex.
estimator->getRho(1.0);
282 result = simplex.
value;
289 chi2 = gandalf(data.begin(), data.end()) / gandalf.
estimator->getRho(1.0);
290 result = gandalf.
value;
302 W += hit->getWeight();
305 const int ndf = data.size() - result.
getN();
307 DEBUG(
"Final values" << endl
308 <<
FIXED(9,3) << chi2 <<
'/' << ndf << endl
312 h0.Fill(
log10((
double) timer.usec_wall));
313 h1.Fill(chi2 / (
double) (W - result.
getN()));
316 H2[i->first]->Fill((i->second.tx - result.
string [i->first].tx) * 1.0e3,
317 (i->second.ty - result.
string [i->first].ty) * 1.0e3);
321 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.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
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
set_variable E_E log10(E_{fit}/E_{#mu})"
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 depth dependend 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...
Acoustic transmission identifier.
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.
Data structure for tripod.
then fatal Not enough tripods
Acoustic transmission identifier.
Data structure for optical module.