65{
   69 
   71 
   73 
   74  string       inputFile;
   76  string       detectorFile;
   77  string       pmtFile;
   79  bool         reverse;
   80  bool         overwriteDetector;
   82  bool         fitAngle;
   83  bool         fitNoise;
   84  bool         fitModel;
   85  bool         fixQE;
   89 
   90  try { 
   91 
   93 
  108 
  109    JParser<> zap(
"Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
 
  110    
  112    zap[
'f'] = 
make_field(inputFile,         
"input file (output from JMergeCalibrateK40).");
 
  114    zap[
'a'] = 
make_field(detectorFile,      
"detector file.");
 
  115    zap[
'P'] = 
make_field(pmtFile,           
"specify PMT file name that can be given to JTriggerEfficiency.") = 
"";
 
  117                          "Fix time offset(s) of PMT(s) of certain module(s), e.g."
  118                          "\n-! \"808969848 0 808982077 23\" will fix time offsets of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
  119                          "\nSame PMT offset can be fixed for all optical modules, e.g."
  120                          "\n-! \"-1 0 -1 23\" will fix time offsets of PMTs 0 and 23 of all optical modules.") = 
JPARSER::initialised();
 
  121    zap[
'r'] = 
make_field(reverse,           
"reverse TDC constraints due to option -! <TDC>.");
 
  122    zap[
'A'] = 
make_field(overwriteDetector, 
"overwrite detector file provided through '-a' with fitted time offsets.");
 
  123    zap[
'w'] = 
make_field(writeFits,         
"write fit results to ROOT file; -ww also write fitted TTS to PMT file.");
 
  124    zap[
'D'] = 
make_field(fitAngle,          
"fit angular distribution; fix normalisation.");
 
  125    zap[
'B'] = 
make_field(fitNoise,          
"fit background.");
 
  126    zap[
'M'] = 
make_field(fitModel,          
"fit angular distribution as well as normalisation; fix PMT QEs = 1.0.");
 
  127    zap[
'Q'] = 
make_field(fixQE,             
"fix PMT QEs = 1.0.");
 
  131 
  132    zap(argc, argv);
  133  }
  134  catch(const exception &error) {
  135    FATAL(error.what() << endl);
 
  136  }
  137 
  138 
  139  if ((fitModel ? 1 : 0)  +
  140      (fitAngle ? 1 : 0)  +
  141      (fitNoise ? 1 : 0)  +
  142      (fixQE    ? 1 : 0)  >  1) {
  143    FATAL(
"Use either option -M, -D, -B or -Q" << endl);
 
  144  } 
  145 
  151 
  152  if (reverse) {
  154  }
  155 
  156  for (JTDC_t::const_iterator i = TDC.begin(); i != TDC.end(); ++i) {
  157    DEBUG(
"PMT " << setw(10) << i->first << 
' ' << setw(2) << i->second << 
" constrain t0." << endl);
 
  158  }
  159 
  160  try {
  162  }
  163  catch(const exception &error) {
  164    FATAL(error.what() << endl);
 
  165  }
  166 
  168 
  169  try {
  171  }
  174  }
  175 
  176 
  178 
  179  if (pmtFile != "") {
  180    try {
  181      parameters.
load(pmtFile.c_str());
 
  182    }
  183    catch(const exception& error) {}
  184  }
  185  
  187 
  188 
  189  TFile* in = TFile::Open(inputFile.c_str(), "exist");
  190 
  191  if (in == NULL || !in->IsOpen()) {
  192    FATAL(
"File: " << inputFile << 
" not opened." << endl);
 
  193  }
  194 
  195 
  197 
  198  TH1D h0("chi2", NULL, 500,  0.0,   5.0);
  199  TH1D hn("hn",   NULL, 501, -0.5, 500.0);
  200  TH1D hr("rate", NULL, 500,  0.0,  25.0);
  201  TH1D h1("p1",   NULL, 500, -5.0,  +5.0);
  202  TH1D h2("p2",   NULL, 500, -5.0,  +5.0);
  203  TH1D h3("p3",   NULL, 500, -5.0,  +5.0);
  204  TH1D h4("p4",   NULL, 500, -5.0,  +5.0);
  205  TH1D hc("cc",   NULL, 500, -0.1,  +0.1);
  206  TH1D hb("bc",   NULL, 500, -0.1,  +0.1);
  207 
  210 
  211  TH2D H2("detector", NULL,
  212          string.size() + 0, -0.5, string.size() - 0.5,
  214 
  215  for (Int_t i = 1; i <= H2.GetXaxis()->GetNbins(); ++i) {
  216    H2.GetXaxis()->SetBinLabel(i, 
MAKE_CSTRING(
string.at(i-1)));
 
  217  }
  218  for (Int_t i = 1; i <= H2.GetYaxis()->GetNbins(); ++i) {
  220  }
  221 
  222  TH2D* HN = (TH2D*) H2.Clone("iterations");
  223 
  225 
  226  for (JDetector::iterator module = 
detector.begin(); module != 
detector.end(); ++module) {
 
  227 
  228    if (module->getFloor() == 0) { 
  229      continue;
  230    }
  231 
  233 
  234    NOTICE(
"Module " << setw(10) << module->getID() << 
' ' << 
getLabel(module->getLocation()) << 
" !" << 
distance(range.first, range.second) << endl);
 
  235 
  237 
  238    if (h2d == NULL || h2d->GetEntries() == 0) {
  239      
  240      NOTICE(
"No data for module " << module->getID() << 
" -> set QEs to 0." << endl);
 
  241 
  244      }
  245 
  246      continue;
  247    }
  248 
  250 
  252 
  255 
  256    for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
  257 
  259      
  261 
  262      double V = 0.0;                                                  
  263      double W = 0.0;                                                  
  264 
  265      for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
  266 
  267        const double x = h2d->GetXaxis()->GetBinCenter(ix);
 
  268        const double y = h2d->GetYaxis()->GetBinCenter(iy);
 
  269 
  270        if (X(x) && Y(y)) {
  271 
  272          double value = h2d->GetBinContent(ix,iy);
  273          double error = h2d->GetBinError  (ix,iy);
  274 
  275          buffer.push_back(
rate_type(y, value, error));
 
  276 
  277          double width = h2d->GetYaxis()->GetBinWidth(iy);
  278 
  279          value *= width;
  280          error *= width;
  281 
  282          V += value;
  283          W += error * error;
  284        }
  285      }
  286 
  287      W = sqrt(W);
  288 
  289      if (V  <=  0.0             - STDEV*W) {
  290        count[0][
pair.first]  += 1;
 
  291        count[0][
pair.second] += 1;
 
  292      }
  293 
  294      if (V  <=  MINIMAL_RATE_HZ + STDEV*W) {
  295        count[1][
pair.first]  += 1;
 
  296        count[1][
pair.second] += 1;
 
  297      }
  298    }
  299 
  301 
  302      if (count[0][pmt] >= MAXIMAL_COUNTS) {                           
  303 
  304        WARNING(
"PMT " << setw(10) << module->getID() << 
'.' << 
FILL(2,
'0') << pmt << 
FILL() << 
" some rates negative -> fit background" << endl);
 
  305 
  306        if (fit.value.parameters[pmt].status) {
  307          model.parameters[pmt].bg.set();
 
  308        }
  309      }
  310      
  311      if (count[1][pmt] == NUMBER_OF_PMTS) {                           
  312 
  313        WARNING(
"PMT " << setw(10) << module->getID() << 
'.' << 
FILL(2,
'0') << pmt << 
FILL() << 
" all rates to low -> disable" << endl);
 
  314 
  315        model.parameters[pmt].disable();
 
  316      }
  317    }
  318 
  319    DEBUG(
"Start value:" << endl << 
model << endl);
 
  320 
  321    try {
  322 
  324 
  326 
  328 
  329        ERROR(
"Fit result " << setw(10) << module->getID() << 
" NDF  " << setw(5) << 
result.ndf << 
" -> skip" << endl);
 
  330 
  331        continue;
  332      }
  333 
  334      bool refit  = false;
  335 
  337 
  338        if (fit.value.parameters[pmt].status) {
  339 
  340          if (fit.value.parameters[pmt].QE() <= QE_MIN  ||
  341              fit.value.parameters[pmt].QE() <= 0.0 + STDEV * fit.error.parameters[pmt].QE()) {
  342 
  343            WARNING(
"PMT " << setw(10) << module->getID() << 
'.' << 
FILL(2,
'0') << pmt << 
FILL() << 
' ' 
  344                    << "QE = "
  345                    << 
FIXED(5,3) << fit.value.parameters[pmt].QE() << 
" +/- "  
  346                    << 
FIXED(5,3) << fit.error.parameters[pmt].QE() << 
" " 
  347                    << " -> disable" << (!refit ? " and refit" : "") << endl);
  348 
  349            fit.value.parameters[pmt].disable();
  350 
  351            refit = true;
  352          }
  353 
  354          if (fit.value.parameters[pmt].t0.atLimit(T0_NS)) {
  355 
  356            WARNING(
"PMT " << setw(10) << module->getID() << 
'.' << 
FILL(2,
'0') << pmt << 
FILL() << 
' ' 
  357                    << "t0 at limit "
  358                    << 
FIXED(5,3) << fit.value.parameters[pmt].t0() << 
" +/- "  
  359                    << 
FIXED(5,3) << fit.error.parameters[pmt].t0() << 
" " 
  360                    << " -> disable" << (!refit ? " and refit" : "") << endl);
  361 
  362            fit.value.parameters[pmt].disable();
  363 
  364            refit = true;
  365          }
  366        }
  367      }
  368 
  369      if (refit) {
  370 
  372          if (fit.value.parameters[pmt].status) {
  374          }
  375        }
  376 
  377        refit  = false;
  379      }
  380 
  381      NOTICE(
"Fit result " << setw(10) << module->getID() << 
" chi2 / NDF  " << 
FIXED(10,2) << 
result.chi2 << 
" / " << setw(5) << 
result.ndf << 
' ' << setw(5) << fit.numberOfIterations << endl);
 
  382 
  385      }
  386 
  388 
  389      if (writeFits) {
  390 
  392        hn.Fill(fit.numberOfIterations);
  393        hr.Fill(fit.value.R );
  394        h1.Fill(fit.value.p1);
  395        h2.Fill(fit.value.p2);
  396        h3.Fill(fit.value.p3);
  397        h4.Fill(fit.value.p4);
  398        hc.Fill(fit.value.cc);
  399        hb.Fill(fit.value.bc);
  400 
  401        TH1D h1t(
MAKE_CSTRING(module->getID() << 
".1t0"),  NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
 
  402        TH1D h1s(
MAKE_CSTRING(module->getID() << 
".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
 
  403        TH1D h1q(
MAKE_CSTRING(module->getID() << 
".1QE"),  NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
 
  404 
  406          h1t.SetBinContent(pmt + 1, fit.value.parameters[pmt].t0 ());
  407          h1t.SetBinError  (pmt + 1, fit.error.parameters[pmt].t0 ()  + numeric_limits<double>::epsilon());
  408          h1s.SetBinContent(pmt + 1, fit.value.parameters[pmt].TTS());
  409          h1s.SetBinError  (pmt + 1, fit.error.parameters[pmt].TTS()  + numeric_limits<double>::epsilon());
  410          h1q.SetBinContent(pmt + 1, fit.value.parameters[pmt].QE ());
  411          h1q.SetBinError  (pmt + 1, fit.error.parameters[pmt].QE ()  + numeric_limits<double>::epsilon());
  412        }
  413 
  414        out << h1t << h1s << h1q;
  415 
  416        for (int ix = 1; ix <= h2d->GetXaxis()->GetNbins(); ++ix) {
  417 
  419 
  420          for (int iy = 1; iy <= h2d->GetYaxis()->GetNbins(); ++iy) {
  421          
  422            const double dt_ns = h2d->GetYaxis()->GetBinCenter(iy);
  423              
  424            h2d->SetBinContent(ix, iy, fit.value.getValue(
pair, dt_ns));
 
  425            h2d->SetBinError  (ix, iy, 0.0);
  426          }
  427        }
  428 
  430        h2d->Write();
  431 
  432        const double x = 
string.getIndex(module->getString());
 
  433        const double y = 
module->getFloor();
 
  434        
  436        HN->Fill(x, y, fit.numberOfIterations);
  437      }
  438 
  439      const double t0 = (fit.value.hasFixedTimeOffset() ? fit.value.getFixedTimeOffset() : 0.0);
  440 
  442      
  444 
  446            
  447        if (P > 0.0) 
  448          data.QE = fit.value.parameters[pmt].QE / P;
 
  449        else
  451 
  452        if (writeFits > 1) {
  453          data.TTS_ns = fit.value.parameters[pmt].TTS();
 
  454        }
  455        
  456        module->getPMT(pmt).addT0(fit.value.parameters[pmt].t0.get() - t0);
  457      }
  458    }  
  459    catch(const exception& error) {
  460 
  461      ERROR(
"Module " << setw(10) << module->getID() << 
' ' << error.what() << 
" -> set QEs to 0." << endl);
 
  462 
  465      }
  466    }
  467  }
  468 
  469 
  471 
  472  {
  474    JSTDObjectWriter  <JMeta> writer(meta);
  475 
  476    writer << reader;
  477  }
  478 
  479  for (vector<JMeta>::const_reverse_iterator i = meta.rbegin(); i != meta.rend(); ++i) {
  482  }
  483 
  484  if (overwriteDetector) {
  485 
  486    NOTICE(
"Store calibration data on file " << detectorFile << endl);
 
  487 
  489  }
  490 
  491  if (pmtFile != "") {
  492    parameters.
store(pmtFile.c_str());
 
  493  }
  494 
  495 
  496  for (vector<JMeta>::const_iterator i = meta.begin(); i != meta.end(); ++i) {
  498  }
  499 
  502  }
  503 
  504  if (writeFits) {
  505    out << h0 << hn << hr << h1 << h2 << h3 << h4 << hc << hb << H2 << *HN;
  506  }
  507 
  508  out.Close();
  509 
  510  return 0;
  511}
#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
 
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 double TEROSTAT_R1
scaling factor
 
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.
 
static double TEROSTAT_DZ
maximal PMT inclination
 
static double BELL_SHAPE
Bell shape.
 
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)...