10 #include "TFitResult.h"
46 static const char*
const WILDCARD =
"%";
57 int 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();
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);
419 parametersMap.
store(pmtFile.c_str());
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
Recording of objects on file according a format that follows from the file name extension.
int main(int argc, char **argv)
General purpose messaging.
#define DEBUG(A)
Message macros.
PMT analogue signal processor.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
TString replace(const TString &target, const TRegexp ®exp, const T &replacement)
Replace regular expression in input by given replacement.
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
#define MAKE_STRING(A)
Make string.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Auxiliary class to define a range between two values.
Scanning of objects from a single file according a format that follows from the extension of each fil...
ROOT TTree parameter settings of various packages.
Auxiliary class for map of PMT parameters.
Data structure for PMT parameters.
double sigma_ns
time-over-threshold standard deviation of threshold-band hits [ns]
double gainSpread
gain spread [unit]
double mean_ns
mean time-over-threshold of threshold-band hits [ns]
double saturation
saturation [ns]
Utility class to parse parameter values.
Utility class to parse command line options.
Object reading from a list of files.
virtual bool hasNext() override
Check availability of next element.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
static const char *const _2SToT
Histogram naming.
static const double FITTOT_GAINSPREAD_MAX
Default maximal gain spread.
static const std::string FITTOT_FNAME
static const std::string FITTOT_SUFFIX
static const double FITTOT_GAIN_MAX
Default maximal gain.
static const double FITTOT_GAIN_MIN
Default minimal gain.
static const double FITTOT_GAINSPREAD_MIN
Default minimal gain spread.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
const double TIME_OVER_THRESHOLD_NS
Specification for time-over-threshold corresponding to a one photo-electron pulse.
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool setParLimits(TF1 &f1, const Int_t index, Double_t xmin, Double_t xmax)
Set fit parameter limits.
KM3NeT DAQ data structures and auxiliaries.
static const char WILDCARD
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Auxiliary data structure for sequence of same character.
Auxiliary data structure for floating point format specification.
Type definition of range.
Fit parameters for two-fold coincidence rate due to K40.
Double_t gainSpread
PMT gain spread.
const Double_t getModelParameter(const int i) const
Get model parameter.
Parametrisation of time-over-threshold distribution.
const JPMTAnalogueSignalProcessor & getCPU() const
Access method for the analogue signal processor.
virtual double getNPE(const double tot_ns) const override
Get number of photo-electrons.
void store(const char *file_name) const
Store to output file.
void load(const char *file_name)
Load from input file.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...