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());
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) {
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()),
227 h1.FindBin(
range.getUpperLimit()));
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);
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()) {
282 if (!
range.is_valid()) {
283 range.setRange(ToTmode - LeftMargin,
284 ToTmode + RightMargin);
287 range.setRange(std::max(h1.GetXaxis()->GetXmin(),
range.getLowerLimit()),
288 std::min(h1.GetXaxis()->GetXmax(),
range.getUpperLimit()));
299 gainLimits.getLowerLimit(), gainLimits.getUpperLimit());
301 gainSpreadLimits.getLowerLimit(), gainSpreadLimits.getUpperLimit());
303 const double initGain = fit.getCPU().getNPE(ToTmode);
304 const double initGainSpread = initGain / 3.0;
306 if (!gainLimits(initGain)) {
308 WARNING(
"Invalid initial gain for PMT " <<
id <<
"; set default values." << endl);
310 parameters.
gain = gainLimits .constrain(initGain);
311 parameters.
gainSpread = gainSpreadLimits.constrain(initGainSpread);
321 fit.gainSpread = initGainSpread;
323 const TFitResultPtr
result = fit(h1, option);
325 if (result.Get() == NULL) {
326 FATAL(
"Invalid TFitResultPtr " << h1.GetName() << endl);
329 if (
int(result) < 0 || !result->IsValid()) {
331 WARNING(
"Fit failure for PMT " <<
id <<
"; set initialization values." << endl);
333 h1.GetListOfFunctions()->Clear();
335 parameters.
gain = initGain;
345 if ((ymin <
range.getLowerLimit() || ymax >
range.getUpperLimit()) &&
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);
376 gain.Fill(fit.gain, fit.gainSpread);
377 chi2.Fill(TMath::Log10(reducedChi2));
379 NOTICE(
"PMT " <<
id <<
" chi2 / NDF " << result->Chi2() <<
' ' << result->Ndf() << endl);
381 parameters.
gain = fit.gain;
403 buffer.push_back(*meta);
409 parametersMap.comment.add(*
i);
419 parametersMap.store(pmtFile.c_str());
Utility class to parse command line options.
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
TString replace(const TString &target, const TRegexp ®exp, const T &replacement)
Replace regular expression in input by given replacement.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
double gainSpread
gain spread [unit]
static const double FITTOT_GAIN_MAX
Default maximal gain.
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.
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.
then fatal Wrong number of parameters fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER set_variable DETID getDetector D $DETECTOR_ID O string $DIR JPlotPMTParameters a $DETECTOR P $INPUT_FILE o $OUTPUT_FILE d $DEBUG for KEY in TTS_ns QE gain gainSpread
Fit parameters for two-fold coincidence rate due to K40.
#define MAKE_STRING(A)
Make string.
static const char WILDCARD
static const std::string FITTOT_SUFFIX
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Type definition of range.
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
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.
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.
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
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.
static const double FITTOT_GAINSPREAD_MIN
Default minimal gain spread.
Object reading from a list of files.
Data structure for PMT parameters.
double saturation
saturation [ns]
do set_variable DETECTOR_TXT $WORKDIR detector
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
static const std::string FITTOT_FNAME
#define DEBUG(A)
Message macros.