Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JFrodo.cc File Reference

Fit time-slewing histogram in output of JCalibrateMuon.cc in calibration mode. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TMath.h"
#include "TFitResult.h"
#include "JROOT/JMinimizer.hh"
#include "JTools/JRange.hh"
#include "JDetector/JCalibration.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.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

Fit time-slewing histogram in output of JCalibrateMuon.cc in calibration mode.

Author
mdejong

Definition in file JFrodo.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 66 of file JFrodo.cc.

67{
68 using namespace std;
69 using namespace JPP;
70
71 typedef JRange<double> JRange_t;
72
73 vector<string> inputFile;
74 string outputFile;
75 string function;
76 JRange_t T_ns;
77 bool writeFits;
78 JPMTParameters parameters;
79 JRange_t X;
80 string option;
81 int debug;
82
83 try {
84
85 JProperties properties = parameters.getProperties();
86
87 JParser<> zap("Program to fit time-slewing histogram in output of JCalibrateMuon.");
88
89 zap['f'] = make_field(inputFile, "input files (output from JCalibrateMuon).");
90 zap['o'] = make_field(outputFile, "output file.") = "frodo.root";
91 zap['F'] = make_field(function, "fit function") = gauss_t, emg_t;
92 zap['T'] = make_field(T_ns, "ROOT fit range around maximum [ns].") = JRange_t(-7.5,+7.5);
93 zap['w'] = make_field(writeFits, "write fit results.");
94 zap['P'] = make_field(properties) = JPARSER::initialised();
95 zap['x'] = make_field(X, "ROOT fit range for time-over threshold") = JRange_t(0.0, 256.0);
96 zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
97 zap['d'] = make_field(debug) = 1;
98
99 zap(argc, argv);
100 }
101 catch(const exception &error) {
102 FATAL(error.what() << endl);
103 }
104
105
106 if (!T_ns.is_valid()) {
107 FATAL("Invalid fit range [ns] " << T_ns << endl);
108 }
109
110 cpu.setPMTParameters(parameters);
111
112 if (option.find('S') == string::npos) { option += 'S'; }
113 if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
114
115 // rebin time-over-threshold
116
117 const Double_t x[] = {
118 -0.5, 3.5, 6.5, 9.5, 12.5, 15.5, 18.5,
119 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
120 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
121 40.5, 42.5, 44.5, 46.5, 48.5,
122 50.5, 52.5, 54.5, 56.5, 58.5,
123 60.5, 65.5,
124 70.5, 75.5,
125 80.5, 85.5,
126 90.5, 95.5,
127 100.5, 120.5, 140.5, 160.5, 180.5, 200.5, 250.5
128 };
129
130 const int N = sizeof(x)/sizeof(x[0]);
131
132 TH2D* h2 = NULL;
133
134 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
135
136 NOTICE("Processing " << *i << endl);
137
138 TFile in(i->c_str(), "exist");
139
140 if (!in.IsOpen()) {
141 FATAL("File " << *i << " not opened." << endl);
142 }
143
144 TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
145
146 if (p == NULL) {
147 FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
148 }
149
150 if (h2 == NULL)
151 h2 = (TH2D*) p->Clone();
152 else
153 h2->Add(p);
154
155 h2->SetDirectory(0);
156
157 in.Close();
158 }
159
160 if (h2 == NULL) {
161 FATAL("No histogram <" << h2_t << ">." << endl);
162 }
163
164 // histogram fit results
165
166 TH1D h0("h0", NULL, N-1, x);
167
168 TFile out(outputFile.c_str(), "recreate");
169
170 for (int i = 1; i <= h0.GetXaxis()->GetNbins(); ++i) {
171
172 const Int_t imin = h2->GetXaxis()->FindBin(h0.GetXaxis()->GetBinLowEdge(i));
173 const Int_t imax = h2->GetXaxis()->FindBin(h0.GetXaxis()->GetBinUpEdge (i));
174
175 TH1D h1(MAKE_CSTRING(i << ".1D"), NULL, h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax());
176
177 // start values
178
179 Double_t ymax = 0.0;
180 Double_t x0 = 0.0; // [ns]
181 Double_t sigma = 3.5; // [ns]
182
183 for (int iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {
184
185 Double_t x = h1.GetBinCenter (iy);
186 Double_t y = 0.0;
187
188 for (int ix = imin; ix != imax; ++ix) {
189 y += h2->GetBinContent(ix,iy);
190 }
191
192 h1.SetBinContent(iy, y);
193
194 if (y > ymax) {
195 ymax = y;
196 x0 = x;
197 }
198 }
199
200
201 TF1* f1 = NULL;
202
203 if (function == gauss_t) {
204
205 f1 = new TF1(function.c_str(), "[0]*TMath::Gaus(x,[1],[2]) + [3]");
206
207 f1->SetParameter(0, ymax);
208 f1->SetParameter(1, x0);
209 f1->SetParameter(2, sigma);
210 f1->SetParameter(3, 0.0);
211
212 f1->SetParLimits(3, 0.0, ymax);
213
214 } else if (function == emg_t) {
215
216 f1 = new TF1(function.c_str(), "[0]*TMath::Exp(0.5*[3]*(2.0*[1]+[3]*[2]*[2]-2.0*x))*TMath::Erfc(([1]+[3]*[2]*[2]-x)/(TMath::Sqrt(2.0)*[2])) +[4]");
217
218 f1->SetParameter(0, ymax);
219 f1->SetParameter(1, x0 - 0.5*sigma);
220 f1->SetParameter(2, sigma);
221 f1->SetParameter(3, 0.04);
222 f1->SetParameter(4, 0.0);
223
224 f1->SetParLimits(4, 0.0, ymax);
225
226 } else {
227
228 FATAL("Unknown fit function " << function << endl);
229 }
230
231 // fit
232
233 Double_t xmin = x0 + T_ns.getLowerLimit();
234 Double_t xmax = x0 + T_ns.getUpperLimit();
235
236 if (xmin < h1.GetXaxis()->GetXmin()) { xmin = h1.GetXaxis()->GetXmin(); }
237 if (xmax > h1.GetXaxis()->GetXmax()) { xmax = h1.GetXaxis()->GetXmax(); }
238
239 TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);
240
241 if (result.Get() == NULL) {
242 FATAL("Invalid TFitResultPtr " << h1.GetName() << endl);
243 }
244
245 if (debug >= notice_t || !result->IsValid()) {
246 cout << "Bin: "
247 << setw(4) << i << ' '
248 << FIXED(7,3) << f1->GetParameter(1) << " +/- "
249 << FIXED(7,3) << f1->GetParError (1) << ' '
250 << FIXED(7,3) << result->Chi2() << '/'
251 << result->Ndf() << ' '
252 << (result->IsValid() ? "" : "failed") << endl;
253 }
254
255 if (result->IsValid()) {
256 h0.SetBinContent(i, f1->GetParameter(1));
257 h0.SetBinError (i, f1->GetParError (1));
258 }
259
260 if (writeFits) {
261 h1.Write();
262 }
263
264 delete f1;
265 }
266
267 // Fit
268
269 TF1 f1("f1", getRisetime, x[0], x[N-1], 3);
270
271 f1.SetParameter(0, 0.0);
272 f1.SetParameter(1, 2.0e-2);
273 f1.FixParameter(2, -7.0);
274
275 TFitResultPtr result = h0.Fit(&f1, "S", "same", X.constrain(x[0]), X.constrain(x[N-1]));
276
277 if (debug >= notice_t || !result->IsValid()) {
278 cout << "Time-slewing correction" << endl;
279 cout << FIXED(7,3) << result->Chi2() << '/'
280 << result->Ndf() << ' '
281 << (result->IsValid() ? "" : "failed") << endl;
282
283 for (int i = 0; i != f1.GetNpar(); ++i) {
284 cout << "parameter[" << i << "] = " << FIXED(8,4) << f1.GetParameter(i) << " +/- " << FIXED(8,4) << f1.GetParError (i) << endl;
285 }
286 }
287
288 cout << "// Produced by JFrodo.cc; to be included on JGetRiseTime.hh." << endl;
289
290 const double t0 = f1.Eval(TIME_OVER_THRESHOLD_NS);
291
292 for (int i = 0; i != 256; ++i) {
293 cout << "this->push_back(" << FIXED(6,2) << f1.Eval((Double_t) i) - t0 << ");" << endl;
294 }
295
296 h0.Write();
297
298 out.Close();
299}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
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::JPMTParameters
Data structure for PMT parameters.
Definition JPMTParameters.hh:29
JDETECTOR::JPMTParameters::getProperties
JProperties getProperties(const JEquationParameters &equation=JPMTParameters::getEquationParameters())
Get properties of this class.
Definition JPMTParameters.hh:218
JEEP::JProperties
Utility class to parse parameter values.
Definition JProperties.hh:501
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
JTOOLS::JRange::is_valid
bool is_valid() const
Check validity of range.
Definition JRange.hh:311
JTOOLS::JRange::constrain
T constrain(argument_type x) const
Constrain value to range.
Definition JRange.hh:350
JTOOLS::JRange::getLowerLimit
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
GAUSS_HERMITE::sigma
const double sigma[]
Definition JQuadrature.cc:74
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::TIME_OVER_THRESHOLD_NS
const double TIME_OVER_THRESHOLD_NS
Specification for time-over-threshold corresponding to a one photo-electron pulse.
Definition JDetector/JCalibration.hh:30
JEEP::debug_t
@ debug_t
debug
Definition JMessage.hh:29
JEEP::notice_t
@ notice_t
notice
Definition JMessage.hh:32
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
JRECONSTRUCTION::h2_t
const char *const h2_t
Name of histogram with results from JMuonCompass.cc.
Definition JMuonCompass.hh:12
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
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Generated by doxygen 1.12.0