66 tripods_container tripods;
67 hydrophones_container hydrophones;
75 JParser<> zap(
"Application to fit position calibration model to acoustic data.");
77 zap[
'f'] =
make_field(inputFile,
"output of JAcousticEventBuilder[.sh]");
79 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
87 zap[
'u'] =
make_field(unify,
"unify weighing of pings");
92 catch(
const exception &error) {
93 FATAL(error.what() << endl);
112 for (JDetector::const_iterator i =
detector.begin(); i !=
detector.end(); ++i) {
113 receivers[i->getID()] = i->getLocation();
116 for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
117 emitters[i->getID()] =
JEmitter(i->getID(),
118 i->getUTMPosition() -
detector.getUTMPosition());
137 simplex.debug =
debug;
138 gandalf.debug =
debug;
141 JKatoomba_t::setOption(
false);
148 TH1D h0(
"chi2/NDF", NULL, 5000, 0.0, 100.0);
149 TH1D
h1(
"h1", NULL, 51, -0.5, 50.5);
159 for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
160 HA->GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
161 HB->GetXaxis()->SetBinLabel(i,
MAKE_CSTRING(geometry.at(i-1).first));
170 while (inputFile.hasNext()) {
172 STATUS(
"input " << setw(6) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
174 const JEvent* evt = inputFile.next();
180 if (oid != evt->
getOID()) {
181 FATAL(
"Invalid detector identifier " << evt->
getOID() <<
" != " << oid << endl);
184 zbuf.push_back(*evt);
188 sort(zbuf.begin(), zbuf.end());
200 for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() +
parameters.Tmax_s; ) {}
209 numberOfPings[i->getID()] += 1;
213 DEBUG(
"Number of pings " << setw(2) << i->first <<
' ' << setw(3) << i->second << endl);
216 int minimum_number_of_pings = numeric_limits<int>::max();
219 minimum_number_of_pings = min(minimum_number_of_pings, i->second);
226 const JEmitter& emitter = emitters[i->getID()];
228 const double signal = (unify ? (double) minimum_number_of_pings / (
double) numberOfPings[i->getID()] : 1.0);
234 for (
JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
235 buffer[hit->getID()].push_back(*hit);
240 if (receivers.has(ps->first) && geometry.hasLocation(receivers[ps->first])) {
244 if (ps->second.begin()->getQ() >=
parameters.Qmin) {
248 receivers[ps->first],
256 DEBUG(
"hit: " << *hit << endl);
263 DEBUG(
"prefit:" << endl
268 HA[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
276 double chi2_old = evaluator(result, data.begin(), data.end());
286 const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->sigma;
298 iter_swap(out, --__end);
300 result = estimator(data.begin(), __end);
302 double chi2_new = evaluator(result, data.begin(), __end);
306 DEBUG(
"Remove outlier " << __end->getLocation() <<
' ' << xmax << endl);
321 DEBUG(
"hit: " << *hit << endl);
324 DEBUG(
"prefit:" << endl
328 HB[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
332 double chi2 = numeric_limits<double>::max();
338 result = estimator(data.begin(), data.end());
339 chi2 = evaluator(result, data.begin(), data.end());
346 chi2 = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
347 result = simplex.value;
354 chi2 = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
355 result = gandalf.value;
365 W += hit->getWeight();
368 const int ndf = data.size() - result.
getN();
370 DEBUG(
"result:" << endl
371 <<
FIXED(9,3) << chi2 <<
'/' << ndf << endl
374 h0.Fill(chi2 / (W - result.
getN()));
378 const JEvt evt =
getEvt(
JHead(oid, data.begin()->getValue(), data.rbegin()->getValue(), ndf, (W - result.
getN()), chi2),
result);
386 const double toe = result.
emitter[
JEKey(i->getID(), i->getCounter())].t1;
388 for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
Utility class to parse command line options.
size_t getN() const
Get number of fit parameters.
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.
Template specialisation of fit function of acoustic model based on linear approximation.
Template specialisation of fit function of acoustic model based on JSimplex minimiser.
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)...
Auxiliary data structure for floating point format specification.
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.
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
const std::string & getOID() const
Get detector identifier.
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
double getQ(const double D_m, const double f_kHz, const double d_m)
Get relative quality for given frequency at given distance.
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Implementation for velocity of sound.
JACOUSTICS::JModel::emitter_type emitter
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
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.
const JLimit & getLimit() const
Get limit.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
do set_variable DETECTOR_TXT $WORKDIR detector
do if[[!-f $ACOUSTICS_WORKDIR/${KEY}.txt]]