JEditPMTParameters.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>

#include "JLang/JLangToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTReadoutAddress.hh"
#include "JDetector/JPMTPhysicalAddress.hh"
#include "JDetector/JModuleAddressMap.hh"
#include "JDetector/JDetectorAddressMap.hh"
#include "JDetector/JDetectorSupportkit.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JSupport/JMeta.hh"
#include "JTools/JRange.hh"

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


namespace {

  using namespace JPP;

  static const char WILDCARD_RING       = '*';   //<! Wild card for ring of PMT in optical module.
  static const int  WILDCARD_POSITION   = -1;    //<! Wild card for position of PMT in ring.
  static const int  WILDCARD            = -1;    //<! Wild card for PMT identifier.


  /**
   * Compare PMT physical addresses taking into account wild cards.
   *
   * \param  first         first  PMT physical address
   * \param  second        second PMT physical address
   * \return               true if addresses are equal; else false
   */
  inline bool compare(const JPMTPhysicalAddress& first, const JPMTPhysicalAddress& second)
  {
    return (((first          == second))                                                                                          ||

            ((first.ring     == WILDCARD_RING     || second.ring     == WILDCARD_RING)     && first.position == second.position)  ||
            ((first.position == WILDCARD_POSITION || second.position == WILDCARD_POSITION) && first.ring     == second.ring)      ||

            ((first.ring     == WILDCARD_RING     || second.ring     == WILDCARD_RING)  && 
             (first.position == WILDCARD_POSITION || second.position == WILDCARD_POSITION)));
  }


  /**
   * Compare PMT identifiers taking into account wild cards.
   *
   * \param  first         first  PMT identifier
   * \param  second        second PMT identifier
   * \return               true if addresses are equal; else false
   */
  inline bool compare(const JPMTIdentifier& first, const JPMTIdentifier& second)
  {
    return (((first          == second))                                                                        ||

            ((first.getID()  == WILDCARD || second.getID()  == WILDCARD) && first.getTDC() == second.getTDC())  ||
            ((first.getTDC() == WILDCARD || second.getTDC() == WILDCARD) && first.getID()  == second.getID())   ||

            ((first.getID()  == WILDCARD || second.getID()  == WILDCARD)  && 
             (first.getTDC() == WILDCARD || second.getTDC() == WILDCARD)));
  }


  /**
   * Empty address.
   */
  struct JEmptyAddress {
    friend inline std::istream& operator>>(std::istream& in, JEmptyAddress& object) { return in; }
    friend inline std::ostream& operator<<(std::ostream& out, const JEmptyAddress& object) { return out; }
  };
  

  /**
   * Auxiliary class to apply modifications to PMT parameters.
   */
  template<class JAddress_t = JEmptyAddress>
  class JModifier {
  public:
    /**
     * Default constructor.
     */
    JModifier()
    {}


    /**
     * Apply modification to given parameters.
     *
     * \param  parameters  PMT parameters
     * \return             true if valid action; else false
     */
    bool apply(JPMTParameters& parameters) const
    {
      using namespace std;
      
      try {

        if        (this->action == "set") {

          parameters.getProperties().getValue<double>(this->key)  = this->value;

        } else if (this->action == "add") {

          parameters.getProperties().getValue<double>(this->key) += this->value;

        } else if (this->action == "sub") {

          parameters.getProperties().getValue<double>(this->key) -= this->value;

        } else if (this->action == "mul") {

          parameters.getProperties().getValue<double>(this->key) *= this->value;

        } else if (this->action == "div") {

          parameters.getProperties().getValue<double>(this->key) /= this->value;

        } else {
        
          return false;
        }
      }
      catch(const std::exception& error) {
        cerr << error.what() << endl;
        return false;
      }
      
      return true;
    }      


    /**
     * Read modifier from input.
     *
     * \param  in          input stream
     * \param  modifier    modifier
     * \return             input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JModifier& modifier)
    {
      return in >> modifier.address
                >> modifier.action
                >> modifier.key
                >> modifier.value;
    }


    /**
     * Write modifier to output.
     *
     * \param  out         output stream
     * \param  modifier    modifier
     * \return             output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JModifier& modifier)
    {
      return out << modifier.address << ' '
                 << modifier.action  << ' '
                 << modifier.key     << ' '
                 << modifier.value;
    }
  

    JAddress_t  address;
    std::string action;
    std::string key;
    double      value;
  };
}


/**
 * \file
 *
 * Auxiliary program to edit PMT parameters map.
 *
 * Syntax:
 * <pre>
 *     -@ "(set|add|sub|mul|div) <key> <value>"
 *     -A "<PMT physical address> (set|add|sub|mul|div) <key> <value>"
 *     -M "<PMT identifier>       (set|add|sub|mul|div) <key> <value>"
 * </pre>
 * In this, the PMT physical address corresponds to the data structure JDETECTOR::JPMTPhysicalAddress and
 * the PMT identifier to JDETECTOR::JPMTIdentifier.\n
 * The key corresponds to one of the data members of the JDETECTOR::JPMTParameters data structure.\n
 * The option <tt>-\@</tt> corresponds to the default PMT values.
 *
 * Note that in the absence of option <tt>-a</tt>, the detector identifier should be specified
 * using option <tt>-D</tt> so to obtain the correct PMT address mapping.
 *
 * The default value for option <tt>-E</tt> (expectation value for npe given two-fold coincidence) 
 * is taken from Analysis e-log entry <a href="https://elog.km3net.de/Analysis/519">519</a>.\n
 * The (default) values for option <tt>-T</tt> (time-over-threshold range) should correspond to 
 * the minimal and maximal value at option <tt>-t</tt> of application JCalibrateK40.cc.\n
 * The values for option <tt>-Q</tt> (QE range) can be used to take into account 
 * the light scaling factor applied in the simulation of the detector response.
 *
 * Multiple options <tt>-\@</tt>, <tt>-A</tt> and <tt>-M</tt> will be processed in order of appearance.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JRange<double>                 JRange_t;

  string               detectorFile;
  int                  detectorID;
  JPMTParametersMap    parameters;
  vector< JModifier<> >                     hdr;
  vector< JModifier<JPMTPhysicalAddress> >  mod;
  vector< JModifier<JPMTIdentifier> >       daq;
  JProxy<double>       mu  ("?", 2.05644e-01);
  JProxy<JRange_t>     T_ns("?", JRange_t(4, 250));
  JRange_t             QE;
  string               outputFile;
  bool                 squash;
  int                  debug;

  try {

    JParser<> zap("Auxiliary program to edit PMT parameters map.");

    zap['a'] = make_field(detectorFile,     "detector file.")                                                   = "";
    zap['D'] = make_field(detectorID,       "detector identifier (in absence of detector file).")               = 0;
    zap['P'] = make_field(parameters,       "PMT simulation data (or corresponding file name)")                 = JPARSER::initialised();
    zap['@'] = make_field(hdr,              "PMT parameter modifier for default values.")                       = JPARSER::initialised();
    zap['A'] = make_field(mod,              "PMT parameter modifier by physical address (e.g. B1).")            = JPARSER::initialised();
    zap['M'] = make_field(daq,              "PMT parameter modifier by DAQ address (e.g. <module> <channel>.")  = JPARSER::initialised();
    zap['E'] = make_field(mu,               "expectation value for npe given two-fold coincidence (" << mu  .getOption() << " -> " << mu  .getCustom() << ")")  = 0.0;
    zap['T'] = make_field(T_ns,             "time-over-threshold rang ("                             << T_ns.getOption() << " -> " << T_ns.getCustom() << ")")  = JRange_t();
    zap['Q'] = make_field(QE,               "QE range.")                                                        = JRange_t();
    zap['o'] = make_field(outputFile,       "output file.");
    zap['q'] = make_field(squash,           "squash meta data");
    zap['d'] = make_field(debug,            "debug level")                                                      = 2;

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


  if (squash) {
    parameters.comment.clear();
  }

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


  for (vector< JModifier<> >::const_iterator i = hdr.begin(); i != hdr.end(); ++i) {

    DEBUG("Modifying default PMT parameters " << i->action << ' ' << i->key << ' ' << i->value << endl);

    if (!i->apply(parameters.getDefaultPMTParameters())) {
      ERROR("No valid action: " << *i << endl);
    }
  }

  if (detectorFile != "") {

    // Setting default PMT parameters for given detector.
    
    JDetector detector;

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

    if        (detectorID == 0) {

      detectorID = detector.getID();
      
    } else if (detectorID != detector.getID()) {

      FATAL("Inconsistent detector identifier " << detectorID << " != " << detector.getID() << endl);
    }
      
    
    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {

      for (unsigned int pmt = 0; pmt != module->size(); ++pmt) {

        const JPMTIdentifier id(module->getID(), pmt);

        if (parameters.find(id) == parameters.end()) {

          DEBUG("Setting default parameters for PMT " << id << endl);
          
          parameters[id] = parameters.getDefaultPMTParameters();
        }
      }
    }
  }

  if (!mod.empty()) {

    if (!hasDetectorAddressMap(detectorID)) {
      FATAL("Invalid detector identifier " << detectorID << endl);
    }
    
    const JDetectorAddressMap& demo = getDetectorAddressMap(detectorID);
    
    for (JPMTParametersMap::iterator ps = parameters.begin(); ps != parameters.end(); ++ps) {

      const JPMTPhysicalAddress& address = demo.get(ps->first);

      for (vector< JModifier<JPMTPhysicalAddress> >::const_iterator i = mod.begin(); i != mod.end(); ++i) {
        
        if (compare(i->address, address)) {

          DEBUG("Modifying parameters for PMT " << ps->first << ' ' << i->action << ' ' << i->key << ' ' << i->value << endl);

          if (!i->apply(ps->second)) {
            ERROR("No valid action: " << *i << endl);
          }
        }
      }
    }
  }
  
  if (!daq.empty()) {

    for (JPMTParametersMap::iterator ps = parameters.begin(); ps != parameters.end(); ++ps) {

      for (vector< JModifier<JPMTIdentifier> >::const_iterator i = daq.begin(); i != daq.end(); ++i) {
        
        if (compare(ps->first, i->address)) {

          DEBUG("Modifying parameters for PMT " << ps->first << ' ' << i->action << ' ' << i->key << ' ' << i->value << endl);

          if (!i->apply(ps->second)) {
            ERROR("No valid action: " << *i << endl);
          }
        }
      }
    }
  }


  if        (mu > 0.0) {
    
    DEBUG("Correct measured QE for two-hit probability " << mu << endl);

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

    FATAL("Invalid expection value for two-hit probability " << mu << endl);
  }
  
  
  if (T_ns != JRange_t()) {

    DEBUG("Correct measured QE for time-over-threshold range " << T_ns << endl);

    const int NPE = 1;

    for (JPMTParametersMap::iterator i = parameters.begin(); i != parameters.end(); ++i) {

      const JPMTAnalogueSignalProcessor cpu(i->second);
      
      i->second.QE *= (cpu.getIntegralOfChargeProbability(i->second.threshold,
                                                          cpu.getNPE(T_ns.getUpperLimit()),
                                                          NPE)
                       /
                       cpu.getIntegralOfChargeProbability(cpu.getNPE(T_ns.getLowerLimit()),
                                                          cpu.getNPE(T_ns.getUpperLimit()),
                                                          NPE));
    }
  }
  
  
  if (QE != JRange_t()) {

    for (JPMTParametersMap::iterator i = parameters.begin(); i != parameters.end(); ++i) {
      i->second.QE = QE.constrain(i->second.QE);
    }
  }


  if (outputFile != "") {

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

    ofstream out(outputFile.c_str());

    out << parameters << endl;

    out.close();
  }
}