Jpp  15.0.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
software/JCompass/JCompass.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 #include <vector>
4 #include <map>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 
11 #include "JDB/JSupport.hh"
12 #include "JDB/JAHRS.hh"
14 #include "JDB/JAHRSToolkit.hh"
15 
18 #include "JSupport/JMeta.hh"
19 
20 #include "JROOT/JManager.hh"
21 #include "JFit/JSimplex.hh"
22 #include "JLang/JComparator.hh"
23 #include "JSystem/JStat.hh"
24 
25 #include "JDetector/JDetector.hh"
29 #include "JDetector/JCompass.hh"
30 
31 #include "JCompass/JHit.hh"
32 #include "JCompass/JModel.hh"
33 #include "JCompass/JEvt.hh"
34 #include "JCompass/JEvtToolkit.hh"
35 #include "JCompass/JSupport.hh"
36 
37 #include "Jeep/JProperties.hh"
38 #include "Jeep/JContainer.hh"
39 #include "Jeep/JPrint.hh"
40 #include "Jeep/JParser.hh"
41 #include "Jeep/JMessage.hh"
42 
43 
44 /**
45  * \file
46  *
47  * Program to calibrate in situ AHRS.
48  * \author mdejong
49  */
50 int main(int argc, char **argv)
51 {
52  using namespace std;
53  using namespace JPP;
54 
55  JMultipleFileScanner_t inputFile;
56  counter_type numberOfEvents;
58  string detectorFile;
59  long long int Tmax_s = 600; // time window data collection [s]
60  int mestimator = EM_LORENTZIAN; // M-estimator fit
61  double sigma_deg = 1.0; // resolution [deg]
62  double stdev = numeric_limits<double>::max(); // number of standard deviations
63  int numberOfOutliers = 0;
64  string ahrsFile;
65  bool overwriteDetector;
66  int debug;
67 
68  try {
69 
70  JProperties properties;
71 
72  properties.insert(gmake_property(Tmax_s));
73  properties.insert(gmake_property(mestimator));
74  properties.insert(gmake_property(sigma_deg));
75  properties.insert(gmake_property(stdev));
76  properties.insert(gmake_property(numberOfOutliers));
77 
78  JParser<> zap("Program to calibrate in situ AHRS.");
79 
80  zap['f'] = make_field(inputFile, "output of JConvertDB -q ahrs");
81  zap['n'] = make_field(numberOfEvents) = JLimit::max();
82  zap['a'] = make_field(detectorFile);
83  zap['@'] = make_field(properties) = JPARSER::initialised();
84  zap['c'] = make_field(ahrsFile, "output of JAHRSCalibration");
85  zap['A'] = make_field(overwriteDetector);
86  zap['o'] = make_field(outputFile) = "compass.root";
87  zap['d'] = make_field(debug) = 2;
88 
89  zap(argc, argv);
90  }
91  catch(const exception &error) {
92  FATAL(error.what() << endl);
93  }
94 
95 
97 
98  try {
99  load(detectorFile, detector);
100  }
101  catch(const JException& error) {
102  FATAL(error);
103  }
104 
106 
107  const JModuleRouter router(detector);
108  const JStringRouter string(detector);
109 
110 
111  map<int, JAverage<JQuaternion3D> > compass; // output compass calibration
112 
113 
114  const JAHRSCalibration_t calibration(ahrsFile.c_str());
115  const JAHRSValidity is_valid;
116 
117 
118  JSimplex<JModel> simplex;
119 
121 
122  simplex.debug = debug;
123 
124  const JChi2 getChi2(mestimator);
125 
126 
127  typedef JManager<int, TH1D> JManager_t;
128 
129  JManager_t H0(new TH1D("%.twist", NULL, 100, 0.0, 5.0));
130  JManager_t H1(new TH1D("%.swing", NULL, 250, 0.0, 2.5));
131  JManager_t HN(new TH1D("%.count", NULL, 100, -0.5, 99.5));
132 
133  TH2D h2("h2", NULL,
135  range.getLength(), range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5);
136 
137  for (Int_t i = 1; i <= h2.GetXaxis()->GetNbins(); ++i) {
138  h2.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
139  }
140  for (Int_t i = 1; i <= h2.GetYaxis()->GetNbins(); ++i) {
141  h2.GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i));
142  }
143 
144  TH2D* h1 = (TH2D*) h2.Clone("h1");
145 
146 
147  outputFile.open();
148 
149  outputFile.put(JMeta(argc, argv));
150 
151  counter_type counter = 0;
152 
153  for (JMultipleFileScanner_t::const_iterator file_name = inputFile.begin(); file_name != inputFile.end(); ++file_name) {
154 
155  STATUS("processing file " << *file_name << endl);
156 
157  map<int, vector<JAHRS> > data; // AHRS data per string
158 
159  for (JMultipleFileScanner<JAHRS> in(*file_name); in.hasNext() && counter != numberOfEvents; ++counter) {
160 
161  const JAHRS* parameters = in.next();
162 
163  if (is_valid(*parameters) && router.hasModule(parameters->DOMID)) {
164  data[router.getModule(parameters->DOMID).getString()].push_back(*parameters);
165  }
166  }
167 
168  for (map<int, vector<JAHRS> >::iterator i = data.begin(); i != data.end(); ++i) {
169 
170  sort(i->second.begin(), i->second.end(), make_comparator(&JAHRS::UNIXTIME));
171 
172  for (vector<JAHRS>::const_iterator p = i->second.begin(); p != i->second.end(); ) {
173 
174  long long int t1 = p->UNIXTIME;
175  long long int t2 = t1;
176 
177  vector<JHit> buffer; // calibrated quaternion data
178 
179  for ( ; p != i->second.end() && p->UNIXTIME < t1 + Tmax_s * 1000; t2 = (p++)->UNIXTIME) {
180 
181  if (calibration.has(p->DOMID)) {
182 
183  const JModule& module = router.getModule(p->DOMID);
184 
185  if (module.getFloor() != 0) {
186 
187  const JCompass compass(*p, calibration.get(p->DOMID));
188 
189  const JQuaternion3D Q = module.getQuaternion() * compass.getQuaternion();
190 
191  buffer.push_back(JHit(p->DOMID, module.getZ(), Q, sigma_deg));
192  }
193  }
194  }
195 
196  if (buffer.size() > JModel::NUMBER_OF_PARAMETERS) {
197 
198  for (vector<JHit>::const_iterator hit = buffer.begin(); hit != buffer.end(); ++hit) {
199 
200  const JLocation& location = router.getModule(hit->getID());
201 
202  h1->Fill((double) string.getIndex(location.getString()), (double) location.getFloor());
203  }
204 
205  JModel result(buffer.begin(), buffer.end()); // prefit
206 
207  vector<JHit>::iterator __end = buffer.end();
208 
209  for (int ns = 0; ns != numberOfOutliers; ++ns) { // outlier removal
210 
211  double xmax = 0.0;
212  vector<JHit>::iterator out = __end;
213 
214  for (vector<JHit>::iterator hit = buffer.begin(); hit != __end; ++hit) {
215 
216  const JQuaternion3D Q1 = result(hit->getZ()); // fitted
217  const JQuaternion3D Q2 = hit->getQuaternion(); // measured
218 
219  const double x = getAngle(Q1, Q2);
220 
221  if (x > xmax) {
222  xmax = x;
223  out = hit;
224  }
225  }
226 
227  if (xmax > stdev * sigma_deg) {
228 
229  const JLocation& location = router.getModule(out->getID());
230 
231  h2.Fill((double) string.getIndex(location.getString()), (double) location.getFloor());
232 
233  if (debug >= debug_t) {
234 
235  const JQuaternion3D Q1 = result(out->getZ()); // fitted
236  const JQuaternion3D Q2 = out->getQuaternion(); // measured
237 
240 
241  cout << "remove " << location << ' '
242  << FIXED(5,2) << getAngle(Q1,Q2) << ' '
243  << FIXED(5,2) << getAngle(q1.twist, q2.twist) << ' '
244  << FIXED(5,2) << getAngle(q1.swing, q2.swing) << endl;
245  }
246 
247  swap(*out, *--__end);
248 
249  result = JModel(buffer.begin(), __end); // refit
250 
251  } else {
252 
253  break;
254  }
255  }
256 
257  simplex.value = result; // start value
258 
259  simplex.step.resize(4);
260 
261  simplex.step[0] = JModel(JQuaternion3X(5.0e-1 * PI / 180.0), JQuaternion3D::getIdentity());
262  simplex.step[1] = JModel(JQuaternion3Y(5.0e-1 * PI / 180.0), JQuaternion3D::getIdentity());
263  simplex.step[2] = JModel(JQuaternion3Z(5.0e-1 * PI / 180.0), JQuaternion3D::getIdentity());
264  simplex.step[3] = JModel(JQuaternion3D::getIdentity(), JQuaternion3Z(5.0e-2 * PI / 180.0));
265 
266  const double chi2 = simplex(getChi2, buffer.begin(), __end);
267  const int ndf = distance(buffer.begin(), __end) * 4 - simplex.step.size();
268 
269  result = simplex.value; // final value
270 
271 
272  outputFile.put(getEvt(JHead(t1, t2, i->first, ndf, chi2), result));
273 
274 
275  for (vector<JHit>::const_iterator hit = buffer.begin(); hit != buffer.end(); ++hit) {
276 
277  const JQuaternion3D Q1 = result(hit->getZ()); // fitted
278  const JQuaternion3D Q2 = hit->getQuaternion(); // measured
279 
282 
283  compass[hit->getID()].put(Q1 * Q2.getConjugate());
284 
285  H0[hit->getID()]->Fill(getAngle(q1.twist, q2.twist));
286  H1[hit->getID()]->Fill(getAngle(q1.swing, q2.swing));
287  }
288 
290 
291  for (vector<JHit>::const_iterator hit = buffer.begin(); hit != __end; ++hit) {
292  count[hit->getID()] += 1;
293  }
294 
295  for (map<int, int>::const_iterator i = count.begin(); i != count.end(); ++i) {
296  HN[i->first]->Fill(i->second);
297  }
298  }
299  }
300  }
301  }
302 
303 
304  h2.Divide(h1);
305 
306  outputFile.put(h2);
307 
308  for (JManager_t* p : { &H0, &H1, &HN }) {
309  for (JManager_t::iterator i = p->begin(); i != p->end(); ++i) {
310  outputFile.put(*(i->second));
311  }
312  }
313 
314  outputFile.close();
315 
316 
317  if (overwriteDetector) {
318 
319  NOTICE("Store calibration data on file " << detectorFile << endl);
320 
321  if (detector.setToLatestVersion()) {
322  NOTICE("Set detector version to " << detector.getVersion() << endl);
323  }
324 
325  detector.comment.add(JMeta(argc, argv));
326 
327  for (map<int, JAverage<JQuaternion3D> >::const_iterator i = compass.begin(); i != compass.end(); ++i) {
328 
329  JModule& module = detector[router.getIndex(i->first)];
330 
331  JQuaternion3D Q(i->second * module.getQuaternion());
332 
333  module.setQuaternion(Q.normalise());
334  }
335 
336  store(detectorFile, detector);
337  }
338 }
static int debug
debug level (default is off).
Definition: JMessage.hh:45
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:1500
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
General exception.
Definition: JException.hh:23
Q(UTCMax_s-UTCMin_s)-livetime_s
debug
Definition: JMessage.hh:29
int main(int argc, char *argv[])
Definition: Main.cc:15
This class represents a rotation around the x-axis.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Definition: JComparator.hh:185
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:68
JQuaternion3D getConjugate() const
Get conjugate of this quaternion.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
do rm f tmp H1
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
const int getIndex(const JObjectID &id) const
Get index of module.
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.
stdev
Utility class to parse parameter values.
Definition: JProperties.hh:496
JEvt getEvt(const JHead &header, const JModel &model)
Get event.
mestimator
*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
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:151
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
Long64_t counter_type
Type definition for counter.
Dynamic ROOT object management.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Data structure for detector geometry and calibration.
Template spacialisation for averaging quaternions.
Definition: JEigen3D.hh:98
Utility class to parse parameter values.
Acoustics hit.
ROOT TTree parameter settings.
const JQuaternion3D & getQuaternion() const
Get quaternion.
Model for fit to acoustics data.
Auxiliary data structure for chi2 evaluation.
JQuaternion3D swing
rotation around perpendicular axis
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
I/O formatting auxiliaries.
Detector file.
Definition: JHead.hh:196
JQuaternion3D twist
rotation around parallel axis
void setQuaternion(const JQuaternion3D &quaternion)
Set quaternion.
Logical location of module.
Definition: JLocation.hh:37
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
return result
Definition: JPolint.hh:727
This class represents a rotation around the z-axis.
bool is_valid(const json &js)
Check validity of JSon data.
Data structure for compass in three dimensions.
Definition: JCompass.hh:49
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
ROOT I/O of application specific meta data.
#define NOTICE(A)
Definition: JMessage.hh:64
static const double PI
Mathematical constants.
int getIndex()
Get index for user I/O manipulation.
Definition: JManip.hh:26
Compass event fit.
ROOT TTree parameter settings.
int debug
debug level
Definition: JSirene.cc:63
General purpose messaging.
Monte Carlo run header.
Definition: JHead.hh:1113
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Direct access to string in detector data structure.
#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.
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
int getString() const
Get string number.
Definition: JLocation.hh:134
Auxiliary base class for list of file names.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
std::vector< int > count
Definition: JAlgorithm.hh:180
Router for mapping of string identifier to index.
General purpose class for object reading from a list of file names.
This class represents a rotation around the y-axis.
Utility class to parse command line options.
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
bool hasModule(const JObjectID &id) const
Has module.
Auxiliary class to map module identifier to AHRS calibration.
Compass event data types.
do set_variable DETECTOR_TXT $WORKDIR detector
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:70
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
Definition: JCanberra.sh:41
double getZ() const
Get z position.
Definition: JVector3D.hh:115
Container I/O.
Auxiliary data structure to check validity of AHRS data.
Definition: JAHRSToolkit.hh:20
then usage $script< detector file >< tripodfile >< stage > input file nInput files correspond to the output of JAcousticsEventBuilder[.sh] nFirst stage eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY eval JPrintDetector a $DETECTOR O CAN source JAcoustics sh $DETECTOR_ID typeset A CONFIGURATION for key in Tmax_s
File status.