Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
software/JCompass/JCompass.cc File Reference

Program to fit AHRS data. More...

#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#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 "JROOT/JManager.hh"
#include "JFit/JSimplex.hh"
#include "JLang/JComparator.hh"
#include "JSystem/JStat.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JStringRouter.hh"
#include "JDetector/JCompass.hh"
#include "JCompass/JHit.hh"
#include "JCompass/JModel.hh"
#include "JCompass/JEvt.hh"
#include "JCompass/JEvtToolkit.hh"
#include "JCompass/JSupport.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JContainer.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

Program to fit AHRS data.

Author
mdejong

Definition in file software/JCompass/JCompass.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 50 of file software/JCompass/JCompass.cc.

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
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
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
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.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
Template spacialisation for averaging quaternions.
Definition: JEigen3D.hh:98
Acoustics hit.
const JQuaternion3D & getQuaternion() const
Get quaternion.
Model for fit to acoustics data.
Auxiliary data structure for chi2 evaluation.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
Detector file.
Definition: JHead.hh:196
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.
#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
int debug
debug level
Definition: JSirene.cc:63
Monte Carlo run header.
Definition: JHead.hh:1113
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
#define FATAL(A)
Definition: JMessage.hh:67
Data structure for unit quaternion in three dimensions.
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.
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
Auxiliary class to map module identifier to AHRS calibration.
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
Auxiliary data structure to check validity of AHRS data.
Definition: JAHRSToolkit.hh:20