JFitK40.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <map>
#include <cmath>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF2.h"
#include "TMath.h"
#include "TFitResult.h"
#include "Math/MinimizerOptions.h"

#include "JTools/JConstants.hh"
#include "JDAQ/JDAQ.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"

#include "JROOT/JRootToolkit.hh"
#include "JTools/JRange.hh"
#include "JTools/JQuantile.hh"
#include "JSupport/JMeta.hh"

#include "JCalibrate/JCalibrateK40.hh"
#include "JCalibrate/JFitK40.hh"

#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

namespace {

  double  STDEV            =  2.0;        //!< number of standard deviations
  double  QE_MIN_VALUE     =  0.01;       //!< minimal QE value (below, set to zero)
  double  QE_MAX_ERROR     =  0.5;        //!< maximal QE error (below, set to zero)
  double  MINIMAL_RATE_HZ  =  1.0e-2;     //!< minimal coincidence rate (bin-by-bin)

  /**
   * Check validity of QE measurement.
   *
   * The QE measurement is valid when:
   *  - fitted value is at least STDEV standard deviations larger than QE_MIN_VALUE; and
   *  - fitted error is smaller than QE_MAX_ERROR.
   *
   * \param  qe_value      QE value
   * \param  qe_error      QE error
   * \return               true if valid; else false
   */
  inline bool is_valid(const double qe_value, const double qe_error)
  {
    return (qe_value > QE_MIN_VALUE  +  STDEV * qe_error  &&
            qe_error < QE_MAX_ERROR);
  }
}


/**
 * \file
 *
 * Auxiliary program to fit PMT parameters from JMergeCalibrateK40.cc output.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JRange<double>      JRange_t;
  typedef multimap<int, int>  JTDC_t;
  typedef pair<JTDC_t::const_iterator, JTDC_t::const_iterator>  range_type;

  const double      epsilon = 1.0e-10;    // minimal error for ROOT plotting
  JFitK40Parameters fitk40;               // setting of internal fit parameters

  string    inputFile;
  string    outputFile;
  string    detectorFile;
  string    pmtFile;
  JTDC_t    TDC;
  bool      reverse;
  bool      overwriteDetector;
  bool      writeFits;
  bool      fitAngDist;
  bool      fitModel;
  double    mu;
  JRange_t  X;
  JRange_t  Y;
  string    option;
  int       debug;

  try { 

    JProperties properties;

    properties.insert(gmake_property(STDEV));
    properties.insert(gmake_property(QE_MIN_VALUE));
    properties.insert(gmake_property(QE_MAX_ERROR));
    properties.insert(gmake_property(MINIMAL_RATE_HZ));
    properties.insert(gmake_property(fitk40.Rate_Hz));
    properties.insert(gmake_property(fitk40.p1));
    properties.insert(gmake_property(fitk40.p2));
    properties.insert(gmake_property(fitk40.p3));
    properties.insert(gmake_property(fitk40.p4));
 
    JParser<> zap("Auxiliary program to fit PMT parameters from JMergeCalibrateK40 output.");
    
    zap['@'] = make_field(properties)        = JPARSER::initialised();
    zap['f'] = make_field(inputFile,         "input file (output from JMergeCalibrateK40).");
    zap['o'] = make_field(outputFile,        "output file.")                                                   = "fit.root";
    zap['a'] = make_field(detectorFile,      "detector file.");
    zap['P'] = make_field(pmtFile,           "specify PMT file name that can be given to JTriggerEfficiency.") = "";
    zap['!'] = make_field(TDC,
                          "Fix time offset(s) of PMT(s) of certain module(s), e.g."
                          "\n-! \"808969848 0 808982077 23\" will fix offset of PMT 0 of module 808969848 and of PMT 23 of module 808982077."
                          "\nSame PMT offset can be fixed for all optical modules, e.g."
                          "\n-! \"-1 0 -1 23\" will fix offsets of PMTs 0 and 23 for all optical modules.")
                                                                                                               = JPARSER::initialised();
    zap['r'] = make_field(reverse,           "reverse TDC constraints due to option -! <TDC>.");
    zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
    zap['w'] = make_field(writeFits,         "write fit results.");
    zap['D'] = make_field(fitAngDist,        "fit angular distribution; fix normalisation.");
    zap['M'] = make_field(fitModel,          "fit angular distribution; fix PMT QEs = 1.0.");
    zap['E'] = make_field(mu,                "expectation value for npe given two-fold coincidence")           = 0.0;
    zap['x'] = make_field(X,                 "ROOT fit range (PMT pairs).")                                    = JRange_t();
    zap['y'] = make_field(Y,                 "ROOT fit range (time residual).")                                = JRange_t();
    zap['O'] = make_field(option,            "ROOT fit option, see TH1::Fit.")                                 = "";
    zap['d'] = make_field(debug,             "debug.")                                                         = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(1000000);

  //----------------------------------------------------------
  // load detector file and PMT parameters
  //----------------------------------------------------------

  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }


  JPMTParametersMap parameters;

  if (pmtFile != "") {
    try {
      parameters.load(pmtFile.c_str());
    }
    catch(const JException& error) {
      //ERROR(error.what());
    }
  }
  
  parameters.comment.add(JMeta(argc, argv));

  {
    //----------------------------------------------------------
    // module identifier -1 corresponds to any module
    //----------------------------------------------------------
        
    const range_type range = TDC.equal_range(-1);       // module wild card

    if (range.first != range.second) {

      const JTDC_t buffer(range.first, range.second);   // temporary copy

      TDC.erase(range.first, range.second);             // remove entries
      
      for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {

        for (JTDC_t::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {

          //----------------------------------------------------------
          // PMT identifier -1 corresponds to any PMT in module
          //----------------------------------------------------------

          if (i->second == -1) {

            for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
              TDC.insert(make_pair(module->getID(), pmt));
            }

          } else { 

            TDC.insert(make_pair(module->getID(), i->second));
          }
        }
      }
    }
  }

  if (reverse) {

    //----------------------------------------------------------
    // reverse TDC constraints
    //----------------------------------------------------------

    JTDC_t buffer;

    for (JTDC_t::const_iterator __p = TDC.begin(); __p != TDC.end(); ) {
      
      JTDC_t::const_iterator __q = __p;
      
      for ( ; __q != TDC.end() && __q->first == __p->first; ++__q) {}
      
      for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
        
        JTDC_t::const_iterator __i = __p;
        
        for ( ; __i != __q && __i->second != i; ++__i) {}

        if (__i == __q) {
          buffer.insert(std::make_pair(__p->first,i));
        }
      }
      
      __p = __q;
    }
    
    TDC.swap(buffer);
  }

  for (JTDC_t::const_iterator tdc = TDC.begin(); tdc != TDC.end(); ++tdc) {
    DEBUG("PMT " << setw(10) << tdc->first << ' ' << setw(2) << tdc->second << " constrain t0." << endl);
  }

  for (JTDC_t::const_iterator tdc = TDC.begin(); tdc != TDC.end(); ++tdc) {
    if (tdc->first < 0) {
      FATAL("Illegal module identifier: " << tdc->first << endl);
    }
    if (tdc->second < 0 || tdc->second >= NUMBER_OF_PMTS) {
      FATAL("Illegal TDC (= readout channel) identifier: " << tdc->second << " {0, .., " << NUMBER_OF_PMTS - 1 << "}" << endl);
    }
  }


  if (option.find('O') == string::npos) { option += '0'; }
  if (option.find('S') == string::npos) { option += 'S'; }
  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }


  TFile* in = TFile::Open(inputFile.c_str(), "exist");

  if (in == NULL || !in->IsOpen()) {
    FATAL("File: " << inputFile << " not opened." << endl);
  }


  // Histograms with global fit results.
  TH1D hc("chi2", NULL, 500, 0.0, 5.0);

  TH1D p1("p1", NULL, 501, -5.0, +5.0);
  TH1D p2("p2", NULL, 501, -5.0, +5.0);
  TH1D p3("p3", NULL, 501, -5.0, +5.0);
  TH1D p4("p4", NULL, 501, -5.0, +5.0);
  TH1D bg("bg", NULL, 501, -0.1, +0.1);
  TH1D cc("cc", NULL, 501, -0.1, +0.1);


  TFile out(outputFile.c_str(), "recreate");

  //----------------------------------------------------------
  // loop over modules in detector file to fit histogram
  //----------------------------------------------------------

  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {

    NOTICE("Module " << setw(10) << module->getID() << endl);

    //----------------------------------------------------------
    // get histogram for this module
    //----------------------------------------------------------

    TH2D* h2s = (TH2D*) in->Get(MAKE_CSTRING(module->getID() << _2R));

    if (h2s == NULL || h2s->GetEntries() == 0) {
      
      NOTICE("No data for module " << module->getID() << "; set corresponding QEs to 0." << endl);

      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
      }

      continue;
    }

    //----------------------------------------------------------
    // get constraints
    //----------------------------------------------------------

    const range_type range = TDC.equal_range(module->getID());

    //----------------------------------------------------------
    // set-up fitting procedure
    //----------------------------------------------------------

    JFitK40 fit(*module,
                h2s->GetXaxis()->GetXmin(),
                h2s->GetXaxis()->GetXmax(),
                h2s->GetYaxis()->GetXmin(),
                h2s->GetYaxis()->GetXmax(),
                range.first == range.second);

    fit.setModelParameters(fitk40.getModelParameters());

    for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
      fixParameter(fit, fit.getModelParameter(i->second, &JPMTParameters_t::t0));
    }

    if (fitModel) {
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        fixParameter(fit, fit.getModelParameter(pmt, &JPMTParameters_t::QE));
      }

    } else {

      fixParameter(fit, fit.getModelParameter(&JFitK40::Rate_Hz));
      
      if (!fitAngDist) {
        fixParameter(fit, fit.getModelParameter(&JFitK40::p1));
        fixParameter(fit, fit.getModelParameter(&JFitK40::p2));
        fixParameter(fit, fit.getModelParameter(&JFitK40::p3));
        fixParameter(fit, fit.getModelParameter(&JFitK40::p4));
      }
    }

    //----------------------------------------------------------
    // check validity of data
    //----------------------------------------------------------

    vector<int> buffer(NUMBER_OF_PMTS, 1);

    for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {

      Double_t value = 0.0;
      Double_t error = 0.0;

      for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
        value += h2s->GetBinContent(ix,iy);
        error += h2s->GetBinError(ix,iy) * h2s->GetBinError(ix,iy);
      }

      error = sqrt(error);

      if (value  <  MINIMAL_RATE_HZ + STDEV*error) {

        const JFitK40::pair_type pair = fit.getPair(ix-1);

        buffer[pair.first]  += 1;
        buffer[pair.second] += 1;
      }
    }

    try {

      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
      
        if (buffer[pmt] == NUMBER_OF_PMTS) {

          WARNING("PMT " << setw(10) << module->getID() << ' ' << setw(2) << pmt << " too low a rate; switch off." << endl);

          fit.disablePMT(pmt);
        }
      }
    }
    catch(const JException& error) {

      NOTICE("Skip fit for module " << module->getID() << "; set corresponding QEs to 0." << endl);

      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
      }

      continue;
    }

    // constrain fit to specified range(s)

    if (X != JRange_t()) {
      h2s->GetXaxis()->SetRangeUser(max(X.getLowerLimit(), h2s->GetXaxis()->GetXmin()),
                                    min(X.getUpperLimit(), h2s->GetXaxis()->GetXmax()));
    }
    if (Y != JRange_t()) {
      h2s->GetYaxis()->SetRangeUser(max(Y.getLowerLimit(), h2s->GetYaxis()->GetXmin()),
                                    min(Y.getUpperLimit(), h2s->GetYaxis()->GetXmax()));
    }


    TFitResultPtr result = fit(*h2s, option);


    bool refit = false;

    try {

      const JFitK40Parameters values(fit.GetParameters());
      const JFitK40Parameters errors(fit.GetParErrors());
    
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {

        if (!result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE))  &&  !is_valid(values.getQE(pmt), errors.getQE(pmt))) {

          WARNING("PMT " << setw(10) << module->getID() << ' ' << setw(2) << pmt << ' '
                  << "measured QE " 
                  << FIXED(5,2) << values.getQE(pmt) << " +/- " 
                  << FIXED(5,2) << errors.getQE(pmt) 
                  << " compatible with zero -> set QE = 0 and t0 = 0; -> refit." << endl);


          fit.disablePMT(pmt);
          
          refit = true;
        }
      }
    }
    catch(const JException& error) {
      
      NOTICE("Skip fit for module " << module->getID() << "; set corresponding QEs to 0." << endl);

      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        parameters[JPMTIdentifier(module->getID(), pmt)].QE = 0.0;
      }
      
      continue;
    }

    
    if (refit) {
      result = fit(*h2s, option);
    }


    const JFitK40Parameters values(fit.GetParameters());
    const JFitK40Parameters errors(fit.GetParErrors());

    hc.Fill(TMath::Log10(result->Chi2() / result->Ndf()));

    // To force drawing of fit functions bin-by-bin compatible to the corresponding histogram data,
    // the function values will be stored as histograms, rather than as associated function objects.    
    
    TH2D h2f(MAKE_CSTRING(module->getID() << _2F),
             NULL,
             h2s->GetXaxis()->GetNbins(), h2s->GetXaxis()->GetXmin (), h2s->GetXaxis()->GetXmax (),
             h2s->GetYaxis()->GetNbins(), h2s->GetYaxis()->GetXmin (), h2s->GetYaxis()->GetXmax ());
    
    for (int ix = 1; ix <= h2f.GetXaxis()->GetNbins(); ++ix) {
      for (int iy = 1; iy <= h2f.GetYaxis()->GetNbins(); ++iy) {

        const Double_t x = h2f.GetXaxis()->GetBinCenter(ix);
        const Double_t y = h2f.GetYaxis()->GetBinCenter(iy);

        h2f.SetBinContent(ix, iy, fit.Eval(x,y));
        h2f.SetBinError  (ix, iy, 0.0);
      }
    }
     

    if (debug >= JEEP::notice_t) {

      cout << "Fit result chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;

      cout << "Rate_Hz= " << FIXED(8,4) << values.Rate_Hz << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::Rate_Hz)) ? " *** fixed *** " : "") << endl;

      cout << "p1=      " << FIXED(8,4) << values.p1      << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p1))      ? " *** fixed *** " : "") << endl;
      cout << "p2=      " << FIXED(8,4) << values.p2      << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2))      ? " *** fixed *** " : "") << endl;
      cout << "p3=      " << FIXED(8,4) << values.p3      << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p2))      ? " *** fixed *** " : "") << endl;
      cout << "p4=      " << FIXED(8,4) << values.p4      << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::p3))      ? " *** fixed *** " : "") << endl;

      cout << "Background (correlated):   " << FIXED(8,5) << values.bg << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::bg)) ? " *** fixed *** " : "") << endl;
      cout << "Background (uncorrelated): " << FIXED(8,5) << values.cc << (result->IsParameterFixed(fit.getModelParameter(&JFitK40::cc)) ? " *** fixed *** " : "") << endl;

      cout << " PMT t0 [ns] TTS [ns]   R" << endl;

      JQuantile Q[3];
      
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {

        cout << ' ' 
             << setw(2)    << pmt             << ' '
             << FIXED(7,2) << fit.getT0 (pmt) << ' '
             << FIXED(7,2) << fit.getTTS(pmt) << ' '
             << FIXED(7,3) << fit.getQE (pmt);

        cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::QE))  ? " *** fixed QE  *** " : "");
        cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::TTS)) ? " *** fixed TTS *** " : "");
        cout << (result->IsParameterFixed(fit.getModelParameter(pmt, &JPMTParameters_t::t0))  ? " *** fixed t0  *** " : "");
        cout << (fit.is_disabled  (pmt) ? " (disabled)" : "");
        cout << (fit.is_average_t0(pmt) ? " (averaged)" : "");
        cout << endl;

        Q[0].put(fit.getT0 (pmt));
        Q[1].put(fit.getTTS(pmt));
        Q[2].put(fit.getQE (pmt));
      }

      cout << "<>  "
           << FIXED(7,2) << Q[0].getMean()  << ' '
           << FIXED(7,2) << Q[1].getMean()  << ' '
           << FIXED(7,3) << Q[2].getMean()  << endl;
      cout << "+/- "
           << FIXED(7,2) << Q[0].getSTDev() << ' '
           << FIXED(7,2) << Q[1].getSTDev() << ' '
           << FIXED(7,3) << Q[2].getSTDev() << endl;
    }

    h2s->Write();
    h2f .Write();

    if (writeFits) {

      // Histograms with fit results.

      p1.Fill(values.p1);
      p2.Fill(values.p2);
      p3.Fill(values.p3);
      p4.Fill(values.p4);
      bg.Fill(values.bg);
      cc.Fill(values.cc);

      TH1D h1t(MAKE_CSTRING(module->getID() << ".1t0"),  NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
      TH1D h1s(MAKE_CSTRING(module->getID() << ".1TTS"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
      TH1D h1q(MAKE_CSTRING(module->getID() << ".1QE"),  NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
      
      h1t.Sumw2();
      h1s.Sumw2();
      h1q.Sumw2();

      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {      

        h1t.SetBinContent(pmt + 1, values.getT0 (pmt));
        h1s.SetBinContent(pmt + 1, values.getTTS(pmt));
        h1q.SetBinContent(pmt + 1, values.getQE (pmt));

        h1t.SetBinError  (pmt + 1, max(errors.getT0 (pmt), epsilon));
        h1s.SetBinError  (pmt + 1, max(errors.getTTS(pmt), epsilon));
        h1q.SetBinError  (pmt + 1, max(errors.getQE (pmt), epsilon));
      }
      
      out << h1t << h1s << h1q;
    }


    // save output

    for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
      
      parameters[JPMTIdentifier(module->getID(), pmt)].QE = values.getQE(pmt);
        
      if (overwriteDetector) {
        module->getPMT(pmt).addT0(fit.getT0(pmt));
      }
    }
  }


  if (writeFits) {
    out << hc << p1 << p2 << p3 << p4 << bg << cc;
  }

  out.Close();


  if (overwriteDetector) {

    NOTICE("Store calibration data on file " << detectorFile << endl);

    detector.comment.add(JMeta(argc, argv));

    store(detectorFile, detector);
  }

  if (pmtFile != "") {

    //----------------------------------------------------------
    // at this stage, fitted QE corresponds to hit probability
    //----------------------------------------------------------

    try {
      parameters.convertHitProbabilityToQE(mu);
    }
    catch(const JException& error) {
      FATAL(error.what());
    }

    ofstream out(pmtFile.c_str());

    out << parameters << endl;

    out.close();
  }

  return 0;
}