66int main(
int argc,
char **argv)
87 JParser<> zap(
"Program to fit time-slewing histogram in output of JCalibrateMuon.");
89 zap[
'f'] =
make_field(inputFile,
"input files (output from JCalibrateMuon).");
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.");
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.") =
"";
101 catch(
const exception &error) {
102 FATAL(error.what() << endl);
107 FATAL(
"Invalid fit range [ns] " << T_ns << endl);
112 if (option.find(
'S') == string::npos) { option +=
'S'; }
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,
127 100.5, 120.5, 140.5, 160.5, 180.5, 200.5, 250.5
130 const int N =
sizeof(x)/
sizeof(x[0]);
134 for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
136 NOTICE(
"Processing " << *i << endl);
138 TFile in(i->c_str(),
"exist");
141 FATAL(
"File " << *i <<
" not opened." << endl);
144 TH2D* p =
dynamic_cast<TH2D*
>(in.Get(
h2_t));
147 FATAL(
"File " << *i <<
" has no histogram <" <<
h2_t <<
">." << endl);
151 h2 = (TH2D*) p->Clone();
161 FATAL(
"No histogram <" <<
h2_t <<
">." << endl);
166 TH1D h0(
"h0", NULL, N-1, x);
170 for (
int i = 1; i <= h0.GetXaxis()->GetNbins(); ++i) {
172 const Int_t imin = h2->GetXaxis()->FindBin(h0.GetXaxis()->GetBinLowEdge(i));
173 const Int_t imax = h2->GetXaxis()->FindBin(h0.GetXaxis()->GetBinUpEdge (i));
175 TH1D h1(
MAKE_CSTRING(i <<
".1D"), NULL, h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax());
181 Double_t sigma = 3.5;
183 for (
int iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {
185 Double_t x = h1.GetBinCenter (iy);
188 for (
int ix = imin; ix != imax; ++ix) {
189 y += h2->GetBinContent(ix,iy);
192 h1.SetBinContent(iy, y);
203 if (function == gauss_t) {
205 f1 =
new TF1(function.c_str(),
"[0]*TMath::Gaus(x,[1],[2]) + [3]");
207 f1->SetParameter(0, ymax);
208 f1->SetParameter(1, x0);
209 f1->SetParameter(2, sigma);
210 f1->SetParameter(3, 0.0);
212 f1->SetParLimits(3, 0.0, ymax);
214 }
else if (function == emg_t) {
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]");
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);
224 f1->SetParLimits(4, 0.0, ymax);
228 FATAL(
"Unknown fit function " << function << endl);
236 if (xmin < h1.GetXaxis()->GetXmin()) { xmin = h1.GetXaxis()->GetXmin(); }
237 if (xmax > h1.GetXaxis()->GetXmax()) { xmax = h1.GetXaxis()->GetXmax(); }
239 TFitResultPtr
result = h1.Fit(f1, option.c_str(),
"same", xmin, xmax);
241 if (
result.Get() == NULL) {
242 FATAL(
"Invalid TFitResultPtr " << h1.GetName() << endl);
247 << setw(4) << i <<
' '
248 <<
FIXED(7,3) << f1->GetParameter(1) <<
" +/- "
249 <<
FIXED(7,3) << f1->GetParError (1) <<
' '
252 << (
result->IsValid() ?
"" :
"failed") << endl;
256 h0.SetBinContent(i, f1->GetParameter(1));
257 h0.SetBinError (i, f1->GetParError (1));
269 TF1 f1(
"f1", getRisetime, x[0], x[N-1], 3);
271 f1.SetParameter(0, 0.0);
272 f1.SetParameter(1, 2.0e-2);
273 f1.FixParameter(2, -7.0);
278 cout <<
"Time-slewing correction" << endl;
281 << (
result->IsValid() ?
"" :
"failed") << endl;
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;
288 cout <<
"// Produced by JFrodo.cc; to be included on JGetRiseTime.hh." << endl;
292 for (
int i = 0; i != 256; ++i) {
293 cout <<
"this->push_back(" <<
FIXED(6,2) << f1.Eval((Double_t) i) - t0 <<
");" << endl;