55 int main(
int argc,
char **argv)
69 tripods_container tripods;
70 hydrophones_container hydrophones;
74 disable_container disable;
79 JParser<> zap(
"Application to fit position calibration model to acoustic data.");
81 zap[
'f'] =
make_field(inputFile,
"output of JAcousticEventBuilder[.sh]");
83 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
91 zap[
'u'] =
make_field(unify,
"unify weighing of pings");
97 catch(
const exception &error) {
98 FATAL(error.what() << endl);
118 for (JDetector::const_iterator i =
detector.begin(); i !=
detector.end(); ++i) {
119 receivers[i->getID()] = i->getLocation();
122 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
123 emitters[i->getID()] =
JEmitter(i->getID(),
124 i->getUTMPosition() -
detector.getUTMPosition());
147 JKatoomba_t::setOption(
false);
154 TH1D h0(
"chi2/NDF", NULL, 50000, 0.0, 1000.0);
155 TH1D
h1(
"h1", NULL, 51, -0.5, 50.5);
159 range.getLength(),
range.getLowerLimit() - 0.5,
range.getLowerLimit() + 0.5));
163 range.getLength(),
range.getLowerLimit() - 0.5,
range.getLowerLimit() + 0.5));
165 for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
166 HA->GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
167 HB->GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
183 limit.setLowerLimit(limit.getLowerLimit() - in.
skip(limit.getLowerLimit()));
185 const size_t N = zbuf.size();
187 for ( ; in.
hasNext() && number_of_events != limit; ++number_of_events) {
189 STATUS(
"input " << setw(8) << number_of_events <<
'\r');
DEBUG(endl);
197 if (oid != evt->
getOID()) {
198 FATAL(
"Invalid detector identifier " << evt->
getOID() <<
" != " << oid << endl);
201 zbuf.push_back(*evt);
210 if (q != zbuf.begin() && q != zbuf.end()) {
216 if (p->empty() || q->empty() || *q < *p) {
228 int counter[] = { 0, 0 };
234 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() +
parameters.Tmax_s; ) {}
243 numberOfPings[i->getID()] += 1;
247 DEBUG(
"Number of pings " << setw(2) << i->first <<
' ' << setw(3) << i->second << endl);
250 int minimum_number_of_pings = numeric_limits<int>::max();
253 minimum_number_of_pings = min(minimum_number_of_pings, i->second);
260 const JEmitter& emitter = emitters[i->getID()];
262 const double signal = (unify ? (double) minimum_number_of_pings / (
double) numberOfPings[i->getID()] : 1.0);
268 for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
271 buffer[hit->getID()].push_back(*hit);
277 if (receivers.has(ps->first) && geometry.
hasLocation(receivers[ps->first])) {
281 if (ps->second.begin()->getQ() >=
parameters.Qmin) {
285 receivers[ps->first],
293 DEBUG(
"hit: " << *hit << endl);
300 DEBUG(
"prefit:" << endl
305 HA[hit->getID()]->Fill(geometry.
getIndex(hit->getString()), hit->getFloor(), 1.0);
313 double chi2_old = evaluator(result, data.begin(), data.end());
323 const double x = fabs(hit->getValue() - estimator.
getToA(result, *hit)) / hit->sigma;
335 iter_swap(out, --__end);
337 result = estimator(data.begin(), __end);
339 double chi2_new = evaluator(result, data.begin(), __end);
343 DEBUG(
"Remove outlier " << __end->getLocation() <<
' ' << xmax << endl);
358 DEBUG(
"hit: " << *hit << endl);
361 DEBUG(
"prefit:" << endl
365 HB[hit->getID()]->Fill(geometry.
getIndex(hit->getString()), hit->getFloor(), 1.0);
369 double chi2 = numeric_limits<double>::max();
375 result = estimator(data.begin(), data.end());
376 chi2 = evaluator(result, data.begin(), data.end());
383 chi2 = simplex(data.begin(), data.end()) / simplex.
estimator->getRho(1.0);
384 result = simplex.
value;
391 chi2 = gandalf(data.begin(), data.end()) / gandalf.
estimator->getRho(1.0);
392 result = gandalf.
value;
402 W += hit->getWeight();
405 const int ndf = data.size() - result.
getN();
407 DEBUG(
"result:" << endl
408 <<
FIXED(9,3) << chi2 <<
'/' << ndf << endl
411 h0.Fill(chi2 / (W - result.
getN()));
417 const JEvt evt =
getEvt(
JHead(oid, data.begin()->getValue(), data.rbegin()->getValue(), ndf, (W - result.
getN()), chi2),
result);
425 const double toe = result.
emitter[
JEKey(i->getID(), i->getCounter())].t1;
427 for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
438 WARNING(endl <<
"Event not written: " << chi2 <<
'/' << ndf << endl);
446 STATUS(
"Number of events written / rejected: " << counter[0] <<
" / " << counter[1] << endl);
static int debug
debug level (default is off).
Utility class to parse command line options.
size_t getN() const
Get number of fit parameters.
int main(int argc, char *argv[])
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
static JDetectorMechanics getMechanics
Function object to get string mechanics.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
General purpose class for hash map of unique elements.
ROOT TTree parameter settings.
Recording of objects on file according a format that follows from the file name extension.
Template specialisation of fit function of acoustic model based on linear approximation.
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
JEvt getEvt(const JHead &header, const JModel &model)
Get event.
*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
#define MAKE_CSTRING(A)
Make C-string.
then for HISTOGRAM in h0 h1
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Long64_t counter_type
Type definition for counter.
Dynamic ROOT object management.
Auxiliary data structure for floating point format specification.
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Data structure for detector geometry and calibration.
Data structure for hydrophone.
size_t getNumberOfEmitters(T __begin, T __end)
Get number of emitters.
Model for fit to acoustics data.
Auxiliary class for defining the range of iterations of objects.
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...
virtual bool hasNext() override
Check availability of next element.
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
const std::string & getOID() const
Get detector identifier.
static int debug
debug level
double getQ(const double D_m, const double f_kHz, const double d_m)
Get relative quality for given frequency at given distance.
General purpose messaging.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Implementation for velocity of sound.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
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
virtual const pointer_type & next() override
Get next element.
double getToA(const JModel &model, const JHit< JPDF_t > &hit) const
Get estimated time-of-arrival for given hit.
bool hasLocation(const JLocation &location) const
Check if this detector has given location.
Custom probability density function of time-of-arrival.
const JLimit & getLimit() const
Get limit.
Fit functions of acoustic model.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for tripod.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Acoustic transmission identifier.
virtual skip_type skip(const skip_type ns) override
Skip items.
do if[[!-f $ACOUSTICS_WORKDIR/${KEY}.txt]]
Data structure for optical module.
then usage $script< string identifier >< detectorfile > event file(toashort file)+" "\nNote that the event files and toashort files should be one-to-one related." fi if (( $