Jpp  test_elongated_shower_pde
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 calibrate in situ AHRS. More...

#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#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 "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 calibrate in situ AHRS.

Author
mdejong

Definition in file software/JCompass/JCompass.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

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

53 {
54  using namespace std;
55  using namespace JPP;
56 
57  JMultipleFileScanner_t inputFile;
58  counter_type numberOfEvents;
60  string detectorFile;
61  long long int Tmax_s = 600; // time window data collection [s]
62  int mestimator = EM_LORENTZIAN; // M-estimator fit
63  double sigma_deg = 1.0; // resolution [deg]
64  double stdev = numeric_limits<double>::max(); // number of standard deviations
65  int numberOfOutliers = 0;
66  string ahrsFile;
67  bool overwriteDetector;
68  int debug;
69 
70  try {
71 
72  JProperties properties;
73 
74  properties.insert(gmake_property(Tmax_s));
75  properties.insert(gmake_property(mestimator));
76  properties.insert(gmake_property(sigma_deg));
77  properties.insert(gmake_property(stdev));
78  properties.insert(gmake_property(numberOfOutliers));
79 
80  JParser<> zap("Program to calibrate in situ AHRS.");
81 
82  zap['f'] = make_field(inputFile, "output of JConvertDB -q ahrs");
83  zap['n'] = make_field(numberOfEvents) = JLimit::max();
84  zap['a'] = make_field(detectorFile);
85  zap['@'] = make_field(properties) = JPARSER::initialised();
86  zap['c'] = make_field(ahrsFile, "output of JAHRSCalibration");
87  zap['A'] = make_field(overwriteDetector);
88  zap['o'] = make_field(outputFile) = "compass.root";
89  zap['d'] = make_field(debug) = 2;
90 
91  zap(argc, argv);
92  }
93  catch(const exception &error) {
94  FATAL(error.what() << endl);
95  }
96 
97 
99 
100  try {
101  load(detectorFile, detector);
102  }
103  catch(const JException& error) {
104  FATAL(error);
105  }
106 
108 
109  const JModuleRouter router(detector);
110  const JStringRouter string(detector);
111 
112 
113  map<int, JAverage<JQuaternion3D> > compass; // output compass calibration
114 
115 
116  const JAHRSCalibration_t calibration(ahrsFile.c_str());
117  const JAHRSValidity is_valid;
118 
119 
120  JSimplex<JModel> simplex;
121 
123 
124  simplex.debug = debug;
125 
126  const JChi2 getChi2(mestimator);
127 
128 
129  typedef JManager<int, TH1D> JManager_t;
130 
131  JManager_t H0(new TH1D("%.twist", NULL, 100, 0.0, 5.0));
132  JManager_t H1(new TH1D("%.swing", NULL, 250, 0.0, 2.5));
133  JManager_t HN(new TH1D("%.count", NULL, 100, -0.5, 99.5));
134 
135  TH2D h2("h2", NULL,
137  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5);
138 
139  for (Int_t i = 1; i <= h2.GetXaxis()->GetNbins(); ++i) {
140  h2.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(string.at(i-1)));
141  }
142  for (Int_t i = 1; i <= h2.GetYaxis()->GetNbins(); ++i) {
143  h2.GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i));
144  }
145 
146  TH2D* h1 = (TH2D*) h2.Clone("h1");
147 
148 
149  outputFile.open();
150 
151  outputFile.put(JMeta(argc, argv));
152 
153  counter_type counter = 0;
154 
155  for (JMultipleFileScanner_t::const_iterator file_name = inputFile.begin(); file_name != inputFile.end(); ++file_name) {
156 
157  STATUS("processing file " << *file_name << endl);
158 
159  map<int, vector<JAHRS> > data; // AHRS data per string
160 
161  for (JMultipleFileScanner<JAHRS> in(*file_name); in.hasNext() && counter != numberOfEvents; ++counter) {
162 
163  const JAHRS* parameters = in.next();
164 
165  if (is_valid(*parameters) && router.hasModule(parameters->DOMID)) {
166  data[router.getModule(parameters->DOMID).getString()].push_back(*parameters);
167  }
168  }
169 
170  for (map<int, vector<JAHRS> >::iterator i = data.begin(); i != data.end(); ++i) {
171 
172  sort(i->second.begin(), i->second.end(), make_comparator(&JAHRS::UNIXTIME));
173 
174  for (vector<JAHRS>::const_iterator p = i->second.begin(); p != i->second.end(); ) {
175 
176  long long int t1 = p->UNIXTIME;
177  long long int t2 = t1;
178 
179  vector<JHit> buffer; // calibrated quaternion data
180 
181  for ( ; p != i->second.end() && p->UNIXTIME < t1 + Tmax_s * 1000; t2 = (p++)->UNIXTIME) {
182 
183  if (calibration.has(p->DOMID)) {
184 
185  const JModule& module = router.getModule(p->DOMID);
186 
187  if (module.getFloor() != 0 && !module.has(COMPASS_DISABLE)) {
188 
189  const JCompass compass(*p, calibration.get(p->DOMID));
190 
191  const JQuaternion3D Q = module.getQuaternion() * compass.getQuaternion();
192 
193  buffer.push_back(JHit(p->DOMID, module.getZ(), Q, sigma_deg));
194  }
195  }
196  }
197 
198  if (buffer.size() > JModel::NUMBER_OF_PARAMETERS) {
199 
200  for (vector<JHit>::const_iterator hit = buffer.begin(); hit != buffer.end(); ++hit) {
201 
202  const JLocation& location = router.getModule(hit->getID());
203 
204  h1->Fill((double) string.getIndex(location.getString()), (double) location.getFloor());
205  }
206 
207  JModel result(buffer.begin(), buffer.end()); // prefit
208 
209  vector<JHit>::iterator __end = buffer.end();
210 
211  for (int ns = 0; ns != numberOfOutliers; ++ns) { // outlier removal
212 
213  double xmax = 0.0;
214  vector<JHit>::iterator out = __end;
215 
216  for (vector<JHit>::iterator hit = buffer.begin(); hit != __end; ++hit) {
217 
218  const JQuaternion3D Q1 = result(hit->getZ()); // fitted
219  const JQuaternion3D Q2 = hit->getQuaternion(); // measured
220 
221  const double x = getAngle(Q1, Q2);
222 
223  if (x > xmax) {
224  xmax = x;
225  out = hit;
226  }
227  }
228 
229  if (xmax > stdev * sigma_deg) {
230 
231  const JLocation& location = router.getModule(out->getID());
232 
233  h2.Fill((double) string.getIndex(location.getString()), (double) location.getFloor());
234 
235  if (debug >= debug_t) {
236 
237  const JQuaternion3D Q1 = result(out->getZ()); // fitted
238  const JQuaternion3D Q2 = out->getQuaternion(); // measured
239 
242 
243  cout << "remove " << location << ' '
244  << FIXED(5,2) << getAngle(Q1,Q2) << ' '
245  << FIXED(5,2) << getAngle(q1.twist, q2.twist) << ' '
246  << FIXED(5,2) << getAngle(q1.swing, q2.swing) << endl;
247  }
248 
249  swap(*out, *--__end);
250 
251  result = JModel(buffer.begin(), __end); // refit
252 
253  } else {
254 
255  break;
256  }
257  }
258 
259  simplex.value = result; // start value
260 
261  simplex.step.resize(4);
262 
263  simplex.step[0] = JModel(JQuaternion3X(5.0e-1 * PI / 180.0), JQuaternion3D::getIdentity());
264  simplex.step[1] = JModel(JQuaternion3Y(5.0e-1 * PI / 180.0), JQuaternion3D::getIdentity());
265  simplex.step[2] = JModel(JQuaternion3Z(5.0e-1 * PI / 180.0), JQuaternion3D::getIdentity());
266  simplex.step[3] = JModel(JQuaternion3D::getIdentity(), JQuaternion3Z(5.0e-2 * PI / 180.0));
267 
268  const double chi2 = simplex(getChi2, buffer.begin(), __end);
269  const int ndf = distance(buffer.begin(), __end) * 4 - simplex.step.size();
270 
271  result = simplex.value; // final value
272 
273 
274  outputFile.put(getEvt(JHead(t1, t2, i->first, ndf, chi2), result));
275 
276 
277  for (vector<JHit>::const_iterator hit = buffer.begin(); hit != buffer.end(); ++hit) {
278 
279  const JQuaternion3D Q1 = result(hit->getZ()); // fitted
280  const JQuaternion3D Q2 = hit->getQuaternion(); // measured
281 
284 
285  compass[hit->getID()].put(Q1 * Q2.getConjugate());
286 
287  H0[hit->getID()]->Fill(getAngle(q1.twist, q2.twist));
288  H1[hit->getID()]->Fill(getAngle(q1.swing, q2.swing));
289  }
290 
292 
293  for (vector<JHit>::const_iterator hit = buffer.begin(); hit != __end; ++hit) {
294  count[hit->getID()] += 1;
295  }
296 
297  for (map<int, int>::const_iterator i = count.begin(); i != count.end(); ++i) {
298  HN[i->first]->Fill(i->second);
299  }
300  }
301  }
302  }
303  }
304 
305 
306  h2.Divide(h1);
307 
308  outputFile.put(h2);
309 
310  for (JManager_t* p : { &H0, &H1, &HN }) {
311  for (JManager_t::iterator i = p->begin(); i != p->end(); ++i) {
312  outputFile.put(*(i->second));
313  }
314  }
315 
316  outputFile.close();
317 
318 
319  if (overwriteDetector) {
320 
321  NOTICE("Store calibration data on file " << detectorFile << endl);
322 
323  if (detector.setToLatestVersion()) {
324  NOTICE("Set detector version to " << detector.getVersion() << endl);
325  }
326 
327  detector.comment.add(JMeta(argc, argv));
328 
329  for (map<int, JAverage<JQuaternion3D> >::const_iterator i = compass.begin(); i != compass.end(); ++i) {
330 
331  JModule& module = detector[router.getIndex(i->first)];
332 
333  JQuaternion3D Q(i->second * module.getQuaternion());
334 
335  module.setQuaternion(Q.normalise());
336  }
337 
338  store(detectorFile, detector);
339  }
340 }
JDetector detector
Definition: JRunAnalyzer.hh:23
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
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: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
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
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.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Template spacialisation for averaging quaternions.
Definition: JEigen3D.hh:98
Acoustics hit.
const JQuaternion3D & getQuaternion() const
Get quaternion.
static const int COMPASS_DISABLE
Enable (disable) use of compass if this status bit is 0 (1);.
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:224
void setQuaternion(const JQuaternion3D &quaternion)
Set quaternion.
Logical location of module.
Definition: JLocation.hh:37
bool has(const int bit) const
Test PMT status.
Definition: JStatus.hh:120
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
return result
Definition: JPolint.hh:743
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
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:68
Monte Carlo run header.
Definition: JHead.hh:1164
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.
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.
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
Auxiliary class to map module identifier to AHRS calibration.
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:56
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:42
double getZ() const
Get z position.
Definition: JVector3D.hh:115
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