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());
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()),
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()) {
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);
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;
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);
403 buffer.push_back(*meta);
409 parametersMap.comment.add(*
i);
419 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)
macros 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.
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.
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.
const JPMTAnalogueSignalProcessor & getCPU() const
Access method for the analogue signal processor.
static const char WILDCARD
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.
then fatal The output file must have the wildcard in the name
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.
then fatal The output file must have the wildcard in the e g root fi 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
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.
double saturation
saturation [ns]
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.
static const std::string FITTOT_FNAME
#define DEBUG(A)
Message macros.