Jpp - 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 fit AHRS data.
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 fit AHRS data.");
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 
105  const JModuleRouter router(detector);
106  const JStringRouter string(detector);
107 
108 
109  map<int, JAverage<JQuaternion3D> > compass; // output compass calibration
110 
111 
112  const JAHRSCalibration_t calibration(ahrsFile.c_str());
113  const JAHRSValidity is_valid;
114 
115 
116  JSimplex<JModel> simplex;
117 
119 
120  simplex.debug = debug;
121 
122  const JChi2 getChi2(mestimator);
123 
124 
125  typedef JManager<int, TH1D> JManager_t;
126 
127  JManager_t H0(new TH1D("%.twist", NULL, 100, 0.0, 5.0));
128  JManager_t H1(new TH1D("%.swing", NULL, 250, 0.0, 2.5));
129  JManager_t HN(new TH1D("%.count", NULL, 100, -0.5, 99.5));
130 
131  TH2D h2("h2", NULL,
134 
135  for (Int_t i = 1; i <= h2.GetXaxis()->GetNbins(); ++i) {
136  h2.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
137  }
138  for (Int_t i = 1; i <= h2.GetYaxis()->GetNbins(); ++i) {
139  h2.GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i));
140  }
141 
142  TH2D* h1 = (TH2D*) h2.Clone("h1");
143 
144 
145  outputFile.open();
146 
147  outputFile.put(JMeta(argc, argv));
148 
149  counter_type counter = 0;
150 
151  for (JMultipleFileScanner_t::const_iterator file_name = inputFile.begin(); file_name != inputFile.end(); ++file_name) {
152 
153  STATUS("processing file " << *file_name << endl);
154 
155  map<int, vector<JAHRS> > data; // AHRS data per string
156 
157  for (JMultipleFileScanner<JAHRS> in(*file_name); in.hasNext() && counter != numberOfEvents; ++counter) {
158 
159  const JAHRS* parameters = in.next();
160 
161  if (is_valid(*parameters) && router.hasModule(parameters->DOMID)) {
162  data[router.getModule(parameters->DOMID).getString()].push_back(*parameters);
163  }
164  }
165 
166  for (map<int, vector<JAHRS> >::iterator i = data.begin(); i != data.end(); ++i) {
167 
168  sort(i->second.begin(), i->second.end(), make_comparator(&JAHRS::UNIXTIME));
169 
170  for (vector<JAHRS>::const_iterator p = i->second.begin(); p != i->second.end(); ) {
171 
172  long long int t1 = p->UNIXTIME;
173  long long int t2 = t1;
174 
175  vector<JHit> buffer; // calibrated quaternion data
176 
177  for ( ; p != i->second.end() && p->UNIXTIME < t1 + Tmax_s * 1000; t2 = (p++)->UNIXTIME) {
178 
179  if (calibration.has(p->DOMID)) {
180 
181  const JModule& module = router.getModule(p->DOMID);
182 
183  if (module.getFloor() != 0) {
184 
185  const JCompass compass(*p, calibration.get(p->DOMID));
186 
187  const JQuaternion3D Q = module.getQuaternion() * compass.getQuaternion();
188 
189  buffer.push_back(JHit(p->DOMID, module.getZ(), Q, sigma_deg));
190  }
191  }
192  }
193 
194  if (buffer.size() > JModel::NUMBER_OF_PARAMETERS) {
195 
196  for (vector<JHit>::const_iterator hit = buffer.begin(); hit != buffer.end(); ++hit) {
197 
198  const JLocation& location = router.getModule(hit->getID());
199 
200  h1->Fill((double) string.getIndex(location.getString()), (double) location.getFloor());
201  }
202 
203  JModel result(buffer.begin(), buffer.end()); // prefit
204 
205  vector<JHit>::iterator __end = buffer.end();
206 
207  for (int ns = 0; ns != numberOfOutliers; ++ns) { // outlier removal
208 
209  double xmax = 0.0;
210  vector<JHit>::iterator out = __end;
211 
212  for (vector<JHit>::iterator hit = buffer.begin(); hit != __end; ++hit) {
213 
214  const JQuaternion3D Q1 = result(hit->getZ()); // fitted
215  const JQuaternion3D Q2 = hit->getQuaternion(); // measured
216 
217  const double x = getAngle(Q1, Q2);
218 
219  if (x > xmax) {
220  xmax = x;
221  out = hit;
222  }
223  }
224 
225  if (xmax > stdev * sigma_deg) {
226 
227  const JLocation& location = router.getModule(out->getID());
228 
229  h2.Fill((double) string.getIndex(location.getString()), (double) location.getFloor());
230 
231  if (debug >= debug_t) {
232 
233  const JQuaternion3D Q1 = result(out->getZ()); // fitted
234  const JQuaternion3D Q2 = out->getQuaternion(); // measured
235 
238 
239  cout << "remove " << location << ' '
240  << FIXED(5,2) << getAngle(Q1,Q2) << ' '
241  << FIXED(5,2) << getAngle(q1.twist, q2.twist) << ' '
242  << FIXED(5,2) << getAngle(q1.swing, q2.swing) << endl;
243  }
244 
245  swap(*out, *--__end);
246 
247  result = JModel(buffer.begin(), __end); // refit
248 
249  } else {
250 
251  break;
252  }
253  }
254 
255  simplex.value = result; // start value
256 
257  simplex.step.resize(4);
258 
259  simplex.step[0] = JModel(JQuaternion3X(5.0e-1 * PI / 180.0), JQuaternion3D::getIdentity());
260  simplex.step[1] = JModel(JQuaternion3Y(5.0e-1 * PI / 180.0), JQuaternion3D::getIdentity());
261  simplex.step[2] = JModel(JQuaternion3Z(5.0e-1 * PI / 180.0), JQuaternion3D::getIdentity());
262  simplex.step[3] = JModel(JQuaternion3D::getIdentity(), JQuaternion3Z(5.0e-2 * PI / 180.0));
263 
264  const double chi2 = simplex(getChi2, buffer.begin(), __end);
265  const int ndf = distance(buffer.begin(), __end) * 4 - simplex.step.size();
266 
267  result = simplex.value; // final value
268 
269 
270  outputFile.put(getEvt(JHead(t1, t2, i->first, ndf, chi2), result));
271 
272 
273  for (vector<JHit>::const_iterator hit = buffer.begin(); hit != buffer.end(); ++hit) {
274 
275  const JQuaternion3D Q1 = result(hit->getZ()); // fitted
276  const JQuaternion3D Q2 = hit->getQuaternion(); // measured
277 
280 
281  compass[hit->getID()].put(Q1 * Q2.getConjugate());
282 
283  H0[hit->getID()]->Fill(getAngle(q1.twist, q2.twist));
284  H1[hit->getID()]->Fill(getAngle(q1.swing, q2.swing));
285  }
286 
288 
289  for (vector<JHit>::const_iterator hit = buffer.begin(); hit != __end; ++hit) {
290  count[hit->getID()] += 1;
291  }
292 
293  for (map<int, int>::const_iterator i = count.begin(); i != count.end(); ++i) {
294  HN[i->first]->Fill(i->second);
295  }
296  }
297  }
298  }
299  }
300 
301 
302  h2.Divide(h1);
303 
304  outputFile.put(h2);
305 
306  for (JManager_t* p : { &H0, &H1, &HN }) {
307  for (JManager_t::iterator i = p->begin(); i != p->end(); ++i) {
308  outputFile.put(*(i->second));
309  }
310  }
311 
312  outputFile.close();
313 
314 
315  if (overwriteDetector) {
316 
317  NOTICE("Store calibration data on file " << detectorFile << endl);
318 
319  if (detector.setToLatestVersion()) {
320  NOTICE("Set detector version to " << detector.getVersion() << endl);
321  }
322 
323  detector.comment.add(JMeta(argc, argv));
324 
325  for (map<int, JAverage<JQuaternion3D> >::const_iterator i = compass.begin(); i != compass.end(); ++i) {
326 
327  JModule& module = detector[router.getIndex(i->first)];
328 
329  JQuaternion3D Q(i->second * module.getQuaternion());
330 
331  module.setQuaternion(Q.normalise());
332  }
333 
334  store(detectorFile, detector);
335  }
336 }
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
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:57
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:80
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.
Utility class to parse parameter values.
Definition: JProperties.hh:496
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
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:69
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
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
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
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.
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:184
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:79
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:38
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
File status.