81   using namespace KM3NETDAQ;
 
   94   bool      overwriteDetector;
 
  117     JParser<> zap(
"Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
 
  120     zap[
'f'] = 
make_field(inputFile,         
"input file (output from JMergeCalibrateK40).");
 
  122     zap[
'a'] = 
make_field(detectorFile,      
"detector file.");
 
  123     zap[
'P'] = 
make_field(pmtFile,           
"specify PMT file name that can be given to JTriggerEfficiency.") = 
"";
 
  125                           "Fix time offset(s) of PMT(s) of certain module(s), e.g." 
  126                           "\n-! \"808969848 0 808982077 23\" will fix offset of PMT 0 of module 808969848 and of PMT 23 of module 808982077." 
  127                           "\nSame PMT offset can be fixed for all optical modules, e.g." 
  128                           "\n-! \"-1 0 -1 23\" will fix offsets of PMTs 0 and 23 for all optical modules.")
 
  130     zap[
'r'] = 
make_field(reverse,           
"reverse TDC constraints due to option -! <TDC>.");
 
  131     zap[
'A'] = 
make_field(overwriteDetector, 
"overwrite detector file provided through '-a' with correct time offsets.");
 
  132     zap[
'w'] = 
make_field(writeFits,         
"write fit results.");
 
  133     zap[
'D'] = 
make_field(fitAngDist,        
"fit angular distribution; fix normalisation.");
 
  134     zap[
'M'] = 
make_field(fitModel,          
"fit angular distribution; fix PMT QEs = 1.0.");
 
  137     zap[
'O'] = 
make_field(option,            
"ROOT fit option, see TH1::Fit.")                                 = 
"";
 
  142   catch(
const exception &error) {
 
  143     FATAL(error.what() << endl);
 
  147   ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(1000000);
 
  154   for (JTDC_t::const_iterator 
i = 
TDC.begin(); 
i != 
TDC.end(); ++
i) {
 
  155     DEBUG(
"PMT " << setw(10) << 
i->first << 
' ' << setw(2) << 
i->second << 
" constrain t0." << endl);
 
  187   if (option.find(
'O') == string::npos) { option += 
'0'; }
 
  188   if (option.find(
'S') == string::npos) { option += 
'S'; }
 
  192   TFile* 
in = TFile::Open(inputFile.c_str(), 
"exist");
 
  194   if (
in == NULL || !
in->IsOpen()) {
 
  195     FATAL(
"File: " << inputFile << 
" not opened." << endl);
 
  199   TDirectory::AddDirectory(NULL);
 
  205   TH1D hc(
"chi2", NULL, 500, 0.0, 5.0);
 
  207   TH1D 
p1(
"p1", NULL, 501, -5.0, +5.0);
 
  208   TH1D 
p2(
"p2", NULL, 501, -5.0, +5.0);
 
  209   TH1D 
p3(
"p3", NULL, 501, -5.0, +5.0);
 
  210   TH1D 
p4(
"p4", NULL, 501, -5.0, +5.0);
 
  211   TH1D bg(
"bg", NULL, 501, -0.1, +0.1);
 
  212   TH1D cc(
"cc", NULL, 501, -0.1, +0.1);
 
  219   for (JDetector::iterator module = 
detector.begin(); module != 
detector.end(); ++module) {
 
  221     NOTICE(
"Module " << setw(10) << module->getID() << 
' ' << module->getLocation() << endl);
 
  223     if (module->empty()) {
 
  233     if (h2s == NULL || h2s->GetEntries() == 0) {
 
  235       NOTICE(
"No data for module " << module->getID() << 
"; set corresponding QEs to 0." << endl);
 
  255                 h2s->GetXaxis()->GetXmin(),
 
  256                 h2s->GetXaxis()->GetXmax(),
 
  257                 h2s->GetYaxis()->GetXmin(),
 
  258                 h2s->GetYaxis()->GetXmax(),
 
  259                 range.first == range.second);
 
  261     fit.setModelParameters(fitk40.getModelParameters());
 
  263     for (JTDC_t::const_iterator 
i = range.first; 
i != range.second; ++
i) {
 
  264       fixParameter(fit, fit.getModelParameter(
i->second, &JPMTParameters_t::t0));
 
  269         fixParameter(fit, fit.getModelParameter(pmt, &JPMTParameters_t::QE));
 
  274       fixParameter(fit, fit.getModelParameter(&JFitK40::Rate_Hz));
 
  290     for (
int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
 
  292       Double_t value = 0.0;
 
  293       Double_t error = 0.0;
 
  295       for (
int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
 
  296         value += h2s->GetBinContent(ix,iy);
 
  297         error += h2s->GetBinError(ix,iy) * h2s->GetBinError(ix,iy);
 
  302       if (value  <  MINIMAL_RATE_HZ + STDEV*error) {
 
  306         buffer[pair.first]  += 1;
 
  307         buffer[pair.second] += 1;
 
  317           WARNING(
"PMT " << setw(10) << module->getID() << 
' ' << setw(2) << pmt << 
" too low a rate; switch off." << endl);
 
  325       NOTICE(
"Skip fit for module " << module->getID() << 
"; set corresponding QEs to 0." << endl);
 
  337       h2s->GetXaxis()->SetRangeUser(max(
X.getLowerLimit(), h2s->GetXaxis()->GetXmin()),
 
  338                                     min(
X.getUpperLimit(), h2s->GetXaxis()->GetXmax()));
 
  341       h2s->GetYaxis()->SetRangeUser(max(
Y.getLowerLimit(), h2s->GetYaxis()->GetXmin()),
 
  342                                     min(
Y.getUpperLimit(), h2s->GetYaxis()->GetXmax()));
 
  347     TH1D h1(
MAKE_CSTRING(module->getID() << 
_1D), NULL, h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin(), h2s->GetXaxis()->GetXmax());
 
  349     const Double_t 
sigma[] = {
 
  351       (min(
Y.getUpperLimit(), h2s->GetYaxis()->GetXmax()) -
 
  352        max(
Y.getLowerLimit(), h2s->GetYaxis()->GetXmin())) / sqrt(12.0),
 
  358     for (Int_t ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
 
  362       Double_t value = 0.0;
 
  363       Double_t error = 0.0;
 
  365       for (Int_t iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
 
  367         const Double_t 
y = h2s->GetYaxis()->GetBinCenter(iy);
 
  368         const Double_t z = h2s->GetBinContent(ix, iy);
 
  369         const Double_t 
w = h2s->GetBinError  (ix, iy);
 
  385       h1.SetBinContent(ix, y0);
 
  386       h1.SetBinError  (ix, value > MINIMAL_RATE_HZ + STDEV * error ? sigma[1] * error / value : sigma[0]);
 
  392       TFitResultPtr 
result = 
f1(h1, option.c_str());
 
  394       if (result.Get() == NULL) {
 
  395         FATAL(
"Invalid TFitResultPtr " << h1.GetName() << endl);
 
  400         TH1D h2(
MAKE_CSTRING(module->getID() << 
_1F), NULL, h1.GetXaxis()->GetNbins(), h1.GetXaxis()->GetXmin(), h1.GetXaxis()->GetXmax());
 
  402         for (
int ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
 
  404           const Double_t 
x = h2.GetXaxis()->GetBinCenter(ix);
 
  406           h2.SetBinContent(ix, 
f1.Eval(x));
 
  407           h2.SetBinError  (ix, 0.0);
 
  415         cout << 
"Fit result chi2 / NDF " << result->Chi2() << 
' ' << result->Ndf() << 
' ' << (result->IsValid() ? 
"" : 
"failed") << endl;
 
  420                  << setw(2)    << pmt             << 
' ' 
  421                  << 
FIXED(7,2) << 
f1.getT0 (pmt)  << 
' ';
 
  422             cout << (result->IsParameterFixed(
f1.getModelParameter(pmt, &JPMTParameters_t::t0))  ? 
" *** fixed t0  *** " : 
"");
 
  423             cout << (
f1.is_disabled  (pmt) ? 
" (disabled)" : 
"");
 
  424             cout << (
f1.is_average_t0(pmt) ? 
" (averaged)" : 
"");
 
  429       fit.setModelParameters(
f1.getModelParameters());
 
  433     TFitResultPtr result = fit(*h2s, option);
 
  445         if (!result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE))  &&  !
is_valid(
values.getQE(pmt), errors.getQE(pmt))) {
 
  447           WARNING(
"PMT " << setw(10) << module->getID() << 
' ' << setw(2) << pmt << 
' ' 
  450                   << 
FIXED(5,2) << errors.getQE(pmt) 
 
  451                   << 
" compatible with zero -> set QE = 0 and t0 = 0; -> refit." << endl);
 
  462       NOTICE(
"Skip fit for module " << module->getID() << 
"; set corresponding QEs to 0." << endl);
 
  478       fit.setModelParameters(
f1.getModelParameters());
 
  480       result = fit(*h2s, option);
 
  487     hc.Fill(TMath::Log10(result->Chi2() / result->Ndf()));
 
  491              h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin (), h2s->GetXaxis()->GetXmax (),
 
  492              h2s->GetYaxis()->GetNbins(), h2s->GetYaxis()->GetXmin (), h2s->GetYaxis()->GetXmax ());
 
  494     for (
int ix = 1; ix <= h2f.GetXaxis()->GetNbins(); ++ix) {
 
  495       for (
int iy = 1; iy <= h2f.GetYaxis()->GetNbins(); ++iy) {
 
  497         const Double_t x = h2f.GetXaxis()->GetBinCenter(ix);
 
  498         const Double_t y = h2f.GetYaxis()->GetBinCenter(iy);
 
  500         h2f.SetBinContent(ix, iy, fit.Eval(x,y));
 
  501         h2f.SetBinError  (ix, iy, 0.0);
 
  508       cout << 
"Fit result chi2 / NDF " << result->Chi2() << 
' ' << result->Ndf() << 
' ' << (result->IsValid() ? 
"" : 
"failed") << endl;
 
  510       cout << 
"Rate_Hz= " << 
FIXED(8,4) << 
values.Rate_Hz << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::Rate_Hz)) ? 
" *** fixed *** " : 
"") << endl;
 
  512       cout << 
"p1=      " << 
FIXED(8,4) << 
values.p1      << (result->IsParameterFixed(fit.getModelParameter(&
JFitK40::p1))      ? 
" *** fixed *** " : 
"") << endl;
 
  513       cout << 
"p2=      " << 
FIXED(8,4) << 
values.p2      << (result->IsParameterFixed(fit.getModelParameter(&
JFitK40::p2))      ? 
" *** fixed *** " : 
"") << endl;
 
  514       cout << 
"p3=      " << 
FIXED(8,4) << 
values.p3      << (result->IsParameterFixed(fit.getModelParameter(&
JFitK40::p2))      ? 
" *** fixed *** " : 
"") << endl;
 
  515       cout << 
"p4=      " << 
FIXED(8,4) << 
values.p4      << (result->IsParameterFixed(fit.getModelParameter(&
JFitK40::p3))      ? 
" *** fixed *** " : 
"") << endl;
 
  517       cout << 
"Background (correlated):   " << 
FIXED(8,5) << 
values.bg << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::bg)) ? 
" *** fixed *** " : 
"") << endl;
 
  518       cout << 
"Background (uncorrelated): " << 
FIXED(8,5) << 
values.cc << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::cc)) ? 
" *** fixed *** " : 
"") << endl;
 
  520       cout << 
" PMT t0 [ns] TTS [ns]   R" << endl;
 
  527              << setw(2)    << pmt             << 
' ' 
  528              << 
FIXED(7,2) << fit.getT0 (pmt) << 
' ' 
  529              << 
FIXED(7,2) << fit.getTTS(pmt) << 
' ' 
  530              << 
FIXED(7,3) << fit.getQE (pmt);
 
  532         cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE))  ? 
" *** fixed QE  *** " : 
"");
 
  533         cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::TTS)) ? 
" *** fixed TTS *** " : 
"");
 
  534         cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::t0))  ? 
" *** fixed t0  *** " : 
"");
 
  535         cout << (fit.is_disabled  (pmt) ? 
" (disabled)" : 
"");
 
  536         cout << (fit.is_average_t0(pmt) ? 
" (averaged)" : 
"");
 
  539         Q[0].
put(fit.getT0 (pmt));
 
  540         Q[1].
put(fit.getTTS(pmt));
 
  541         Q[2].
put(fit.getQE (pmt));
 
  578         h1t.SetBinContent(pmt + 1, 
values.getT0 (pmt));
 
  579         h1s.SetBinContent(pmt + 1, 
values.getTTS(pmt));
 
  580         h1q.SetBinContent(pmt + 1, 
values.getQE (pmt));
 
  582         h1t.SetBinError  (pmt + 1, max(errors.getT0 (pmt), 
epsilon));
 
  583         h1s.SetBinError  (pmt + 1, max(errors.getTTS(pmt), 
epsilon));
 
  584         h1q.SetBinError  (pmt + 1, max(errors.getQE (pmt), 
epsilon));
 
  587       out << h1t << h1s << h1q;
 
  606       if (overwriteDetector) {
 
  607         module->getPMT(pmt).addT0(fit.getT0(pmt));
 
  617     out << hc << 
p1 << 
p2 << 
p3 << 
p4 << bg << cc;
 
  629       buffer.push_back(*meta);
 
  638   if (overwriteDetector) {
 
  640     NOTICE(
"Store calibration data on file " << detectorFile << endl);
 
Utility class to parse command line options. 
 
JCombinatorics::pair_type pair_type
 
Q(UTCMax_s-UTCMin_s)-livetime_s
 
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object 
 
bool is_valid(const json &js)
Check validity of JSon data. 
 
Auxiliary class for TDC constraints. 
 
Template specialisation of two-fold coincidence rate due to K40 and other radioactive decays...
 
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)...
 
Template specialisation of two-fold coincidence rate due to K40 and other radioactive decays...
 
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
 
Auxiliary data structure for floating point format specification. 
 
static const char *const _1F
Name extension for 1D time offset fit. 
 
const JPolynome f1(1.0, 2.0, 3.0)
Function. 
 
Type definition of range. 
 
static const char *const _2F
Name extension for 2F rate fitted. 
 
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object 
 
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file. 
 
static const char *const _2R
Name extension for 2D rate measured. 
 
Auxiliary class for map of PMT parameters. 
 
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory. 
 
Fit parameters for single PMT. 
 
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
 
void load(const std::string &file_name, JDetector &detector)
Load detector from input file. 
 
double getSurvivalProbability(const JPMTParameters ¶meters)
Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assu...
 
bool fixParameter(TF1 &f1, const JFitParameter_t ¶meter)
Fix fit parameter. 
 
no fit printf nominal n $STRING awk v X
 
Object reading from a list of files. 
 
Data structure for PMT parameters. 
 
do set_variable DETECTOR_TXT $WORKDIR detector
 
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
 
Fit parameters for two-fold coincidence rate due to K40. 
 
double QE
relative quantum efficiency 
 
static const char *const _1D
Name extension for 1D time offset. 
 
#define DEBUG(A)
Message macros.