Jpp in_tag_pdf_generation
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 )

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