Jpp test-rotations-old
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/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 )

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

36{
37 using namespace std;
38 using namespace JPP;
39
40 vector<string> inputFile;
41 string outputFile;
42 string detectorFile;
43 bool overwriteDetector;
44 string option;
45 double factor;
46 int debug;
47
48 try {
49
50 JParser<> zap("Program to fit histogram in output of JMuonCompass.cc.");
51
52 zap['f'] = make_field(inputFile, "input files (output from JMuonCompass).");
53 zap['o'] = make_field(outputFile, "output file.") = "hobbit.root";
54 zap['a'] = make_field(detectorFile, "detector file.");
55 zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
56 zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
57 zap['E'] = make_field(factor, "scaling factor for error on bin contents") = 1.0;
58 zap['d'] = make_field(debug) = 1;
59
60 zap(argc, argv);
61 }
62 catch(const exception &error) {
63 FATAL(error.what() << endl);
64 }
65
66
67 gErrorIgnoreLevel = kFatal;
68
69
70 JDetector detector;
71
72 try {
73 load(detectorFile, detector);
74 }
75 catch(const JException& error) {
76 FATAL(error);
77 }
78
79 const JModuleRouter router(detector);
80
81 if (option.find('S') == string::npos) { option += 'S'; }
82 if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
83
84 TH2D* h2 = NULL;
85
86 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
87
88 NOTICE("Processing " << *i << endl);
89
90 TFile in(i->c_str(), "exist");
91
92 if (!in.IsOpen()) {
93 FATAL("File " << *i << " not opened." << endl);
94 }
95
96 TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
97
98 if (p == NULL) {
99 FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
100 }
101
102 if (h2 == NULL)
103 h2 = (TH2D*) p->Clone();
104 else
105 h2->Add(p);
106
107 h2->SetDirectory(0);
108
109 in.Close();
110 }
111
112 if (h2 == NULL) {
113 FATAL("No histogram <" << h2_t << ">." << endl);
114 }
115
116 // auxiliary data structure for fit result
117
118 struct result_type {
119
120 result_type() :
121 value(0.0),
122 error(0.0)
123 {}
124
125 result_type(double value,
126 double error) :
127 value(value),
128 error(error)
129 {}
130
131 double value;
132 double error;
133 };
134
135 typedef map<int, result_type> map_type; // index -> fit result
136
137 map_type zmap;
138
139
140 // histogram fit results
141
142 const TAxis* x_axis = h2->GetXaxis();
143 const TAxis* y_axis = h2->GetYaxis();
144
145 TH1D h0("h0", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
146 TH1D hc("hc", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
147 TH1D hq("hq", NULL, x_axis->GetNbins(), -0.5, x_axis->GetNbins() - 0.5);
148
149 const floor_range floors = getRangeOfFloors(detector);
150 const JStringRouter strings(detector);
151
152 TH2D H2("detector", NULL,
153 strings.size(), -0.5, strings.size() - 0.5,
154 floors.getUpperLimit(), 1 - 0.5, floors.getUpperLimit() + 0.5);
155
156 for (Int_t ix = 1; ix <= H2.GetXaxis()->GetNbins(); ++ix) {
157 H2.GetXaxis()->SetBinLabel(ix, MAKE_CSTRING(strings.at(ix-1)));
158 }
159 for (Int_t iy = 1; iy <= H2.GetYaxis()->GetNbins(); ++iy) {
160 H2.GetYaxis()->SetBinLabel(iy, MAKE_CSTRING(iy));
161 }
162
163 TFile out(outputFile.c_str(), "recreate");
164
165 for (Int_t ix = 1; ix <= x_axis->GetNbins(); ++ix) {
166
167 const int index = ix - 1;
168 const int id = detector[index].getID();
169
170 if (detector[index].getFloor() == 0) {
171 continue;
172 }
173
174 DEBUG("Bin " << setw(4) << ix << " -> " << setw(8) << id << endl);
175
176 TH1D h1(MAKE_CSTRING(id << ".1D"), NULL, y_axis->GetNbins(), y_axis->GetXmin(), y_axis->GetXmax());
177
178 // start values
179
180 Double_t ymin = numeric_limits<double>::max();
181 Double_t ymax = numeric_limits<double>::lowest();
182 Double_t x0 = 0.0; // [rad]
183
184 for (Int_t iy = 1; iy <= y_axis->GetNbins(); ++iy) {
185
186 h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
187 h1.SetBinError (iy, h2->GetEntries() * 1.0e-7 * factor);
188
189 const Double_t x = h1.GetBinCenter (iy);
190 const Double_t y = h1.GetBinContent(iy);
191
192 if (y < ymin) {
193 x0 = x;
194 }
195
196 if (y > ymax) { ymax = y; }
197 if (y < ymin) { ymin = y; }
198 }
199
200 Double_t xmin = x0 - 0.4*PI;
201 Double_t xmax = x0 + 0.4*PI;
202
203 if (xmin < y_axis->GetXmin()) { xmin = y_axis->GetXmin(); xmax = xmin + 0.8*PI; }
204 if (xmax > y_axis->GetXmax()) { xmax = y_axis->GetXmax(); xmin = xmax - 0.8*PI; }
205
206 TF1* f1 = new TF1("f1", "[0] - [1]*cos(x - [2])");
207
208 f1->SetParameter(0, 0.5 * (ymin + ymax));
209 f1->SetParameter(1, 0.5 * (ymax - ymin));
210 f1->SetParameter(2, x0);
211
212 f1->SetParError(0, 0.01*(ymin+ymax));
213 f1->SetParError(1, 0.01*(ymin+ymax));
214 f1->SetParError(2, 0.05);
215
216 DEBUG("Start values:" << endl);
217
218 for (int i = 0; i != f1->GetNpar(); ++i) {
219 DEBUG(left << setw(12) << f1->GetParName (i) << ' ' <<
220 SCIENTIFIC(12,5) << f1->GetParameter(i) << endl);
221 }
222
223 // fit
224
225 TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);
226
227 const bool status = result->IsValid();
228
229 if (debug >= notice_t || !result->IsValid()) {
230 cout << "Slice: "
231 << setw(4) << ix << ' '
232 << setw(16) << h1.GetName() << ' '
233 << FIXED(7,3) << f1->GetParameter(2) << " +/- "
234 << FIXED(7,3) << f1->GetParError(2) << ' '
235 << FIXED(9,2) << result->Chi2() << '/'
236 << result->Ndf() << ' '
237 << result->IsValid()
238 << (status ? "" : "failed") << endl;
239 }
240
241 if (status) {
242 zmap[index] = result_type(f1->GetParameter(2), f1->GetParError(2));
243 }
244
245 if (result->Ndf() > 0) {
246 hc.SetBinContent(ix, result->Chi2() / result->Ndf());
247 }
248
249 hq.SetBinContent(ix, status ? 1.0 : 0.0);
250
251 h1.Write();
252
253 delete f1;
254 }
255
256 if (zmap.empty()) {
257 NOTICE("No calibration results." << endl);
258 }
259
260 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
261 h0.SetBinContent(i->first, i->second.value);
262 h0.SetBinError (i->first, i->second.error);
263 }
264
265 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
266
267 const JModule& module = detector[i->first];
268
269 const int ix = i->first + 1;
270 const int iy = module.getFloor();
271
272 H2.SetBinContent(ix, iy, i->second.value);
273 }
274
275 if (overwriteDetector) {
276
277 NOTICE("Store calibration data on file " << detectorFile << endl);
278
279 detector.comment.add(JMeta(argc, argv));
280
281 for (map_type::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
282
283 JModule& module = detector[i->first];
284 const JVector3D center = module.getPosition();
285
286 module.sub(center);
287
288 module.rotate(JRotation3Z(i->second.value));
289
290 module.add(center);
291 }
292
293 store(detectorFile, detector);
294 }
295
296 h0 .Write();
297 hc .Write();
298 hq .Write();
299 h2->Write();
300 H2 .Write();
301
302 out.Close();
303}
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
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:75
JGEOMETRY3D::JVector3D
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
JLANG::JException
General exception.
Definition JException.hh:24
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
GAUSS_LEGENDRE::xmax
const double xmax
Definition JQuadrature.cc:24
GAUSS_LEGENDRE::xmin
const double xmin
Definition JQuadrature.cc:23
JDETECTOR::getRangeOfFloors
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Definition JDetectorToolkit.hh:267
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:508
JDETECTOR::store
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
Definition JDetectorToolkit.hh:605
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:43
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
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.12.0