Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JFootprint.cc File Reference

Auxiliary application to determine tilt angles of seabed based on measured tilt angles. More...

#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF2.h"
#include "TFitResult.h"
#include "JROOT/JMinimizer.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JTools/JRange.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JSuperEvt.hh"
#include "JAcoustics/JSupport.hh"
#include "JLang/JObjectMultiplexer.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

Auxiliary application to determine tilt angles of seabed based on measured tilt angles.

Author
mdejong

Definition in file JFootprint.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 39 of file JFootprint.cc.

40{
41 using namespace std;
42 using namespace JPP;
43
44 typedef JAbstractHistogram<Double_t> JHistogram_t;
45 typedef JRange <Double_t> JRange_t;
46
47 typedef JTYPELIST<JEvt, JSuperEvt>::typelist typelist;
48
49 JMultipleFileScanner<typelist> inputFile;
50 JLimit_t& numberOfEvents = inputFile.getLimit();
51 string outputFile;
52 JHistogram_t X;
53 JHistogram_t Y;
54 double Z;
55 JRange_t x;
56 JRange_t y;
57 string option;
58 bool writeFits;
59 int debug;
60
61 try {
62
63 JParser<> zap("Auxiliary application to determine tilt angles of seabed based on measured tilt angles.");
64
65 zap['f'] = make_field(inputFile, "input file (output of JKatoomba[.sh]/JFremantle[.sh])");
66 zap['n'] = make_field(numberOfEvents) = JLimit::max();
67 zap['o'] = make_field(outputFile) = "footprint.root";
68 zap['x'] = make_field(X, "histogram x abscissa") = JHistogram_t(1000, -25.0, +25.0);
69 zap['y'] = make_field(Y, "histogram y abscissa") = JHistogram_t(1000, -25.0, +25.0);
70 zap['X'] = make_field(x, "fit x-range centered at top [mrad]") = JRange_t(-2.0, +2.0);
71 zap['Y'] = make_field(y, "fit y-range centered at top [mrad]") = JRange_t(-2.0, +2.0);
72 zap['Z'] = make_field(Z, "detector height (for 2nd order tilt correction)") = 0.0;
73 zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "L";
74 zap['w'] = make_field(writeFits, "write fit function(s) to ROOT file " << "(\"<name>" << _F2 << "\")");
75 zap['d'] = make_field(debug) = 2;
76
77 zap(argc, argv);
78 }
79 catch(const exception &error) {
80 FATAL(error.what() << endl);
81 }
82
83 Z *= 0.5; // half height
84
85 if (option.find('O') == string::npos) { option += "O"; }
86 if (option.find("R") == string::npos) { option += "R"; }
87 if (option.find("S") == string::npos) { option += "S"; }
88 //if (option.find('N') == string::npos) { option += "N"; }
89 if (debug == 0 && option.find('Q') == string::npos) { option += "Q"; }
90
91 TH2D H2("%", NULL,
92 X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(),
93 Y.getNumberOfBins(), Y.getLowerLimit(), Y.getUpperLimit());
94
95 JObjectMultiplexer<typelist, JEvt> in(inputFile);
96
97 for (counter_type counter = 0; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
98
99 STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
100
101 const JEvt* evt = in.next();
102
103 if (!evt->empty()) {
104
105 double Tx = 0.0;
106 double Ty = 0.0;
107
108 for (JEvt::const_iterator i = evt->begin(); i != evt->end(); ++i) {
109
110 const double tx = (i->tx + i->tx2 * Z) * 1.0e3; // [mrad]
111 const double ty = (i->ty + i->ty2 * Z) * 1.0e3; // [mrad]
112
113 Tx += tx;
114 Ty += ty;
115 }
116
117 Tx /= evt->size();
118 Ty /= evt->size();
119
120 H2.Fill(Tx, Ty);
121 }
122 }
123 STATUS(endl);
124
125
126 double x0 = 0.0;
127 double y0 = 0.0;
128 double z0 = 0.0;
129
130 for (Int_t ix = 1; ix <= H2.GetXaxis()->GetNbins(); ++ix) {
131 for (Int_t iy = 1; iy <= H2.GetYaxis()->GetNbins(); ++iy) {
132 if (H2.GetBinContent(ix,iy) > z0) {
133 x0 = H2.GetXaxis()->GetBinCenter(ix);
134 y0 = H2.GetYaxis()->GetBinCenter(iy);
135 z0 = H2.GetBinContent(ix,iy);
136 }
137 }
138 }
139
140 STATUS("top: " << FIXED(7,3) << x0 << ' ' << FIXED(7,3) << y0 << ' ' << FIXED(9,1) << z0 << endl);
141
142
143 TF2 f2("f2", "[0] * exp(-0.5 * (x-[1])*(x-[1]) / ([2]*[2])) * exp(-0.5 * (y-[3])*(y-[3]) / ([4]*[4])) + [5] + [6]*x + [7]*y");
144
145 f2.SetParameter(0, z0 * 0.7);
146 f2.SetParameter(1, x0);
147 f2.SetParameter(2, 0.5);
148 f2.SetParameter(3, y0);
149 f2.SetParameter(4, 0.5);
150 f2.SetParameter(5, 0.0);
151 f2.SetParameter(6, 0.0);
152 f2.SetParameter(7, 0.0);
153
154 f2.SetParError(0, 0.1);
155 f2.SetParError(1, 0.02);
156 f2.SetParError(2, 0.02);
157 f2.SetParError(3, 0.02);
158 f2.SetParError(4, 0.02);
159 f2.SetParError(5, 0.001);
160 f2.SetParError(6, 0.001);
161 f2.SetParError(7, 0.001);
162
163 f2.SetRange(x0 + x.getLowerLimit(), y0 + y.getLowerLimit(), x0 + x.getUpperLimit(), y0 + y.getUpperLimit());
164
165 H2.GetXaxis()->SetRangeUser(x0 + x.getLowerLimit(), x0 + x.getUpperLimit());
166 H2.GetYaxis()->SetRangeUser(y0 + y.getLowerLimit(), y0 + y.getUpperLimit());
167
168 TFitResultPtr result = H2.Fit(&f2, option.c_str(), "same");
169
170 H2.GetXaxis()->SetRangeUser(H2.GetXaxis()->GetXmin(), H2.GetXaxis()->GetXmax());
171 H2.GetYaxis()->SetRangeUser(H2.GetYaxis()->GetXmin(), H2.GetYaxis()->GetXmax());
172
173 if (result.Get() == NULL) {
174 FATAL("Invalid TFitResultPtr " << H2.GetName() << endl);
175 }
176
177 if (debug >= status_t || !result->IsValid()) {
178 cout << "chi2/NDF " << FIXED(7,3) << result->Chi2() << '/' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
179 }
180
181 cout << "tilt: " << FIXED(9,6) << f2.GetParameter(1)*1.0e-3 << ' ' << FIXED(9,6) << f2.GetParameter(3)*1.0e-3 << endl;
182
183
184 TFile out(outputFile.c_str(), "recreate");
185
186 out << H2;
187
188 if (writeFits) {
189
190 TH2D* h2 = (TH2D*) H2.Clone(MAKE_CSTRING(H2.GetName() << _F2));
191
192 h2->Reset();
193
194 for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix) {
195 for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {
196
197 const Double_t x = h2->GetXaxis()->GetBinCenter(ix);
198 const Double_t y = h2->GetYaxis()->GetBinCenter(iy);
199
200 h2->SetBinContent(ix, iy, f2.Eval(x,y));
201 h2->SetBinError (ix, iy, 0.0);
202 }
203 }
204 }
205
206 out.Write();
207 out.Close();
208}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JACOUSTICS::JHistogram_t
JAbstractHistogram< double > JHistogram_t
Type definition for scan along axis.
Definition JBillabong.cc:61
JEEP::status_t
@ status_t
status
Definition JMessage.hh:30
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JTOOLS::result
return result
Definition JPolint.hh:862
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::JRange_t
Type definition of range.
Definition JHead.hh:43
JACOUSTICS::JEvt
Acoustic event fit.
Definition JAcoustics/JEvt.hh:310
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Generated by doxygen 1.10.0