10 #include "TFitResult.h"
46 static const char*
const WILDCARD =
"%";
57 int main(
int argc,
char **argv)
61 using namespace KM3NETDAQ;
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());
143 parametersMap.comment.add(
JMeta(argc, argv));
146 if (option.find(
'R') == string::npos) { option +=
'R'; }
147 if (option.find(
'S') == string::npos) { option +=
'S'; }
151 TFile*
in = TFile::Open(inputFile.c_str(),
"exist");
153 if (
in == NULL || !
in->IsOpen()) {
154 FATAL(
"File: " << inputFile <<
" not opened." << endl);
158 TH2D
gain(
"gain", NULL,
161 TH1D
chi2(
"chi2", NULL, 100, 0.0, 5.0);
175 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
177 NOTICE(
"Module " << setw(10) << module->getID() <<
' ' << module->getLocation() << endl);
179 if (module->empty()) {
191 WARNING(
"No histogram " << module->getID() <<
"; skip fit." << endl);
196 const double ny = h2s->GetYaxis()->GetNbins();
197 const double ymin = h2s->GetYaxis()->GetXmin();
198 const double ymax = h2s->GetYaxis()->GetXmax();
204 for (
int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
211 DEBUG(
"Processing histogram " << name << endl);
213 TH1D h1(name.c_str(), NULL, ny, ymin, ymax);
217 for (
int iy = h2s->GetNbinsY(); iy >= 1; --iy) {
219 const double w = h1.GetBinWidth (iy);
220 const double y = h2s->GetBinContent(ix, iy);
222 h1.SetBinContent(iy, y/w);
223 h1.SetBinError (iy, sqrt(y)/w);
231 if (weight <= Wmin) {
233 WARNING(
"Insufficient histogram entries for PMT " <<
id <<
234 "(" << weight <<
" < " << Wmin <<
"); skip fit." << endl);
253 RIGHT(10) <<
"threshold" <<
RIGHT(10) <<
"Mode" <<
RIGHT(10) <<
"Max" << endl);
255 for (
int i = h1.GetBinCenter(ny-1);
256 i > h1.GetBinCenter(h1.FindBin(parameters.
mean_ns + parameters.
sigma_ns));
259 const double x = h1.GetBinCenter (i);
260 const double y = h1.GetBinContent(i);
262 const double gradient = ( (h1.GetBinContent(i-1) - h1.GetBinContent(i+1)) /
263 (h1.GetBinCenter (i+1) - h1.GetBinCenter (i-1)) );
266 FIXED(10,1) << -fabs(gradientThreshold) * h1.Integral() <<
FIXED(10,1) << ToTmode <<
FIXED(10,1) << max << endl);
274 }
else if (gradient < -fabs(gradientThreshold) * h1.Integral()) {
285 range.
setRange(ToTmode - LeftMargin,
286 ToTmode + RightMargin);
298 JFitToT fit(parameters, range);
306 const double initGainSpread = initGain / 3.0;
308 if (!gainLimits(initGain)) {
310 WARNING(
"Invalid initial gain for PMT " <<
id <<
"; set default values." << endl);
325 const TFitResultPtr
result = fit(h1, option);
327 if (result.Get() == NULL) {
328 FATAL(
"Invalid TFitResultPtr " << h1.GetName() << endl);
331 if (
int(result) < 0 || !result->IsValid()) {
333 WARNING(
"Fit failure for PMT " <<
id <<
"; set initialization values." << endl);
335 h1.GetListOfFunctions()->Clear();
337 parameters.
gain = initGain;
348 (option.find(
"0") == string::npos && option.find(
"N") == string::npos)) {
357 JFitToT::getNumberOfModelParameters());
358 f1->SetParameters(fit.GetParameters());
359 f1->SetLineStyle(kDotted);
362 h1.GetListOfFunctions()->Add(f1);
367 const double reducedChi2 = result->Chi2() / result->Ndf();
372 chi2_1d. SetBinContent(ix, reducedChi2);
373 gain_1d. SetBinContent(ix,
values.gain);
374 gain_1d. SetBinError (ix, errors.gain);
375 gainspread_1d.SetBinContent(ix,
values.gainSpread);
376 gainspread_1d.SetBinError (ix, errors.gainSpread);
379 chi2.Fill(TMath::Log10(reducedChi2));
381 NOTICE(
"PMT " <<
id <<
" chi2 / NDF " << result->Chi2() <<
' ' << result->Ndf() << endl);
400 const JMeta& meta = *
in.next();
402 parametersMap.comment.add(meta);
413 parametersMap.store(pmtFile.c_str());
Utility class to parse command line options.
const Double_t getModelParameter(const int i) const
Get model parameter.
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
int main(int argc, char *argv[])
ROOT TTree parameter settings of various packages.
TString replace(const TString &target, const TRegexp ®exp, const T &replacement)
Replace regular expression in input by given replacement.
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
double gainSpread
gain spread [unit]
static const double FITTOT_GAIN_MAX
Default maximal gain.
Double_t gainSpread
PMT gain spread.
const double TIME_OVER_THRESHOLD_NS
Specification for time-over-threshold corresponding to a one photo-electron pulse.
static const double FITTOT_GAIN_MIN
Default minimal gain.
Recording of objects on file according a format that follows from the file name extension.
then echo Enter input within $TIMEOUT_S seconds echo n User name
Utility class to parse parameter values.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
#define MAKE_CSTRING(A)
Make C-string.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary data structure for floating point format specification.
Parametrisation of time-over-threshold distribution.
Data structure for detector geometry and calibration.
Fit parameters for two-fold coincidence rate due to K40.
#define MAKE_STRING(A)
Make string.
const JPMTAnalogueSignalProcessor & getCPU() const
Access method for the analogue signal processor.
static const std::string FITTOT_SUFFIX
Scanning of objects from a single file according a format that follows from the extension of each fil...
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Type definition of range.
I/O formatting auxiliaries.
double mean_ns
mean time-over-threshold of threshold-band hits [ns]
double sigma_ns
time-over-threshold standard deviation of threshold-band hits [ns]
static const double FITTOT_GAINSPREAD_MAX
Default maximal gain spread.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
virtual double getNPE(const double tot_ns) const override
Get number of photo-electrons.
Auxiliary class for map of PMT parameters.
bool setParLimits(TF1 &f1, const Int_t index, Double_t xmin, Double_t xmax)
Set fit parameter limits.
General purpose messaging.
Auxiliary data structure for sequence of same character.
void load(const char *file_name)
Load from input file.
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
then $JPP_DIR examples JDetector JToT o $OUTPUT_FILE n N $NPE P gain
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
Utility class to parse command line options.
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
static const char *const _2SToT
Histogram naming.
PMT analogue signal processor.
static const double FITTOT_GAINSPREAD_MIN
Default minimal gain spread.
Object reading from a list of files.
Data structure for PMT parameters.
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
static const std::string FITTOT_FNAME
std::vector< double > weight
#define DEBUG(A)
Message macros.