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

Program to fit histogram in output of JMuonCompass.cc. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TFitResult.h"
#include "TError.h"
#include "JROOT/JMinimizer.hh"
#include "JSupport/JMeta.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JStringRouter.hh"
#include "JReconstruction/JMuonCompass.hh"
#include "Jeep/JProperties.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 fit histogram in output of JMuonCompass.cc.

Author
mdejong

Definition in file examples/JReconstruction/JMuonCompass.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

< minimal abscissa with respect to minimum [rad]

< maximal abscissa with respect to minimum [rad]

< scaling factor for error

< minimal difference in likelihood

Definition at line 36 of file examples/JReconstruction/JMuonCompass.cc.

37{
38 using namespace std;
39 using namespace JPP;
40
41 struct {
42 double XMIN = -0.2*PI; //!< minimal abscissa with respect to minimum [rad]
43 double XMAX = +0.2*PI; //!< maximal abscissa with respect to minimum [rad]
44 double FACTOR = 1.0; //!< scaling factor for error
45 double YMIN = 10.0; //!< minimal difference in likelihood
46 } parameters;
47
48 vector<string> inputFile;
49 string outputFile;
50 string detectorFile;
51 JCounter overwriteDetector;
52 string option;
53 int debug;
54
55 try {
56
57 JProperties properties;
58
59 properties.insert(gmake_property(parameters.XMIN));
60 properties.insert(gmake_property(parameters.XMAX));
61 properties.insert(gmake_property(parameters.FACTOR));
62 properties.insert(gmake_property(parameters.YMIN));
63
64 JParser<> zap("Program to fit histogram in output of JMuonCompass.cc.");
65
66 zap['@'] = make_field(properties) = JPARSER::initialised();
67 zap['f'] = make_field(inputFile, "input files (output from JMuonCompass).");
68 zap['o'] = make_field(outputFile, "output file.") = "hobbit.root";
69 zap['a'] = make_field(detectorFile, "detector file.");
70 zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with corrected orientations;"\
71 "\nUse option -AA to also correct compass calibrations.");
72 zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
73 zap['d'] = make_field(debug) = 1;
74
75 zap(argc, argv);
76 }
77 catch(const exception &error) {
78 FATAL(error.what() << endl);
79 }
80
81
82 gErrorIgnoreLevel = kFatal;
83
84
85 JDetector detector;
86
87 try {
88 load(detectorFile, detector);
89 }
90 catch(const JException& error) {
91 FATAL(error);
92 }
93
94 const JModuleRouter router(detector);
95
96 if (option.find('S') == string::npos) { option += 'S'; }
97 if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
98
99 TH2D* h2 = NULL;
100
101 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
102
103 NOTICE("Processing " << *i << endl);
104
105 TFile in(i->c_str(), "exist");
106
107 if (!in.IsOpen()) {
108 FATAL("File " << *i << " not opened." << endl);
109 }
110
111 TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
112
113 if (p == NULL) {
114 FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
115 }
116
117 if (h2 == NULL)
118 h2 = (TH2D*) p->Clone();
119 else
120 h2->Add(p);
121
122 h2->SetDirectory(0);
123
124 in.Close();
125 }
126
127 if (h2 == NULL) {
128 FATAL("No histogram <" << h2_t << ">." << endl);
129 }
130
131 // auxiliary data structure for fit result
132
133 struct result_type {
134
135 result_type() :
136 value(0.0),
137 error(0.0)
138 {}
139
140 result_type(double value,
141 double error) :
142 value(value),
143 error(error)
144 {}
145
146 double value;
147 double error;
148 };
149
150 typedef map<int, result_type> map_type; // index -> fit result
151
152 map_type zmap;
153
154
155 // histogram fit results
156
157 const TAxis* x_axis = h2->GetXaxis();
158 const TAxis* y_axis = h2->GetYaxis();
159
160 TH1D h0("h0", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
161
162 const floor_range floors = getRangeOfFloors(detector);
163 const JStringRouter strings(detector);
164
165 TH2D H2d("detector", NULL,
166 strings.size(), -0.5, strings.size() - 0.5,
167 floors.getUpperLimit(), 1 - 0.5, floors.getUpperLimit() + 0.5);
168
169 for (Int_t ix = 1; ix <= H2d.GetXaxis()->GetNbins(); ++ix) {
170 H2d.GetXaxis()->SetBinLabel(ix, MAKE_CSTRING(strings.at(ix-1)));
171 }
172 for (Int_t iy = 1; iy <= H2d.GetYaxis()->GetNbins(); ++iy) {
173 H2d.GetYaxis()->SetBinLabel(iy, MAKE_CSTRING(iy));
174 }
175
176 TH2D* H2c = (TH2D*) H2d.Clone("chi2");
177 TH2D* H2q = (TH2D*) H2d.Clone("status");
178
179
180 TFile out(outputFile.c_str(), "recreate");
181
182 for (Int_t ix = 1; ix <= x_axis->GetNbins(); ++ix) {
183
184 const int index = ix - 1;
185 const JModule& module = detector[index];
186
187 if (module.getFloor() == 0) {
188 continue;
189 }
190
191 DEBUG("Bin " << setw(4) << ix << " -> " << setw(8) << module.getID() << endl);
192
193 TH1D h1(MAKE_CSTRING(module.getID() << ".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
194
195 // start values
196
197 Double_t ymin = numeric_limits<double>::max();
198 Double_t ymax = numeric_limits<double>::lowest();
199 Double_t x0 = 0.0; // [rad]
200
201 for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
202
203 h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
204 h1.SetBinError (iy, h2->GetEntries() * 1.0e-7 * parameters.FACTOR);
205
206 const Double_t x = h1.GetBinCenter (iy);
207 const Double_t y = h1.GetBinContent(iy);
208
209 if (y < ymin) {
210 x0 = x;
211 }
212
213 if (y > ymax) { ymax = y; }
214 if (y < ymin) { ymin = y; }
215 }
216
217 if (ymax - ymin >= parameters.YMIN) {
218
219 Double_t xmin = x0 + parameters.XMIN;
220 Double_t xmax = x0 + parameters.XMAX;
221
222 Double_t Rx = 0.0;
223
224 if (xmin < y_axis->GetXmin()) { Rx = +PI; }
225 if (xmax > y_axis->GetXmax()) { Rx = -PI; }
226
227 // apply offset
228
229 if (Rx != 0.0) {
230
231 NOTICE("Module " << setw(8) << module.getID() << " offset " << Rx/PI << " [pi]" << endl);
232
233 x0 += Rx;
234 xmin += Rx;
235 xmax += Rx;
236
237 h1.Reset();
238
239 const Double_t dx = (y_axis->GetXmax() - y_axis->GetXmin()) / y_axis->GetNbins();
240
241 for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
242
243 Int_t i = iy + (Int_t) (Rx / dx);
244
245 if (i < 1) { i += (Int_t) (2.0*PI / dx); }
246 if (i > y_axis->GetNbins()) { i -= (Int_t) (2.0*PI / dx); }
247
248 h1.SetBinContent(i, h2->GetBinContent(ix,iy));
249 h1.SetBinError (i, h2->GetEntries() * 1.0e-7 * parameters.FACTOR);
250 }
251 }
252
253 if (xmin < y_axis->GetXmin()) { xmin = y_axis->GetXmin(); xmax = xmin + 0.8*PI; }
254 if (xmax > y_axis->GetXmax()) { xmax = y_axis->GetXmax(); xmin = xmax - 0.8*PI; }
255
256 TF1* f1 = new TF1("f1", "[0] - [1]*cos(x - [2])");
257
258 f1->SetParameter(0, 0.5 * (ymin + ymax));
259 f1->SetParameter(1, 0.5 * (ymax - ymin));
260 f1->SetParameter(2, x0);
261
262 f1->SetParError(0, 0.01*(ymin+ymax));
263 f1->SetParError(1, 0.01*(ymin+ymax));
264 f1->SetParError(2, 0.05);
265
266 DEBUG("Start values:" << endl);
267
268 for (int i = 0; i != f1->GetNpar(); ++i) {
269 DEBUG(left << setw(12) << f1->GetParName (i) << ' ' <<
270 SCIENTIFIC(12,5) << f1->GetParameter(i) << endl);
271 }
272
273 // fit
274
275 TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);
276
277 const bool status = result->IsValid();
278
279 // revert offset
280
281 if (Rx != 0.0) {
282 f1->SetParameter(2, f1->GetParameter(2) - Rx);
283 }
284
285 if (debug >= notice_t || !result->IsValid()) {
286 cout << "Slice: "
287 << setw(4) << ix << ' '
288 << setw(16) << h1.GetName() << ' '
289 << FIXED(7,3) << f1->GetParameter(2) << " +/- "
290 << FIXED(7,3) << f1->GetParError(2) << " [rad] "
291 << FIXED(9,2) << result->Chi2() << '/'
292 << result->Ndf() << ' '
293 << result->IsValid()
294 << (status ? "" : "failed") << endl;
295 }
296
297 if (status) {
298 zmap[index] = result_type(f1->GetParameter(2), f1->GetParError(2));
299 }
300
301 struct {
302 const int ix;
303 const int iy;
304 } H2 = {
305 strings.getIndex(module.getString()) + 1,
306 module .getFloor()
307 };
308
309 if (status) {
310 H2d .SetBinContent(H2.ix, H2.iy, f1->GetParameter(2));
311 }
312 if (result->Ndf() > 0) {
313 H2c->SetBinContent(H2.ix, H2.iy, result->Chi2() / result->Ndf());
314 }
315 {
316 H2q->SetBinContent(H2.ix, H2.iy, status ? +1.0 : -1.0);
317 }
318
319 h1.Write();
320
321 delete f1;
322 }
323 }
324
325 if (zmap.empty()) {
326 NOTICE("No calibration results." << endl);
327 }
328
329 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
330 h0.SetBinContent(i->first, i->second.value);
331 h0.SetBinError (i->first, i->second.error);
332 }
333
334 if (overwriteDetector) {
335
336 NOTICE("Store calibration data on file " << detectorFile << endl);
337
338 detector.comment.add(JMeta(argc, argv));
339
340 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
341
342 JModule& module = detector[i->first];
343 const JVector3D center = module.getPosition();
344
345 module.sub(center);
346
347 module.rotate(JRotation3Z(i->second.value));
348
349 module.add(center);
350
351 if (overwriteDetector > 1) {
352 module.setQuaternion(module.getQuaternion() * JQuaternion3Z(i->second.value));
353 }
354 }
355
356 store(detectorFile, detector);
357 }
358
359 h0 .Write();
360 h2 ->Write();
361 H2d .Write();
362 H2c->Write();
363 H2q->Write();
364
365 out.Close();
366}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
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
gmake_property
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Definition JProperties.hh:1139
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition JModuleRouter.hh:38
JDETECTOR::JModule
Data structure for a composite optical module.
Definition JModule.hh:76
JEEP::JProperties
Utility class to parse parameter values.
Definition JProperties.hh:503
JGEOMETRY3D::JVector3D
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
JLANG::JException
General exception.
Definition JException.hh:24
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JPARSER::JCounter
Auxiliary class to handle multiple boolean-like I/O.
Definition JParser.hh:423
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JDETECTOR::getRangeOfFloors
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Definition JDetectorToolkit.hh:295
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:520
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition JDetectorToolkit.hh:617
JEEP::debug_t
@ debug_t
debug
Definition JMessage.hh:29
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::detector
Detector file.
Definition JHead.hh:227
JDETECTOR::JStringRouter
Router for mapping of string identifier to index.
Definition JStringRouter.hh:23
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488
Generated by doxygen 1.10.0