Jpp  19.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
software/JCompass/JBallarat.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <algorithm>
5 #include <utility>
6 #include <vector>
7 
9 
10 #include "JDB/JSupport.hh"
11 #include "JDB/JAHRS.hh"
13 #include "JDB/JAHRSToolkit.hh"
14 
17 #include "JSupport/JMeta.hh"
18 
19 #include "JDetector/JDetector.hh"
22 #include "JDetector/JCompass.hh"
23 
25 #include "JGeometry3D/JEigen3D.hh"
26 
27 #include "JTools/JHashMap.hh"
28 #include "JLang/JVectorize.hh"
29 #include "JLang/JComparator.hh"
30 #include "JMath/JLegendre.hh"
31 
32 #include "JCompass/JEvt.hh"
33 #include "JCompass/JEvtToolkit.hh"
34 #include "JCompass/JSupport.hh"
35 #include "JCompass/JNOAA.hh"
36 #include "JCompass/JPolicy.hh"
38 
39 #include "Jeep/JPrint.hh"
40 #include "Jeep/JParser.hh"
41 #include "Jeep/JMessage.hh"
42 
43 namespace {
44 
45  static const int NUMBER_OF_POINTS = 20; //!< number of points for interpolation
46  static const int NUMBER_OF_DEGREES = 1; //!< number of degrees for interpolation
47 
48 }
49 
50 
51 /**
52  * \file
53  *
54  * Auxiliary program to process AHRS data.
55  *
56  * The AHRS calibration as well as the quaternion calibration are applied to the data.\n
57  * Possible outliers are removed according the specified maximal angle (option <tt>-S <angle_deg></tt>.\n
58  * For this, interpolations of the data are made based on Legendre polynomials with
59  * the number of points and the number of degrees\n
60  * set to NUMBER_OF_POINTS and NUMBER_OF_DEGREES, respectively.\n
61  * A policy for missing compass data is applied in which data from neighbouring compasses are used.\n
62  * Finally, data within a given time window can be averaged (option <tt>-T <T_s></tt>.\n)
63  * \author mdejong
64  */
65 int main(int argc, char **argv)
66 {
67  using namespace std;
68  using namespace JPP;
69 
70  JMultipleFileScanner_t inputFile;
71  counter_type numberOfEvents;
73  string detectorFile;
74  string ahrsFile;
75  double angle_deg;
76  double T_s;
77  int debug;
78 
79  try {
80 
81  JParser<> zap("Auxiliary program to process AHRS data.");
82 
83  zap['f'] = make_field(inputFile, "output of JConvertDB -q ahrs");
84  zap['n'] = make_field(numberOfEvents) = JLimit::max();
85  zap['a'] = make_field(detectorFile);
86  zap['c'] = make_field(ahrsFile, "output of JAHRSCalibration");
87  zap['o'] = make_field(outputFile, "ROOT formatted orientations file");
88  zap['S'] = make_field(angle_deg, "maximal angle w.r.t. fit") = 0.0;
89  zap['T'] = make_field(T_s, "time window for averaging") = 0.0;
90  zap['d'] = make_field(debug) = 1;
91 
92  zap(argc, argv);
93  }
94  catch(const exception &error) {
95  FATAL(error.what() << endl);
96  }
97 
98 
100 
101  try {
102  load(detectorFile, detector);
103  }
104  catch(const JException& error) {
105  FATAL(error);
106  }
107 
108  const JModuleRouter router(detector);
109  const JNOAAFunction1D_t& getMagneticDeclination = (isARCADetector(detector) ? static_cast<const JNOAAFunction1D_t&>(getARCAMagneticDeclination) :
111  static_cast<const JNOAAFunction1D_t&>(getZEROMagneticDeclination));
112  const double meridian = (isARCADetector(detector) ? ARCA_MERIDIAN_CONVERGENCE_ANGLE_RAD :
114  0.0);
115 
116  NOTICE("Magnetic declination " << getMagneticDeclination << endl);
117  NOTICE("Meridian angle [rad] " << FIXED(5,3) << meridian << endl);
118 
119  const JAHRSCalibration_t calibration(ahrsFile.c_str());
120  const JAHRSValidity is_valid;
121 
122 
123  typedef pair<double, JQuaternion3D> element_type;
125  typedef JHashMap<int, buffer_type> map_type; // container for temporary storage of data
126 
127  JLegendre<JQuaternion3D, NUMBER_OF_DEGREES> f1; // interpolator for outlier removal
128 
129 
130  outputFile.open();
131 
132  outputFile.put(JMeta(argc, argv));
133 
134  counter_type counter = 0;
135 
136  for (JMultipleFileScanner_t::const_iterator file_name = inputFile.begin(); file_name != inputFile.end(); ++file_name) {
137 
138  STATUS("processing file " << *file_name << "... " << flush);
139 
140  // collect data
141 
142  map_type data;
143 
144  for (JMultipleFileScanner<JAHRS> in(*file_name); in.hasNext() && counter != numberOfEvents; ++counter) {
145 
146  const JAHRS* parameters = in.next();
147  const int id = parameters->DOMID;
148 
149  if (is_valid(*parameters) && calibration.has(id) && router.hasModule(id)) {
150 
151  const JModule& module = router.getModule(id);
152 
153  if (module.getFloor() != 0 && !module.has(COMPASS_DISABLE)) {
154 
155  JCompass compass(*parameters, calibration.get(id));
156 
157  compass.correct(getMagneticDeclination(parameters->UNIXTIME * 1.0e-3), meridian);
158 
159  const JQuaternion3D Q = compass.getQuaternion();
160 
161  data[id].push_back(element_type(parameters->UNIXTIME * 1.0e-3, Q));
162  }
163  }
164  }
165 
166  // remove outliers
167 
168  if (angle_deg > 0.0) {
169 
170  for (map_type::iterator module = data.begin(); module != data.end(); ++module) {
171 
172  if (!module->second.empty()) {
173 
174  sort(module->second.begin(), module->second.end(), make_comparator(&element_type::first));
175 
176  buffer_type::iterator out = module->second.begin(); // output
177  buffer_type::const_iterator in = module->second.begin(); // input
178  buffer_type::const_iterator p = module->second.begin(); // begin interpolation data
179  buffer_type::const_iterator q = module->second.begin(); // end interpolation data
180 
181  for (int i = 0; i != NUMBER_OF_POINTS && q != module->second.end(); ++i, ++q) {}
182 
183  for (int i = 0; i != NUMBER_OF_POINTS/2 && in != q; ++i, ++in) {
184 
185  f1.set(p, in, q);
186 
187  if (getAngle(in->second, f1(in->first)) <= angle_deg) {
188  *out = *in;
189  ++out;
190  }
191  }
192 
193  for (++p; q++ != module->second.end(); ++p, ++in) {
194 
195  f1.set(p, in, q);
196 
197  if (getAngle(in->second, f1(in->first)) <= angle_deg) {
198  *out = *in;
199  ++out;
200  }
201  }
202 
203  for ( ; in != module->second.end(); ++in) {
204 
205  f1.set(p, in, q);
206 
207  if (getAngle(in->second, f1(in->first)) <= angle_deg) {
208  *out = *in;
209  ++out;
210  }
211  }
212 
213  module->second.erase(out, module->second.end()); // remove leftovers
214  }
215  }
216  }
217 
218 
219  // apply policy for missing modules and write output
220 
221  const array_type<map_type::key_type>& buffer = make_array(data.begin(), data.end(), &map_type::value_type::first);
222 
223  const JPolicy policy(router, buffer.begin(), buffer.end());
224 
225  for (JPolicy::const_iterator module = policy.begin(); module != policy.end(); ++module) {
226 
227  buffer_type buffer;
228 
229  for (JPolicy::mapped_type::const_iterator in = module->second.begin(); in != module->second.end(); ++in) {
230  for (map_type::mapped_type::const_iterator i = data[*in].begin(); i != data[*in].end(); ++i) {
231  buffer.push_back(element_type(i->first, i->second));
232  }
233  }
234 
235  // average
236 
237  if (T_s > 0.0) {
238 
239  sort(buffer.begin(), buffer.end(), make_comparator(&element_type::first));
240 
241  buffer_type::iterator out = buffer.begin(); // output
242  buffer_type::const_iterator p = buffer.begin(); // begin averaging data
243  buffer_type::const_iterator q = buffer.begin(); // end averaging data
244 
245  for (p = buffer.begin(); p != buffer.end(); p = q) {
246 
247  for (q = p; q != buffer.end() && q->first - p->first < T_s; ++q) {}
248 
249  if (distance(p,q) > 1) {
250 
252  const array_type<JQuaternion3D>& q1 = make_array(p, q, &element_type::second);
253 
254  out->first = getAverage(t1.begin(), t1.end());
255  out->second = getAverage(q1.begin(), q1.end());
256 
257  } else {
258 
259  *out = *p;
260  }
261 
262  ++out;
263  }
264 
265  buffer.erase(out, buffer.end()); // remove leftovers
266  }
267 
268  for (buffer_type::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
269  outputFile.put(JOrientation(module->first, i->first, getQuaternion(i->second)));
270  }
271  }
272 
273  STATUS(endl);
274  }
275 
276  JMultipleFileScanner<JMeta> io(inputFile);
277 
278  io >> outputFile;
279 
280  outputFile.close();
281 }
bool isORCADetector(const JDetectorHeader &header)
Check if given detector header is compatible with that of ORCA.
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1711
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
General exception.
Definition: JException.hh:24
static JARCAMagneticDeclination getARCAMagneticDeclination
Function object for magnetic declination at ARCA site.
Definition: JNOAA.hh:730
Q(UTCMax_s-UTCMin_s)-livetime_s
void correct(const double declination, const double meridian)
Correct compass for magnetic declination and meridian convergence angle.
Definition: JCompass.hh:274
int main(int argc, char *argv[])
Definition: Main.cc:15
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for a composite optical module.
Definition: JModule.hh:67
bool isARCADetector(const JDetectorHeader &header)
Check if given detector header is compatible with tat of ARCA.
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 keys.
Definition: JHashMap.hh:72
JQuaternion3D getQuaternion(const JQuaternion &Q)
Get quaternion.
General purpose class for hash map of unique elements.
#define STATUS(A)
Definition: JMessage.hh:63
bool is_valid(const json &js)
Check validity of JSon data.
Detector data structure.
Definition: JDetector.hh:89
static JZEROMagneticDeclination getZEROMagneticDeclination
Function object for zero magnetic declination.
Definition: JNOAA.hh:732
Recording of objects on file according a format that follows from the file name extension.
Router for direct addressing of module data in detector data structure.
void set(const double *pA)
Set parameter values.
Definition: JMathlib.hh:1465
std::iterator_traits< T >::value_type getAverage(T __begin, T __end)
Get average.
Definition: JMath.hh:494
*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
Definition: diff-Tuna.sh:38
Long64_t counter_type
Type definition for counter.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
static const double ARCA_MERIDIAN_CONVERGENCE_ANGLE_RAD
ARCA meridian convergence angle [rad].
Data structure for detector geometry and calibration.
static JORCAMagneticDeclination getORCAMagneticDeclination
Function object for magnetic declination at ORCA site.
Definition: JNOAA.hh:731
std::vector< JHitW0 > buffer_type
hits
Definition: JPerth.cc:70
ROOT TTree parameter settings.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
Calibration.
Definition: JHead.hh:328
static const int COMPASS_DISABLE
Enable (disable) use of compass if this status bit is 0 (1);.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
I/O formatting auxiliaries.
Detector file.
Definition: JHead.hh:226
Auxiliary base class for interpolation of magnetic declination data obtained from website of NOAA...
Definition: JNOAA.hh:25
bool has(const int bit) const
Test PMT status.
Definition: JStatus.hh:120
Auxiliary class to define policy for invalid modules.
Definition: JPolicy.hh:31
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition: JVectorize.hh:54
Template definition for function evaluation of Legendre polynome.
Definition: JLegendre.hh:213
Data structure for compass in three dimensions.
Definition: JCompass.hh:49
ROOT I/O of application specific meta data.
#define NOTICE(A)
Definition: JMessage.hh:64
Compass event fit.
ROOT TTree parameter settings.
Meridian convergence angles for different sites.
Auxiliary methods to convert data members or return values of member methods of a set of objects to a...
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Data structure for unit quaternion in three dimensions.
Direct access to module in detector data structure.
Auxiliary base class for list of file names.
Policy for invalid modules.
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.
then fatal The output file must have the wildcard in the e g root fi eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:48
Utility class to parse command line options.
bool hasModule(const JObjectID &id) const
Has module.
Auxiliary class to map module identifier to AHRS calibration.
Orientation of module.
Compass event data types.
std::map< int, range_type > map_type
do set_variable DETECTOR_TXT $WORKDIR detector
static const double ORCA_MERIDIAN_CONVERGENCE_ANGLE_RAD
ORCA meridian convergence angle [rad].
long long int UNIXTIME
[ms]
Definition: JAHRS.hh:26
Auxiliary data structure for return type of make methods.
Definition: JVectorize.hh:26
int debug
debug level
Auxiliary data structure to check validity of AHRS data.
Definition: JAHRSToolkit.hh:20