57int main(
int argc,
char **argv)
73 double LeftMargin = 10.0;
74 double RightMargin = 10.0;
75 double gradientThreshold = -0.005;
99 JParser<> zap(
"Auxiliary program to fit time-over-threshold distributions.");
101 zap[
'f'] =
make_field(inputFile,
"input file (output from JCalibrateToT).");
103 zap[
'a'] =
make_field(detectorFile,
"detector file.");
104 zap[
'P'] =
make_field(pmtFile,
"specify PMT file name that can be given to JTriggerEfficiency.") =
"";
105 zap[
'w'] =
make_field(writeFits,
"write fit results.");
106 zap[
'O'] =
make_field(option,
"ROOT fit options, see TH1::Fit") =
"";
108 zap[
'x'] =
make_field(ToTfitRange,
"ROOT fit range (time-over-threshold).") =
JRange_t(0.0, 35.0);
114 catch(
const exception &error) {
115 FATAL(error.what() << endl);
135 if (!pmtFile.empty()) {
138 parametersMap.
load(pmtFile.c_str());
144 if (option.find(
'R') == string::npos) { option +=
'R'; }
145 if (option.find(
'S') == string::npos) { option +=
'S'; }
149 TFile* in = TFile::Open(inputFile.c_str(),
"exist");
151 if (in == NULL || !in->IsOpen()) {
152 FATAL(
"File: " << inputFile <<
" not opened." << endl);
156 TH2D gain(
"gain", NULL,
159 TH1D chi2(
"chi2", NULL, 100, 0.0, 5.0);
173 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
175 NOTICE(
"Module " << setw(10) << module->getID() <<
' ' << module->getLocation() << endl);
177 if (module->empty()) {
189 WARNING(
"No histogram " << module->getID() <<
"; skip fit." << endl);
194 const double ny = h2s->GetYaxis()->GetNbins();
195 const double ymin = h2s->GetYaxis()->GetXmin();
196 const double ymax = h2s->GetYaxis()->GetXmax();
198 TH1D chi2_1d (
MAKE_CSTRING(module->getID() <<
".1chi2"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
199 TH1D gain_1d (
MAKE_CSTRING(module->getID() <<
".1gain"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
200 TH1D gainspread_1d(
MAKE_CSTRING(module->getID() <<
".1gainspread"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
202 for (
int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
206 const string name =
MAKE_STRING(module->getID() <<
'.' <<
FILL(2,
'0') <<
209 DEBUG(
"Processing histogram " << name << endl);
211 TH1D h1(name.c_str(), NULL, ny, ymin, ymax);
215 for (
int iy = h2s->GetNbinsY(); iy >= 1; --iy) {
217 const double w = h1.GetBinWidth (iy);
218 const double y = h2s->GetBinContent(ix, iy);
220 h1.SetBinContent(iy, y/w);
221 h1.SetBinError (iy, sqrt(y)/w);
226 const double weight = h1.Integral(h1.FindBin(range.
getLowerLimit()),
229 if (weight <= Wmin) {
231 WARNING(
"Insufficient histogram entries for PMT " <<
id <<
232 "(" << weight <<
" < " << Wmin <<
"); skip fit." << endl);
251 RIGHT(10) <<
"threshold" <<
RIGHT(10) <<
"Mode" <<
RIGHT(10) <<
"Max" << endl);
253 for (
int i = h1.GetBinCenter(ny-1);
254 i > h1.GetBinCenter(h1.FindBin(parameters.
mean_ns + parameters.
sigma_ns));
257 const double x = h1.GetBinCenter (i);
258 const double y = h1.GetBinContent(i);
260 const double gradient = ( (h1.GetBinContent(i-1) - h1.GetBinContent(i+1)) /
261 (h1.GetBinCenter (i+1) - h1.GetBinCenter (i-1)) );
264 FIXED(10,1) << -fabs(gradientThreshold) * h1.Integral() <<
FIXED(10,1) << ToTmode <<
FIXED(10,1) << max << endl);
272 }
else if (gradient < -fabs(gradientThreshold) * h1.Integral()) {
283 range.
setRange(ToTmode - LeftMargin,
284 ToTmode + RightMargin);
296 JFitToT fit(parameters, range);
304 const double initGainSpread = initGain / 3.0;
306 if (!gainLimits(initGain)) {
308 WARNING(
"Invalid initial gain for PMT " <<
id <<
"; set default values." << endl);
323 const TFitResultPtr
result = fit(h1, option);
325 if (
result.Get() == NULL) {
326 FATAL(
"Invalid TFitResultPtr " << h1.GetName() << endl);
331 WARNING(
"Fit failure for PMT " <<
id <<
"; set initialization values." << endl);
333 h1.GetListOfFunctions()->Clear();
335 parameters.
gain = initGain;
346 (option.find(
"0") == string::npos && option.find(
"N") == string::npos)) {
355 JFitToT::getNumberOfModelParameters());
356 f1->SetParameters(fit.GetParameters());
357 f1->SetLineStyle(kDotted);
360 h1.GetListOfFunctions()->Add(f1);
365 const double reducedChi2 =
result->Chi2() /
result->Ndf();
370 chi2_1d. SetBinContent(ix, reducedChi2);
371 gain_1d. SetBinContent(ix, values.
gain);
372 gain_1d. SetBinError (ix, errors.
gain);
373 gainspread_1d.SetBinContent(ix, values.
gainSpread);
374 gainspread_1d.SetBinError (ix, errors.
gainSpread);
377 chi2.Fill(TMath::Log10(reducedChi2));
379 NOTICE(
"PMT " <<
id <<
" chi2 / NDF " <<
result->Chi2() <<
' ' <<
result->Ndf() << endl);
401 const JMeta* meta = in.next();
403 buffer.push_back(*meta);
408 for (vector<JMeta>::const_reverse_iterator i = buffer.rbegin(); i != buffer.rend(); ++i) {
419 parametersMap.
store(pmtFile.c_str());