Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
examples/JDynamics/JBallarat.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <limits>
5#include <map>
6
7#include "TROOT.h"
8#include "TFile.h"
9#include "TH1D.h"
10#include "TGraph.h"
11
12#include "JROOT/JGraph.hh"
13#include "JROOT/JRootToolkit.hh"
14#include "JROOT/JManager.hh"
15
20
22#include "JCompass/JEvt.hh"
23#include "JCompass/JSupport.hh"
24
26
27#include "Jeep/JPrint.hh"
28#include "Jeep/JParser.hh"
29#include "Jeep/JMessage.hh"
30
31
32namespace {
33
36
37 /**
38 * Auxiliary method to get UTC time range of calibration data.
39 *
40 * \param __begin begin of calibration data
41 * \param __end end of calibration data
42 * \return UTC time range
43 */
44 template<class T>
45 static inline JUTCTimeRange getUTCTimeRange(T __begin, T __end)
46 {
48
49 for (T i = __begin; i != __end; ++i) {
50 if (!i->second.empty()) {
51 range.combine(JUTCTimeRange(i->second.getXmin(), i->second.getXmax()));
52 }
53 }
54
55 return range;
56 }
57
58 /**
59 * Get twist angle.
60 *
61 * \param Q quaternion
62 * \return twist angle [rad]
63 */
64 inline double getTwist(const JQuaternion3D& Q)
65 {
66 using namespace JPP;
67
69
70 double phi = (q.twist.getD() > 0.0 ? +2.0 : -2.0) * acos(q.twist.getA());
71
72 while (phi > +PI) { phi -= 2*PI; }
73 while (phi < -PI) { phi += 2*PI; }
74
75 return phi;
76 }
77}
78
79
80/**
81 * \file
82 *
83 * Example program to plot orientation data.
84 * \author mdejong
85 */
86int main(int argc, char **argv)
87{
88 using namespace std;
89 using namespace JPP;
90
92 JLimit_t& numberOfEvents = inputFile.getLimit();
93 string detectorFile;
94 string outputFile;
95 double Tmax_s;
96 int debug;
97
98 try {
99
100 JParser<> zap("Example program to plot orientation data.");
101
102 zap['f'] = make_field(inputFile, "output of JBallarat[.sh]");
103 zap['n'] = make_field(numberOfEvents) = JLimit::max();
104 zap['a'] = make_field(detectorFile);
105 zap['o'] = make_field(outputFile) = "gold.root";
106 zap['T'] = make_field(Tmax_s);
107 zap['d'] = make_field(debug) = 2;
108
109 zap(argc, argv);
110 }
111 catch(const exception &error) {
112 FATAL(error.what() << endl);
113 }
114
115
117
118 try {
119 load(detectorFile, detector);
120 }
121 catch(const JException& error) {
122 FATAL(error);
123 }
124
125 const JModuleRouter router(detector);
126
127 if (!hasDetectorAddressMap(detector.getID())) {
128 FATAL("No detector address map for detector identier " << detector.getID() << endl);
129 }
130
131 const JDetectorBuilder& demo = getDetectorBuilder(detector.getID());
132
133
134 JDynamics dynamics(detector, Tmax_s);
135
136 STATUS("loading input from file(s) " << inputFile << "... " << flush);
137
138 dynamics.load(inputFile);
139
140 STATUS("OK" << endl);
141
142
144
145 for (JDynamics::JOrientation::const_iterator module = dynamics.orientation.begin(); module != dynamics.orientation.end(); ++module) {
146
147 if (router.getModule(module->first).getFloor() != 0) {
148
149 STATUS("Creating graphics for input data " << setw(10) << module->first << flush << '\r');
150
151 JQuaternion3D Q0 = getRotation(demo.getModule(module->first), router.getModule(module->first));
152
153 Q0.conjugate();
154
155 JGraph_t& z0 = Z0[module->first];
156
157 for (JDynamics::JOrientation::function_type::const_iterator i = module->second.begin(); i != module->second.end(); ++i) {
158
159 const JQuaternion3D Q1 = router.getModule(module->first).getQuaternion() * i->getY();
160
161 z0.put(i->getX(), getTwist(Q1 * Q0));
162 }
163 }
164 }
165 STATUS(endl);
166
168
169 for (JDynamics::JOrientation::const_iterator module = dynamics.orientation.begin(); module != dynamics.orientation.end(); ++module) {
170
171 if (router.getModule(module->first).getFloor() != 0 && module->second.size() > 1u) {
172
173 STATUS("Creating graphics for input data " << setw(10) << module->first << flush << '\r');
174
175 TH1D* h1 = H1[module->first] = new TH1D(MAKE_CSTRING("U[" << module->first << "].twist"), NULL, 500, -180.0, +180.0);
176
177 for (JDynamics::JOrientation::function_type::const_iterator i1 = module->second.begin(), i0 = i1++; i1 != module->second.end(); ++i0, ++i1) {
178
179 JQuaternion3D Q0 = i0->getY();
180 JQuaternion3D Q1 = i1->getY();
181
182 h1->Fill(getTwist(Q1 * Q0.conjugate()) * 180.0 / PI);
183 }
184 }
185 }
186 STATUS(endl);
187
189
190 for (JDynamics::JOrientation::const_iterator module = dynamics.orientation.begin(); module != dynamics.orientation.end(); ++module) {
191
192 STATUS("Creating graphics for compiled data " << setw(10) << module->first << flush << '\r');
193
194 if (router.getModule(module->first).getFloor() != 0 && !module->second.empty()) {
195
196 JQuaternion3D Q0 = getRotation(demo.getModule(module->first), router.getModule(module->first));
197
198 Q0.conjugate();
199
200 const Double_t xmin = module->second.getXmin();
201 const Double_t xmax = module->second.getXmax();
202 const Int_t nx = (Int_t) ((xmax - xmin) / Tmax_s);
203
204 TH1D* h0 = H0[module->first] = new TH1D(MAKE_CSTRING("H[" << module->first << "].twist"), NULL, nx, xmin, xmax);
205
206 for (Int_t i = 1; i <= h0->GetXaxis()->GetNbins(); ++i) {
207
208 const Double_t x = h0->GetXaxis()->GetBinCenter(i);
209 const JQuaternion3D Q = module->second(x) * Q0;
210
211 h0->SetBinContent(i, getTwist(Q));
212 }
213 }
214 }
215 STATUS(endl);
216
217 const Double_t xmin = getUTCTimeRange(dynamics.orientation.begin(), dynamics.orientation.end()).getLowerLimit();
218 const Double_t xmax = getUTCTimeRange(dynamics.orientation.begin(), dynamics.orientation.end()).getUpperLimit();
219 const Int_t nx = (Int_t) ((xmax - xmin) / Tmax_s);
220
221 JManager<int, TH1D> X0(new TH1D("X[%].twist", NULL, nx, xmin, xmax));
222
223 for (Int_t i = 1; i <= X0->GetXaxis()->GetNbins(); ++i) {
224
225 STATUS("Creating graphics for detector data " << FIXED(5,1) << (double) (i * 100) / (double) X0->GetXaxis()->GetNbins() << "%" << flush << '\r');
226
227 dynamics.update(X0->GetXaxis()->GetBinCenter(i));
228
229 for (JDetector::const_iterator module = dynamics.begin(); module != dynamics.end(); ++module) {
230
231 if (module->getFloor() != 0) {
232
233 TH1D* x0 = X0[module->getID()];
234
235 const JQuaternion3D Q = getRotation(router.getModule(module->getID()), *module);
236
237 x0->SetBinContent(i, getTwist(Q));
238 }
239 }
240 }
241 STATUS(endl);
242
243
244 TFile out(outputFile.c_str(), "recreate");
245
246 for (map<int, TH1D*>::const_iterator i = H0.begin(); i != H0.end(); ++i) {
247 out << *i->second;
248 }
249
250 for (map<int, TH1D*>::const_iterator i = H1.begin(); i != H1.end(); ++i) {
251 out << *i->second;
252 }
253
254 out << X0;
255
256 for (map<int, JGraph_t>::const_iterator i = Z0.begin(); i != Z0.end(); ++i) {
257 out << JGraph(i->second, MAKE_CSTRING("G[" << i->first << "].twist"));
258 }
259
260 out.Write();
261 out.Close();
262}
Compass event data types.
ROOT TTree parameter settings.
string outputFile
Detector support kit.
Data structure for detector geometry and calibration.
Dynamic detector calibration.
Dynamic ROOT object management.
General purpose messaging.
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
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
Detector data structure.
Definition JDetector.hh:96
int getFloor() const
Get floor number.
Definition JLocation.hh:146
Router for direct addressing of module data in detector data structure.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for unit quaternion in three dimensions.
const JQuaternion3D & getQuaternion() const
Get quaternion.
JQuaternion3D & conjugate()
Conjugate quaternion.
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
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
Definition JManager.hh:304
General purpose class for object reading from a list of file names.
collection_type::const_iterator const_iterator
Definition JPolfit.hh:184
static JRange< T, JComparator_t > DEFAULT_RANGE()
Default range.
Definition JRange.hh:555
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
range_type & combine(const range_type &range)
Combine ranges.
Definition JRange.hh:432
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
UTC time range [s].
int main(int argc, char **argv)
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
JDetectorBuilder & getDetectorBuilder()
Get detector builder.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
bool hasDetectorAddressMap(const int id)
Check if detector address map is available.
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
static const double PI
Mathematical constants.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
JDAQUTCTimeRange getUTCTimeRange()
Get UTC time range.
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Detector file.
Definition JHead.hh:227
Auxiliary interface for building detector.
const JModule & getModule(const int id=-1, const JLocation &location=JLocation()) const
Get module.
const_iterator begin() const
begin of calibration data
Definition JDynamics.hh:210
data_type::const_iterator const_iterator
Definition JDynamics.hh:146
const_iterator end() const
end of calibration data
Definition JDynamics.hh:211
Dynamic detector calibration.
Definition JDynamics.hh:86
JOrientation orientation
orientation calibration
Definition JDynamics.hh:612
bool update(const double t1_s)
Get detector calibrated at given time.
Definition JDynamics.hh:545
void load(JObjectIterator_t &input)
Load calibration data.
Definition JDynamics.hh:530
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
Auxiliary data structure to build TGraph.
Definition JGraph.hh:44
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128