74   typedef JRange<double>      JRange_t;
 
   78   const double      epsilon = 1.0e-10;    
 
   79   JFitK40Parameters fitk40;               
 
   87   bool      overwriteDetector;
 
  111     JParser<> zap(
"Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
 
  114     zap[
'f'] = 
make_field(inputFile,         
"input file (output from JMergeCalibrateK40).");
 
  116     zap[
'a'] = 
make_field(detectorFile,      
"detector file.");
 
  117     zap[
'P'] = 
make_field(pmtFile,           
"specify PMT file name that can be given to JTriggerEfficiency.") = 
"";
 
  119                           "Fix time offset(s) of PMT(s) of certain module(s), e.g." 
  120                           "\n-! \"808969848 0 808982077 23\" will fix offset of PMT 0 of module 808969848 and of PMT 23 of module 808982077." 
  121                           "\nSame PMT offset can be fixed for all optical modules, e.g." 
  122                           "\n-! \"-1 0 -1 23\" will fix offsets of PMTs 0 and 23 for all optical modules.")
 
  124     zap[
'r'] = 
make_field(reverse,           
"reverse TDC constraints due to option -! <TDC>.");
 
  125     zap[
'A'] = 
make_field(overwriteDetector, 
"overwrite detector file provided through '-a' with correct time offsets.");
 
  126     zap[
'w'] = 
make_field(writeFits,         
"write fit results.");
 
  127     zap[
'D'] = 
make_field(fitAngDist,        
"fit angular distribution; fix normalisation.");
 
  128     zap[
'M'] = 
make_field(fitModel,          
"fit angular distribution; fix PMT QEs = 1.0.");
 
  129     zap[
'E'] = 
make_field(mu,                
"expectation value for npe given two-fold coincidence")           = 0.0;
 
  130     zap[
'x'] = 
make_field(X,                 
"ROOT fit range (PMT pairs).")                                    = JRange_t();
 
  131     zap[
'y'] = 
make_field(Y,                 
"ROOT fit range (time residual).")                                = JRange_t();
 
  132     zap[
'O'] = 
make_field(option,            
"ROOT fit option, see TH1::Fit.")                                 = 
"";
 
  137   catch(
const exception &error) {
 
  138     FATAL(error.what() << endl);
 
  142   ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(1000000);
 
  151     load(detectorFile, detector);
 
  153   catch(
const JException& error) {
 
  158   JPMTParametersMap parameters;
 
  162       parameters.load(pmtFile.c_str());
 
  164     catch(
const JException& error) {
 
  169   parameters.comment.add(JMeta(argc, argv));
 
  176     const range_type range = TDC.equal_range(-1);       
 
  178     if (range.first != range.second) {
 
  180       const JTDC_t buffer(range.first, range.second);   
 
  182       TDC.erase(range.first, range.second);             
 
  184       for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
 
  186         for (JTDC_t::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
 
  192           if (i->second == -1) {
 
  195               TDC.insert(make_pair(module->getID(), pmt));
 
  200             TDC.insert(make_pair(module->getID(), i->second));
 
  215     for (JTDC_t::const_iterator __p = TDC.begin(); __p != TDC.end(); ) {
 
  217       JTDC_t::const_iterator __q = __p;
 
  219       for ( ; __q != TDC.end() && __q->first == __p->first; ++__q) {}
 
  223         JTDC_t::const_iterator __i = __p;
 
  225         for ( ; __i != __q && __i->second != i; ++__i) {}
 
  228           buffer.insert(std::make_pair(__p->first,i));
 
  238   for (JTDC_t::const_iterator tdc = TDC.begin(); tdc != TDC.end(); ++tdc) {
 
  239     DEBUG(
"PMT " << setw(10) << tdc->first << 
' ' << setw(2) << tdc->second << 
" constrain t0." << endl);
 
  242   for (JTDC_t::const_iterator tdc = TDC.begin(); tdc != TDC.end(); ++tdc) {
 
  243     if (tdc->first < 0) {
 
  244       FATAL(
"Illegal module identifier: " << tdc->first << endl);
 
  247       FATAL(
"Illegal TDC (= readout channel) identifier: " << tdc->second << 
" {0, .., " << 
NUMBER_OF_PMTS - 1 << 
"}" << endl);
 
  252   if (option.find(
'O') == string::npos) { option += 
'0'; }
 
  253   if (option.find(
'S') == string::npos) { option += 
'S'; }
 
  257   TFile* in = TFile::Open(inputFile.c_str(), 
"exist");
 
  259   if (in == NULL || !in->IsOpen()) {
 
  260     FATAL(
"File: " << inputFile << 
" not opened." << endl);
 
  265   TH1D hc(
"chi2", NULL, 500, 0.0, 5.0);
 
  267   TH1D 
p1(
"p1", NULL, 501, -5.0, +5.0);
 
  268   TH1D p2(
"p2", NULL, 501, -5.0, +5.0);
 
  269   TH1D p3(
"p3", NULL, 501, -5.0, +5.0);
 
  270   TH1D p4(
"p4", NULL, 501, -5.0, +5.0);
 
  271   TH1D bg(
"bg", NULL, 501, -0.1, +0.1);
 
  272   TH1D cc(
"cc", NULL, 501, -0.1, +0.1);
 
  281   for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
 
  283     NOTICE(
"Module " << setw(10) << module->getID() << endl);
 
  291     if (h2s == NULL || h2s->GetEntries() == 0) {
 
  293       NOTICE(
"No data for module " << module->getID() << 
"; set corresponding QEs to 0." << endl);
 
  296         parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
 
  306     const range_type range = TDC.equal_range(module->getID());
 
  313                 h2s->GetXaxis()->GetXmin(),
 
  314                 h2s->GetXaxis()->GetXmax(),
 
  315                 h2s->GetYaxis()->GetXmin(),
 
  316                 h2s->GetYaxis()->GetXmax(),
 
  317                 range.first == range.second);
 
  319     fit.setModelParameters(fitk40.getModelParameters());
 
  321     for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
 
  322       fixParameter(fit, fit.getModelParameter(i->second, &JPMTParameters_t::t0));
 
  327         fixParameter(fit, fit.getModelParameter(pmt, &JPMTParameters_t::QE));
 
  332       fixParameter(fit, fit.getModelParameter(&JFitK40::Rate_Hz));
 
  348     for (
int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
 
  350       Double_t value = 0.0;
 
  351       Double_t error = 0.0;
 
  353       for (
int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
 
  354         value += h2s->GetBinContent(ix,iy);
 
  355         error += h2s->GetBinError(ix,iy) * h2s->GetBinError(ix,iy);
 
  360       if (value  <  MINIMAL_RATE_HZ + STDEV*error) {
 
  364         buffer[
pair.first]  += 1;
 
  365         buffer[
pair.second] += 1;
 
  375           WARNING(
"PMT " << setw(10) << module->getID() << 
' ' << setw(2) << pmt << 
" too low a rate; switch off." << endl);
 
  381     catch(
const JException& error) {
 
  383       NOTICE(
"Skip fit for module " << module->getID() << 
"; set corresponding QEs to 0." << endl);
 
  386         parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
 
  394     if (X != JRange_t()) {
 
  395       h2s->GetXaxis()->SetRangeUser(max(X.getLowerLimit(), h2s->GetXaxis()->GetXmin()),
 
  396                                     min(X.getUpperLimit(), h2s->GetXaxis()->GetXmax()));
 
  398     if (Y != JRange_t()) {
 
  399       h2s->GetYaxis()->SetRangeUser(max(Y.getLowerLimit(), h2s->GetYaxis()->GetXmin()),
 
  400                                     min(Y.getUpperLimit(), h2s->GetYaxis()->GetXmax()));
 
  404     TFitResultPtr 
result = fit(*h2s, option);
 
  411       const JFitK40Parameters values(fit.GetParameters());
 
  412       const JFitK40Parameters errors(fit.GetParErrors());
 
  416         if (!
result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE))  &&  !
is_valid(values.getQE(pmt), errors.getQE(pmt))) {
 
  418           WARNING(
"PMT " << setw(10) << module->getID() << 
' ' << setw(2) << pmt << 
' ' 
  420                   << 
FIXED(5,2) << values.getQE(pmt) << 
" +/- "  
  421                   << 
FIXED(5,2) << errors.getQE(pmt) 
 
  422                   << 
" compatible with zero -> set QE = 0 and t0 = 0; -> refit." << endl);
 
  431     catch(
const JException& error) {
 
  433       NOTICE(
"Skip fit for module " << module->getID() << 
"; set corresponding QEs to 0." << endl);
 
  436         parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
 
  444       result = fit(*h2s, option);
 
  448     const JFitK40Parameters values(fit.GetParameters());
 
  449     const JFitK40Parameters errors(fit.GetParErrors());
 
  458              h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin (), h2s->GetXaxis()->GetXmax (),
 
  459              h2s->GetYaxis()->GetNbins(), h2s->GetYaxis()->GetXmin (), h2s->GetYaxis()->GetXmax ());
 
  461     for (
int ix = 1; ix <= h2f.GetXaxis()->GetNbins(); ++ix) {
 
  462       for (
int iy = 1; iy <= h2f.GetYaxis()->GetNbins(); ++iy) {
 
  464         const Double_t x = h2f.GetXaxis()->GetBinCenter(ix);
 
  465         const Double_t y = h2f.GetYaxis()->GetBinCenter(iy);
 
  467         h2f.SetBinContent(ix, iy, fit.Eval(x,y));
 
  468         h2f.SetBinError  (ix, iy, 0.0);
 
  475       cout << 
"Fit result chi2 / NDF " << 
result->Chi2() << 
' ' << 
result->Ndf() << 
' ' << (
result->IsValid() ? 
"" : 
"failed") << endl;
 
  477       cout << 
"Rate_Hz= " << 
FIXED(8,4) << values.Rate_Hz << (
result->IsParameterFixed(fit.getModelParameter(&JFitK40::Rate_Hz)) ? 
" *** fixed *** " : 
"") << endl;
 
  479       cout << 
"p1=      " << 
FIXED(8,4) << values.p1      << (
result->IsParameterFixed(fit.getModelParameter(&
JFitK40::p1))      ? 
" *** fixed *** " : 
"") << endl;
 
  480       cout << 
"p2=      " << 
FIXED(8,4) << values.p2      << (
result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2))      ? 
" *** fixed *** " : 
"") << endl;
 
  481       cout << 
"p3=      " << 
FIXED(8,4) << values.p3      << (
result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2))      ? 
" *** fixed *** " : 
"") << endl;
 
  482       cout << 
"p4=      " << 
FIXED(8,4) << values.p4      << (
result->IsParameterFixed(fit.getModelParameter(&JFitK40::p3))      ? 
" *** fixed *** " : 
"") << endl;
 
  484       cout << 
"Background (correlated):   " << 
FIXED(8,5) << values.bg << (
result->IsParameterFixed(fit.getModelParameter(&JFitK40::bg)) ? 
" *** fixed *** " : 
"") << endl;
 
  485       cout << 
"Background (uncorrelated): " << 
FIXED(8,5) << values.cc << (
result->IsParameterFixed(fit.getModelParameter(&JFitK40::cc)) ? 
" *** fixed *** " : 
"") << endl;
 
  487       cout << 
" PMT t0 [ns] TTS [ns]   R" << endl;
 
  494              << setw(2)    << pmt             << 
' ' 
  495              << 
FIXED(7,2) << fit.getT0 (pmt) << 
' ' 
  496              << 
FIXED(7,2) << fit.getTTS(pmt) << 
' ' 
  497              << 
FIXED(7,3) << fit.getQE (pmt);
 
  499         cout << (
result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE))  ? 
" *** fixed QE  *** " : 
"");
 
  500         cout << (
result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::TTS)) ? 
" *** fixed TTS *** " : 
"");
 
  501         cout << (
result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::t0))  ? 
" *** fixed t0  *** " : 
"");
 
  502         cout << (fit.is_disabled  (pmt) ? 
" (disabled)" : 
"");
 
  503         cout << (fit.is_average_t0(pmt) ? 
" (averaged)" : 
"");
 
  506         Q[0].put(fit.getT0 (pmt));
 
  507         Q[1].put(fit.getTTS(pmt));
 
  508         Q[2].put(fit.getQE (pmt));
 
  512            << 
FIXED(7,2) << Q[0].getMean()  << 
' ' 
  513            << 
FIXED(7,2) << Q[1].getMean()  << 
' ' 
  514            << 
FIXED(7,3) << Q[2].getMean()  << endl;
 
  516            << 
FIXED(7,2) << Q[0].getSTDev() << 
' ' 
  517            << 
FIXED(7,2) << Q[1].getSTDev() << 
' ' 
  518            << 
FIXED(7,3) << Q[2].getSTDev() << endl;
 
  545         h1t.SetBinContent(pmt + 1, values.getT0 (pmt));
 
  546         h1s.SetBinContent(pmt + 1, values.getTTS(pmt));
 
  547         h1q.SetBinContent(pmt + 1, values.getQE (pmt));
 
  549         h1t.SetBinError  (pmt + 1, max(errors.getT0 (pmt), epsilon));
 
  550         h1s.SetBinError  (pmt + 1, max(errors.getTTS(pmt), epsilon));
 
  551         h1q.SetBinError  (pmt + 1, max(errors.getQE (pmt), epsilon));
 
  554       out << h1t << h1s << h1q;
 
  562       parameters[JPMTIdentifier(module->getID(), pmt)].QE = values.getQE(pmt);
 
  564       if (overwriteDetector) {
 
  565         module->getPMT(pmt).addT0(fit.getT0(pmt));
 
  572     out << hc << 
p1 << p2 << p3 << p4 << bg << cc;
 
  578   if (overwriteDetector) {
 
  580     NOTICE(
"Store calibration data on file " << detectorFile << endl);
 
  582     detector.comment.add(JMeta(argc, argv));
 
  584     store(detectorFile, detector);
 
  594       parameters.convertHitProbabilityToQE(mu);
 
  596     catch(
const JException& error) {
 
  600     ofstream out(pmtFile.c_str());
 
  602     out << parameters << endl;