Jpp  master_rocky-40-g5f0272dcd
the software that should make you happy
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 
12 
13 #include "JDB/JSupport.hh"
14 #include "JDB/JAHRS.hh"
16 #include "JDB/JAHRSToolkit.hh"
17 
20 #include "JSupport/JMeta.hh"
21 
22 #include "JROOT/JManager.hh"
23 #include "JFit/JSimplex.hh"
24 #include "JLang/JComparator.hh"
25 #include "JSystem/JStat.hh"
26 
27 #include "JDetector/JDetector.hh"
31 #include "JDetector/JCompass.hh"
32 
33 #include "JCompass/JHit.hh"
34 #include "JCompass/JModel.hh"
35 #include "JCompass/JEvt.hh"
36 #include "JCompass/JEvtToolkit.hh"
37 #include "JCompass/JSupport.hh"
38 
39 #include "Jeep/JProperties.hh"
40 #include "Jeep/JContainer.hh"
41 #include "Jeep/JPrint.hh"
42 #include "Jeep/JParser.hh"
43 #include "Jeep/JMessage.hh"
44 
45 
46 /**
47  * \file
48  *
49  * Program to calibrate in situ AHRS.
50  * \author mdejong
51  */
52 int main(int argc, char **argv)
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 
107  const floor_range range = getRangeOfFloors(detector);
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,
136  string.size(), -0.5, string.size() - 0.5,
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-1));
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 
291  map<int, int> count;
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 }
Compass event fit.
Compass event data types.
ROOT TTree parameter settings.
Container I/O.
string outputFile
ROOT TTree parameter settings.
Data structure for detector geometry and calibration.
Recording of objects on file according a format that follows from the file name extension.
Dynamic ROOT object management.
General purpose messaging.
#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:69
ROOT I/O of application specific meta data.
Direct access to module in detector data structure.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:72
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
File status.
Direct access to string in detector data structure.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Monte Carlo run header.
Definition: JHead.hh:1236
Data structure for compass in three dimensions.
Definition: JCompass.hh:51
JQuaternion3D getQuaternion() const
Get quaternion.
Definition: JCompass.hh:201
Detector data structure.
Definition: JDetector.hh:96
Logical location of module.
Definition: JLocation.hh:40
int getFloor() const
Get floor number.
Definition: JLocation.hh:146
int getString() const
Get string number.
Definition: JLocation.hh:135
Router for direct addressing of module data in detector data structure.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
bool hasModule(const JObjectID &id) const
Has module.
const int getIndex(const JObjectID &id) const
Get index of module.
Data structure for a composite optical module.
Definition: JModule.hh:75
bool has(const int bit) const
Test PMT status.
Definition: JStatus.hh:120
Utility class to parse parameter values.
Definition: JProperties.hh:501
JModel_t value
Definition: JSimplex.hh:240
std::vector< JModel_t > step
Definition: JSimplex.hh:241
Data structure for unit quaternion in three dimensions.
const JQuaternion3D & getQuaternion() const
Get quaternion.
void setQuaternion(const JQuaternion3D &quaternion)
Set quaternion.
JQuaternion3D getConjugate() const
Get conjugate of this quaternion.
JQuaternion3D & normalise()
Normalise quaternion.
double getZ() const
Get z position.
Definition: JVector3D.hh:115
General exception.
Definition: JException.hh:24
Utility class to parse command line options.
Definition: JParser.hh:1698
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition: JManager.hh:47
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.
T getLength() const
Get length (difference between upper and lower limit).
Definition: JRange.hh:289
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
static const int COMPASS_DISABLE
Enable (disable) use of compass if this status bit is 0 (1);.
const double xmax
Definition: JQuadrature.cc:24
JEvt getEvt(const JHead &header, const JModel &model)
Get event.
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
@ debug_t
debug
Definition: JMessage.hh:29
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:56
@ EM_LORENTZIAN
Definition: JMEstimator.hh:187
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
int getIndex()
Get index for user I/O manipulation.
Definition: JManip.hh:26
static const double PI
Mathematical constants.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Long64_t counter_type
Type definition for counter.
std::vector< size_t > ns
bool is_valid(const json &js)
Check validity of JSon data.
Definition: JSTDTypes.hh:14
int main(int argc, char **argv)
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
Calibration.
Definition: JHead.hh:330
Detector file.
Definition: JHead.hh:227
Acoustics hit.
Model for fit to acoustics data.
Auxiliary data structure for chi2 evaluation.
Auxiliary class to map module identifier to AHRS calibration.
Auxiliary data structure to check validity of AHRS data.
Definition: JAHRSToolkit.hh:20
Router for mapping of string identifier to index.
static int debug
debug level (default is off).
Definition: JMessage.hh:45
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
JQuaternion3D swing
rotation around perpendicular axis
JQuaternion3D twist
rotation around parallel axis
This class represents a rotation around the x-axis.
This class represents a rotation around the y-axis.
This class represents a rotation around the z-axis.
Template spacialisation for averaging quaternions.
Definition: JEigen3D.hh:98
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:68
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:72
Auxiliary base class for list of file names.