64{
   68 
   70 
   72 
   73  string       inputFile;
   75  string       detectorFile;
   76  string       pmtFile;
   78  bool         reverse;
   79  bool         overwriteDetector;
   81  bool         fitAngle;
   82  bool         fitNoise;
   83  bool         fitModel;
   84  bool         fixQE;
   88 
   89  try { 
   90 
   92 
  103 
  104    JParser<> zap(
"Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
 
  105    
  107    zap[
'f'] = 
make_field(inputFile,         
"input file (output from JMergeCalibrateK40).");
 
  109    zap[
'a'] = 
make_field(detectorFile,      
"detector file.");
 
  110    zap[
'P'] = 
make_field(pmtFile,           
"specify PMT file name that can be given to JTriggerEfficiency.") = 
"";
 
  112                          "Fix time offset(s) of PMT(s) of certain module(s), e.g."
  113                          "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
  114                          "\nSame PMT offset can be fixed for all optical modules, e.g."
  115                          "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") = 
JPARSER::initialised();
 
  116    zap[
'r'] = 
make_field(reverse,           
"reverse TDC constraints due to option -! <TDC>.");
 
  117    zap[
'A'] = 
make_field(overwriteDetector, 
"overwrite detector file provided through '-a' with fitted time offsets.");
 
  118    zap[
'w'] = 
make_field(writeFits,         
"write fit results to ROOT file; -ww also write fitted TTS to PMT file.");
 
  119    zap[
'D'] = 
make_field(fitAngle,          
"fit angular distribution; fix normalisation.");
 
  120    zap[
'B'] = 
make_field(fitNoise,          
"fit background.");
 
  121    zap[
'M'] = 
make_field(fitModel,          
"fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
 
  122    zap[
'Q'] = 
make_field(fixQE,             
"fix PMT QEs = 1.0.");
 
  126 
  127    zap(argc, argv);
  128  }
  129  catch(const exception &error) {
  130    FATAL(error.what() << endl);
 
  131  }
  132 
  133 
  134  if ((fitModel ? 1 : 0)  +
  135      (fitAngle ? 1 : 0)  +
  136      (fitNoise ? 1 : 0)  +
  137      (fixQE    ? 1 : 0)  >  1) {
  138    FATAL(
"Use either option -M, -D, -B or -Q" << endl);
 
  139  } 
  140 
  146 
  147  if (reverse) {
  149  }
  150 
  151  for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
  152    DEBUG(
"PMT " << setw(10) << i->first << 
' ' << setw(2) << i->second << 
" constrain t0." << endl);
 
  153  }
  154 
  155  try {
  157  }
  158  catch(const exception &error) {
  159    FATAL(error.what() << endl);
 
  160  }
  161 
  163 
  164  try {
  166  }
  169  }
  170 
  171 
  173 
  174  if (pmtFile != "") {
  175    try {
  176      parameters.
load(pmtFile.c_str());
 
  177    }
  178    catch(const exception& error) {}
  179  }
  180  
  182 
  183 
  184  TFile* in = TFile::Open(inputFile.c_str(), "exist");
  185 
  186  if (in == NULL || !in->IsOpen()) {
  187    FATAL(
"File: " << inputFile << 
" not opened." << endl);
 
  188  }
  189 
  190 
  192 
  193  TH1D h0("chi2", NULL, 500,  0.0,   5.0);
  194  TH1D hn("hn",   NULL, 501, -0.5, 500.0);
  195  TH1D hr("rate", NULL, 500,  0.0,  25.0);
  196  TH1D h1("p1",   NULL, 500, -5.0,  +5.0);
  197  TH1D h2("p2",   NULL, 500, -5.0,  +5.0);
  198  TH1D h3("p3",   NULL, 500, -5.0,  +5.0);
  199  TH1D h4("p4",   NULL, 500, -5.0,  +5.0);
  200  TH1D hc("cc",   NULL, 500, -0.1,  +0.1);
  201 
  204 
  205  TH2D H2("detector", NULL,
  206          string.size() + 0, -0.5, string.size() - 0.5,
  208 
  209  for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
  210    H2.GetXaxis()->SetBinLabel(i, 
MAKE_CSTRING(
string.at(i-1)));
 
  211  }
  212  for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
  214  }
  215 
  216  TH2D* HN = (TH2D*) H2.Clone("iterations");
  217 
  219 
  220  for (JDetector::iterator module = 
detector.begin(); module != 
detector.end(); ++module) {
 
  221 
  222    if (module->getFloor() == 0) { 
  223      continue;
  224    }
  225 
  227 
  228    NOTICE(
"Module " << setw(10) << module->getID() << 
' ' << 
getLabel(module->getLocation()) << 
" !" << 
distance(range.first, range.second) << endl);
 
  229 
  231 
  232    if (h2d == NULL || h2d->GetEntries() == 0) {
  233      
  234      NOTICE(
"No data for module " << module->getID() << 
" -> set QEs to 0." << endl);
 
  235 
  238      }
  239 
  240      continue;
  241    }
  242 
  244 
  246 
  249 
  250    for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
  251 
  253      
  255 
  256      double V = 0.0;                                                  
  257      double W = 0.0;                                                  
  258 
  259      for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
  260 
  261        const double x = h2d->GetXaxis()->GetBinCenter(ix);
 
  262        const double y = h2d->GetYaxis()->GetBinCenter(iy);
 
  263 
  264        if (X(x) && Y(y)) {
  265 
  266          double value = h2d->GetBinContent(ix,iy);
  267          double error = h2d->GetBinError  (ix,iy);
  268 
  269          buffer.push_back(
rate_type(y, value, error));
 
  270 
  271          double width = h2d->GetYaxis()->GetBinWidth(iy);
  272 
  273          value *= width;
  274          error *= width;
  275 
  276          V += value;
  277          W += error * error;
  278        }
  279      }
  280 
  281      W = sqrt(W);
  282 
  283      if (V  <=  0.0             - STDEV*W) {
  284        count[0][
pair.first]  += 1;
 
  285        count[0][
pair.second] += 1;
 
  286      }
  287 
  288      if (V  <=  MINIMAL_RATE_HZ + STDEV*W) {
  289        count[1][
pair.first]  += 1;
 
  290        count[1][
pair.second] += 1;
 
  291      }
  292    }
  293 
  295 
  296      if (count[0][pmt] >= MAXIMAL_COUNTS) {                           
  297 
  298        WARNING(
"PMT " << setw(10) << module->getID() << 
'.' << 
FILL(2,
'0') << pmt << 
FILL() << 
" some rates negative -> fit background" << endl);
 
  299 
  300        if (fit.value.parameters[pmt].status) {
  301          model.parameters[pmt].bg.set();
 
  302        }
  303      }
  304      
  305      if (count[1][pmt] == NUMBER_OF_PMTS) {                           
  306 
  307        WARNING(
"PMT " << setw(10) << module->getID() << 
'.' << 
FILL(2,
'0') << pmt << 
FILL() << 
" all rates to low -> disable" << endl);
 
  308 
  309        model.parameters[pmt].disable();
 
  310      }
  311    }
  312 
  313    DEBUG(
"Start value:" << endl << 
model << endl);
 
  314 
  315    try {
  316 
  318 
  320 
  322 
  323        ERROR(
"Fit result " << setw(10) << module->getID() << 
" NDF  " << setw(5) << 
result.ndf << 
" -> skip" << endl);
 
  324 
  325        continue;
  326      }
  327 
  328      bool refit  = false;
  329 
  331 
  332        if (fit.value.parameters[pmt].status) {
  333 
  334          if (fit.value.parameters[pmt].QE() <= QE_MIN  ||
  335              fit.value.parameters[pmt].QE() <= 0.0 + STDEV * fit.error.parameters[pmt].QE()) {
  336 
  337            WARNING(
"PMT " << setw(10) << module->getID() << 
'.' << 
FILL(2,
'0') << pmt << 
FILL() << 
' ' 
  338                    << "QE = "
  339                    << 
FIXED(5,3) << fit.value.parameters[pmt].QE() << 
" +/- "  
  340                    << 
FIXED(5,3) << fit.error.parameters[pmt].QE() << 
" " 
  341                    << " -> disable" << (!refit ? " and refit" : "") << endl);
  342 
  343            fit.value.parameters[pmt].disable();
  344 
  345            refit = true;
  346          }
  347 
  348          if (fit.value.parameters[pmt].t0.atLimit(T0_NS)) {
  349 
  350            WARNING(
"PMT " << setw(10) << module->getID() << 
'.' << 
FILL(2,
'0') << pmt << 
FILL() << 
' ' 
  351                    << "t0 at limit "
  352                    << 
FIXED(5,3) << fit.value.parameters[pmt].t0() << 
" +/- "  
  353                    << 
FIXED(5,3) << fit.error.parameters[pmt].t0() << 
" " 
  354                    << " -> disable" << (!refit ? " and refit" : "") << endl);
  355 
  356            fit.value.parameters[pmt].disable();
  357 
  358            refit = true;
  359          }
  360        }
  361      }
  362 
  363      if (refit) {
  364 
  366          if (fit.value.parameters[pmt].status) {
  368          }
  369        }
  370 
  371        refit  = false;
  373      }
  374 
  375      NOTICE(
"Fit result " << setw(10) << module->getID() << 
" chi2 / NDF  " << 
FIXED(10,2) << 
result.chi2 << 
" / " << setw(5) << 
result.ndf << 
' ' << setw(5) << fit.numberOfIterations << endl);
 
  376 
  378 
  379      if (writeFits) {
  380 
  382        hn.Fill(fit.numberOfIterations);
  383        hr.Fill(fit.value.R );
  384        h1.Fill(fit.value.p1);
  385        h2.Fill(fit.value.p2);
  386        h3.Fill(fit.value.p3);
  387        h4.Fill(fit.value.p4);
  388        hc.Fill(fit.value.cc);
  389 
  390        TH1D h1t(
MAKE_CSTRING(module->getID() << 
".1t0"),  NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
 
  391        TH1D h1s(
MAKE_CSTRING(module->getID() << 
".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
 
  392        TH1D h1q(
MAKE_CSTRING(module->getID() << 
".1QE"),  NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
 
  393 
  395          h1t.SetBinContent(pmt + 1, fit.value.parameters[pmt].t0 ());
  396          h1t.SetBinError  (pmt + 1, fit.error.parameters[pmt].t0 ()  + numeric_limits<double>::epsilon());
  397          h1s.SetBinContent(pmt + 1, fit.value.parameters[pmt].TTS());
  398          h1s.SetBinError  (pmt + 1, fit.error.parameters[pmt].TTS()  + numeric_limits<double>::epsilon());
  399          h1q.SetBinContent(pmt + 1, fit.value.parameters[pmt].QE ());
  400          h1q.SetBinError  (pmt + 1, fit.error.parameters[pmt].QE ()  + numeric_limits<double>::epsilon());
  401        }
  402 
  403        out << h1t << h1s << h1q;
  404 
  405        for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
  406 
  408 
  409          for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
  410          
  411            const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
  412              
  413            h2d->SetBinContent(ix, iy, fit.value.getValue(
pair, dt_ns));
 
  414            h2d->SetBinError  (ix, iy, 0.0);
  415          }
  416        }
  417 
  419        h2d->Write();
  420 
  421        const double x = 
string.getIndex(module->getString());
 
  422        const double y = 
module->getFloor();
 
  423        
  425        HN->Fill(x, y, fit.numberOfIterations);
  426      }
  427 
  428      const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
  429 
  431      
  433 
  435            
  436        if (P > 0.0) 
  437          data.QE = fit.value.parameters[pmt].QE / P;
 
  438        else
  440 
  441        if (writeFits > 1) {
  442          data.TTS_ns = fit.value.parameters[pmt].TTS();
 
  443        }
  444        
  445        module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0.get() - t0);
  446      }
  447    }  
  448    catch(const exception& error) {
  449 
  450      ERROR(
"Module " << setw(10) << module->getID() << 
' ' << error.what() << 
" -> set QEs to 0." << endl);
 
  451 
  454      }
  455    }
  456  }
  457 
  458 
  460 
  461  {
  463    JSTDObjectWriter  <JMeta> writer(meta);
  464 
  465    writer << reader;
  466  }
  467 
  468  for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
  471  }
  472 
  473  if (overwriteDetector) {
  474 
  475    NOTICE(
"Store calibration data on file " << detectorFile << endl);
 
  476 
  478  }
  479 
  480  if (pmtFile != "") {
  481    parameters.
store(pmtFile.c_str());
 
  482  }
  483 
  484 
  485  for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
  487  }
  488 
  491  }
  492 
  493  if (writeFits) {
  494    out << h0 << hn << hr << h1 << h2 << h3 << h4 << hc << H2 << *HN;
  495  }
  496 
  497  out.Close();
  498 
  499  return 0;
  500}
#define DEBUG(A)
Message macros.
 
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
 
#define MAKE_CSTRING(A)
Make C-string.
 
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
 
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
 
Auxiliary class for map of PMT parameters.
 
Data structure for PMT parameters.
 
Utility class to parse parameter values.
 
Utility class to parse command line options.
 
Object reading from a list of files.
 
static const char *const _2F
Name extension for 2F rate fitted.
 
@ FIT_PMTS_QE_FIXED_t
fit parameters of PMTs with QE fixed
 
@ FIT_PMTS_AND_ANGULAR_DEPENDENCE_t
fit parameters of PMTs and angular dependence of K40 rate
 
@ FIT_MODEL_t
fit parameters of K40 rate and TTSs of PMTs
 
@ FIT_PMTS_AND_BACKGROUND_t
fit parameters of PMTs and background
 
@ FIT_PMTS_t
fit parameters of PMTs
 
static const char *const _2R
Name extension for 2D rate measured.
 
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
 
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
 
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
 
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
 
double getSurvivalProbability(const JPMTParameters ¶meters)
Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assu...
 
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
 
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
 
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
 
KM3NeT DAQ data structures and auxiliaries.
 
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
 
Auxiliary data structure for sequence of same character.
 
Auxiliary data structure for floating point format specification.
 
Type definition of range.
 
Model for fit to acoustics data.
 
Fit parameters for two-fold coincidence rate due to K40.
 
static const JK40Parameters & getInstance()
Get default values.
 
static const JPMTParameters_t & getInstance()
Get default values.
 
Auxiliary class for TDC constraints.
 
range_type equal_range(const int id)
Get range of constraints for given module.
 
void reverse()
Reverse constraints.
 
bool is_valid(const bool option=false) const
Check validity of TDC constrants.
 
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
 
Data structure for measured coincidence rate of pair of PMTs.
 
Router for mapping of string identifier to index.
 
void store(const char *file_name) const
Store to output file.
 
void load(const char *file_name)
Load from input file.
 
Empty structure for specification of parser element that is initialised (i.e. does not require input)...