Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JRose.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 "TGraph.h"
10
12#include "JROOT/JManager.hh"
13#include "JROOT/JGraph.hh"
14#include "JLang/JTitle.hh"
15
16#include "JCompass/JModel.hh"
17#include "JCompass/JEvt.hh"
19#include "JCompass/JSupport.hh"
20
21#include "Jeep/JProperties.hh"
22#include "Jeep/JParser.hh"
23#include "Jeep/JMessage.hh"
24
25
26namespace {
27
28 /**
29 * Get rotation angle of given quaternion.
30 *
31 * \param Q quaternion
32 * \return angle [rad]
33 */
34 inline double getAngle(const JGEOMETRY3D::JQuaternion3D& Q)
35 {
36 using namespace std;
37 using namespace JPP;
38
39 return atan2(sqrt(Q.getB()*Q.getB() +
40 Q.getC()*Q.getC() +
41 Q.getD()*Q.getD()),
42 Q.getA()) * 2.0;
43 }
44
45
46 /**
47 * Get signed rotation angle of given quaternion.
48 *
49 * \param Q quaternion
50 * \param v direction
51 * \return angle [rad]
52 */
53 inline double getAngle(const JGEOMETRY3D::JQuaternion3D& Q, const JGEOMETRY3D::JVector3D& v)
54 {
55 using namespace std;
56 using namespace JPP;
57
58 return getAngle(Q) * (v.getDot(Q) >= 0.0 ? +1.0 : -1.0);
59 }
60
61
62 /**
63 * Auxiliary data structure for organising TGraph's.
64 */
65 struct map_type :
66 public std::map<int, JROOT::JGraph_t>,
67 public JLANG::JTitle
68 {
69 /**
70 * Constructor.
71 *
72 * \param title title
73 */
74 map_type(const std::string& title) :
75 JTitle(title)
76 {}
77 };
78}
79
80/**
81 * \file
82 *
83 * Example program to plot compass fit results.
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 outputFile;
94 double stdev = numeric_limits<double>::max(); // maximal chi2 / NDF
95 int debug;
96
97 try {
98
99 JProperties properties;
100
101 properties.insert(gmake_property(stdev));
102
103 JParser<> zap("Example program to plot compass fit results.");
104
105 zap['f'] = make_field(inputFile, "output of JCompass");
106 zap['n'] = make_field(numberOfEvents) = JLimit::max();
107 zap['@'] = make_field(properties) = JPARSER::initialised();
108 zap['o'] = make_field(outputFile) = "rose.root";
109 zap['d'] = make_field(debug) = 2;
110
111 zap(argc, argv);
112 }
113 catch(const exception &error) {
114 FATAL(error.what() << endl);
115 }
116
117
118 JManager<int, TH1D> H0(new TH1D("H0[%]", NULL, 100, 0.0, +0.200));
119 JManager<int, TH1D> H1(new TH1D("H1[%]", NULL, 100, 0.0, +0.005));
120 JManager<int, TH1D> HC(new TH1D("HC[%]", NULL, 300,-2.0, +1.0));
121
122 map_type G0("Q0.twist");
123 map_type G1("Q1.twist");
124 map_type GA("Q0.swing");
125 map_type GB("Q0.atan2");
126
127
128 for (JModel previous; inputFile.hasNext(); ) {
129
130 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
131
132 const JEvt* evt = inputFile.next();
133 const JModel model = getModel(*evt);
134
135 if (inputFile.getCounter() > 1) {
136 H0[evt->id]->Fill(getAngle(model.Q0, previous.Q0));
137 H1[evt->id]->Fill(getAngle(model.Q1, previous.Q1));
138 HC[evt->id]->Fill(log10(evt->chi2 / evt->ndf));
139 }
140
141 previous = model;
142
143 if (evt->chi2 / evt->ndf <= stdev) {
144
147
148 const double t1 = 0.5e-3 * (evt->UNIXTimeStart + evt->UNIXTimeStop);
149
150 G0[evt->id].put(t1, getAngle(q0.twist, JVector3Z_t));
151 G1[evt->id].put(t1, getAngle(q1.twist, JVector3Z_t));
152 GA[evt->id].put(t1, getAngle(q0.swing));
153 GB[evt->id].put(t1, atan2(q0.swing.getB(), q0.swing.getC())); // conform JAcoustics::JModel::getAngle
154 }
155 }
156 STATUS(endl);
157
158
159 TFile out(outputFile.c_str(), "recreate");
160
161 for (JManager<int, TH1D>* h1 : { &H0, &H1, &HC }) {
162 for (JManager<int, TH1D>::const_iterator i = h1->begin(); i != h1->end(); ++i) {
163 out << *(i->second);
164 }
165 }
166
167 for (map_type* g1 : { &G0, &G1, &GA, &GB }) {
168 for (map_type::const_iterator i = g1->begin(); i != g1->end(); ++i) {
169 out << JGraph(i->second, MAKE_CSTRING("[" << i->first << "]." << g1->getTitle()));
170 }
171 }
172
173 out.Write();
174 out.Close();
175}
Compass event fit.
Compass event data types.
ROOT TTree parameter settings.
string outputFile
Dynamic ROOT object management.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
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
#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
Double_t g1(const Double_t x)
Function.
Definition JQuantiles.cc:25
Double_t G1(const Double_t x)
Integral of method g1.
Definition JQuantiles.cc:37
int main(int argc, char **argv)
Definition JRose.cc:86
Utility class to parse parameter values.
Data structure for unit quaternion in three dimensions.
double getB() const
Get b value.
double getD() const
Get d value.
double getC() const
Get c value.
double getA() const
Get a value.
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
double getDot(const JVector3D &vector) const
Get dot product.
Definition JVector3D.hh:282
Auxiliary class for title.
Definition JTitle.hh:19
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
General purpose class for object reading from a list of file names.
virtual bool hasNext() override
Check availability of next element.
counter_type getCounter() const
Get counter.
virtual const pointer_type & next() override
Get next element.
JModel getModel(const JEvt &evt)
Get model.
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition JGandalf.hh:57
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::map< int, range_type > map_type
Acoustic event fit.
double UNIXTimeStop
stop time
double ndf
weighed number of degrees of freedom
double UNIXTimeStart
start time
Model for fit to acoustics data.
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
JQuaternion3D swing
rotation around perpendicular axis
JQuaternion3D twist
rotation around parallel axis
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
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