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.