Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JBallarat.cc File Reference

Auxiliary program to process AHRS data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <utility>
#include <vector>
#include "km3net-dataformat/definitions/module_status.hh"
#include "JDB/JSupport.hh"
#include "JDB/JAHRS.hh"
#include "JDB/JAHRSCalibration_t.hh"
#include "JDB/JAHRSToolkit.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JCompass.hh"
#include "JGeometry3D/JQuaternion3D.hh"
#include "JGeometry3D/JEigen3D.hh"
#include "JTools/JHashMap.hh"
#include "JLang/JVectorize.hh"
#include "JLang/JComparator.hh"
#include "JMath/JLegendre.hh"
#include "JCompass/JEvt.hh"
#include "JCompass/JEvtToolkit.hh"
#include "JCompass/JSupport.hh"
#include "JCompass/JNOAA.hh"
#include "JCompass/JPolicy.hh"
#include "JCompass/JCompassSupportkit.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Auxiliary program to process AHRS data.

The AHRS calibration as well as the quaternion calibration are applied to the data.
Possible outliers are removed according the specified maximal angle (option -S <angle_deg>.
For this, interpolations of the data are made based on Legendre polynomials with the number of points and the number of degrees
set to NUMBER_OF_POINTS and NUMBER_OF_DEGREES, respectively.
A policy for missing compass data is applied in which data from neighbouring compasses are used.
Finally, data within a given time window can be averaged (option -T <T_s>.
)

Author
mdejong

Definition in file software/JCompass/JBallarat.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 65 of file software/JCompass/JBallarat.cc.

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));
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());
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
251 const array_type<double>& t1 = make_array(p, q, &element_type::first);
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}
string outputFile
#define STATUS(A)
Definition JMessage.hh:63
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Data structure for compass in three dimensions.
Definition JCompass.hh:51
Detector data structure.
Definition JDetector.hh:96
int getFloor() const
Get floor number.
Definition JLocation.hh:146
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
Definition JModule.hh:75
bool has(const int bit) const
Test PMT status.
Definition JStatus.hh:120
Data structure for unit quaternion in three dimensions.
const JQuaternion3D & getQuaternion() const
Get quaternion.
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Object writing to file.
General purpose class for object reading from a list of file names.
virtual bool hasNext() override
Check availability of next element.
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.
JQuaternion3D getQuaternion(const JQuaternion &Q)
Get quaternion.
static JARCAMagneticDeclination getARCAMagneticDeclination
Function object for magnetic declination at ARCA site.
Definition JNOAA.hh:850
static JORCAMagneticDeclination getORCAMagneticDeclination
Function object for magnetic declination at ORCA site.
Definition JNOAA.hh:851
static const double ORCA_MERIDIAN_CONVERGENCE_ANGLE_RAD
ORCA meridian convergence angle [rad].
static const double ARCA_MERIDIAN_CONVERGENCE_ANGLE_RAD
ARCA meridian convergence angle [rad].
static JZEROMagneticDeclination getZEROMagneticDeclination
Function object for zero magnetic declination.
Definition JNOAA.hh:852
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
bool isORCADetector(const JDetectorHeader &header)
Check if given detector header is compatible with that of ORCA.
bool isARCADetector(const JDetectorHeader &header)
Check if given detector header is compatible with tat of ARCA.
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition JVectorize.hh:54
std::iterator_traits< T >::value_type getAverage(T __begin, T __end)
Get average.
Definition JMath.hh:494
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::vector< JHitW0 > buffer_type
hits
Definition JPerth.cc:70
Long64_t counter_type
Type definition for counter.
bool is_valid(const json &js)
Check validity of JSon data.
std::map< int, range_type > map_type
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Calibration.
Definition JHead.hh:330
Detector file.
Definition JHead.hh:227
Auxiliary base class for interpolation of magnetic declination data obtained from website of NOAA.
Definition JNOAA.hh:28
Orientation of module.
Auxiliary class to define policy for invalid modules.
Definition JPolicy.hh:34
Auxiliary class to map module identifier to AHRS calibration.
Auxiliary data structure to check validity of AHRS data.
long long int UNIXTIME
[ms]
Definition JAHRS.hh:26
Auxiliary data structure for return type of make methods.
Definition JVectorize.hh:28
Template definition for function evaluation of Legendre polynome.
Definition JLegendre.hh:269
void set(const double *pA)
Set parameter values.
Definition JMathlib.hh:1465
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
Auxiliary base class for list of file names.
General purpose class for hash map of unique keys.
Definition JHashMap.hh:75