JEditDetector.cc

Go to the documentation of this file.

#include <iostream>
#include <sstream>
#include <string>
#include <vector>

#include "TRandom3.h"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorHeader.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTIdentifier.hh"
#include "JDetector/JPMTPhysicalAddress.hh"
#include "JDetector/JModuleAddressMap.hh"
#include "JDetector/JDetectorAddressMap.hh"
#include "JDetector/JDetectorAddressMapToolkit.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JMath/JConstants.hh"
#include "JMath/JMath.hh"
#include "JTools/JRange.hh"
#include "JLang/JException.hh"
#include "JLang/JToken.hh"
#include "JSupport/JMeta.hh"

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


namespace {

  using namespace JPP;

  /**
   * Wild card for string identifier, module identifier or PMT address.
   */
  static const int         WILDCARD    =  -1;

  static const char        EOL         = ';';                 //!< end-of-line

  static const std::string reset_t     = "reset";             //!< Reset    time offset, position and PMT status bit
  static const std::string set_t       = "set";               //!< Set      time offset, position or  PMT status bit
  static const std::string add_t       = "add";               //!< Add      time offset or position
  static const std::string sub_t       = "sub";               //!< Subtract time offset or position
  static const std::string rot_t       = "rot";               //!< Rotate around z-axis by value [rad]
  static const std::string mul_t       = "mul";               //!< Multiply z-position by (1 + value)
  static const std::string via_t       = "via";               //!< Apply time offset by (floor * value)
  static const std::string assign_t    = "assign";            //!< Assign identifier

  static const std::string rand_t      = "rand";              //!< Random value(s)
  static const std::string randset_t   =  rand_t + set_t;     //!< Set      time offset or position
  static const std::string randadd_t   =  rand_t + add_t;     //!< Add      time offset or position
  static const std::string randsub_t   =  rand_t + sub_t;     //!< Subtract time offset or position
  static const std::string randrot_t   =  rand_t + rot_t;     //!< Rotate around z-axis by value [rad]
  static const std::string randmul_t   =  rand_t + mul_t;     //!< Multiply z-position by (1 + value)
  static const std::string randvia_t   =  rand_t + via_t;     //!< Apply time offset by (floor * value)

  static const std::string RESET_t     = "RESET";             //!< Reset    time offset and quaternion of module
  static const std::string SET_t       = "SET";               //!< Set      time offset or quaternion of module
  static const std::string ADD_t       = "ADD";               //!< Add      time offset or quaternion of module
  static const std::string SUB_t       = "SUB";               //!< Subtract time offset or quaternion of module


  /**
   * Auxiliary class to apply detector modifications.
   *
   * Note that the internal identifier may apply to a module as well as a string.
   */
  class JModifier {
  public:
    /**
     * Default constructor.
     */
    JModifier()
    {}


    /**
     * Check validity.
     *
     * \return                 true if valid modifier; else false
     */
    bool is_valid() const
    {
      return (action != "" && !data.empty());
    }


    /**
     * Apply modification to given module.
     *
     * \param  module      module
     * \return             true if valid action; else false
     */
    bool apply(JModule& module) const
    {
      switch (data.size()) {

      case 0:
        return  apply(module, action);                                                     // reset calibration
          
      case 1:
        return  apply(module, action, data[0]);                                            // time/position/orientation calibration
          
      case 2:
        return  apply(module, action, JVector3D(data[0], data[1], 0.0));                   // 2D position calibration
          
      case 3:
        return  apply(module, action, JVector3D(data[0], data[1], data[2]));               // 3D position calibration

      case 4:
        return  apply(module, action, JQuaternion3D(data[0], data[1], data[2], data[3]));  // 3D quaternion calibration

      default:
        return false;
      }
    }      


    /**
     * Read modifier from input.
     *
     * \param  in          input stream
     * \param  modifier    modifier
     * \return             input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JModifier& modifier)
    {
      using namespace std;

      if (in >> modifier.id >> modifier.action) {

        modifier.data.clear();

        for (double x; in >> x; ) {
          modifier.data.push_back(x);
        }

        in.clear(ios_base::eofbit);
      }

      return in;
    }


    /**
     * 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)
    {
      out << modifier.id;
      out << ' ';
      out << modifier.action;

      for (std::vector<double>::const_iterator i = modifier.data.begin(); i != modifier.data.end(); ++i) {
        out << ' ' << *i;
      }

      return out;
    }
  

    int                 id;
    std::string         action;
    std::vector<double> data;

  private:
    
    /**
     * Apply reset to given module.
     *
     * \param  module      module
     * \param  action      action
     * \return             true if valid action; else false
     */
    static bool apply(JModule& module, const std::string& action)
    {
      if        (action == reset_t) {

        module.set(JModule::getInstance().getPosition());

        for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
          pmt->setCalibration(JCalibration());
          pmt->setStatus(JStatus());
        }

        return true;

      } else if (action == RESET_t) {

        module.setCalibration(JCalibration());
        module.setQuaternion(JQuaternion3D());

        return true;

      } else if (action == SET_t) {

        if (module.getFloor() != 0)
          module.setT0(getAverage(make_array(module.begin(), module.end(), &JPMT::getT0))  -  PIEZO_DELAYTIME_US      * 1.0e+3);
        else
          module.setT0(0.0                                                                 -  HYDROPHONE_DELAYTIME_US * 1.0e+3);

        return true;

      } else {

        return false;
      }
    }

    
    /**
     * Apply time/position/orientation calibration to given module.
     *
     * \param  module      module
     * \param  action      action
     * \param  value       value
     * \return             true if valid action; else false
     */
    static bool apply(JModule& module, const std::string& action, const double value)
    {
      if        (action == set_t) {                          // actions with fixed values

        module.set(value);

      } else if (action == add_t) {

        module.add(value);

      } else if (action == sub_t) {

        module.sub(value);

      } else if (action == rot_t) {

        const JVector3D center = module.getPosition();

        module.sub(center);
        
        module.rotate(JRotation3Z(value));

        module.add(center);

      } else if (action == mul_t) {

        if (module.getFloor() != 0) {

          const JVector3D center(module.getPosition().getX(),
                                 module.getPosition().getY(),
                                 module.getPosition().getZ() * (1.0 + value));
        
          module.set(center);
        }

      } else if (action == via_t) {

        module.add(value * module.getFloor());

      } else if (action == SET_t) {

        module.setT0(value);

      } else if (action == ADD_t) {

        module.addT0(value);

      } else if (action == SUB_t) {

        module.subT0(value);

      } else if (action == randadd_t) {                      // actions with random values

        module.add(gRandom->Gaus(0.0, value));

      } else if (action == randsub_t) {
        
        module.sub(gRandom->Gaus(0.0, value));

      } else if (action == randrot_t){

        const JVector3D center = module.getPosition();

        module.sub(center);

        module.rotate(JRotation3Z(gRandom->Gaus(0.0, value)));

        module.add(center);
      
      } else if (action == randmul_t) {

        if (module.getFloor() != 0) {

          const JVector3D center(module.getPosition().getX(),
                                 module.getPosition().getY(),
                                 module.getPosition().getZ() * gRandom->Gaus(1.0, value));
        
          module.set(center);
        }

      } else if (action == randvia_t) {

        module.add(gRandom->Gaus(0.0, value) * module.getFloor());

      } else if (action == assign_t) {

        module.setID((int) value);

      } else {
        
        return false;
      }
      
      return true;
    }


    /**
     * Apply position calibration to given module.
     *
     * \param  module      module
     * \param  action      action
     * \param  pos         pos
     * \return             true if valid action; else false
     */
    static bool apply(JModule& module, const std::string& action, const JVector3D& pos)
    {
      const JVector3D randpos(gRandom->Gaus(0.0, pos.getX()),
                              gRandom->Gaus(0.0, pos.getY()),
                              gRandom->Gaus(0.0, pos.getZ()));

      if      (action == set_t)                          // actions with fixed values
        module.set(pos);
      else if (action == add_t)
        module.add(pos);
      else if (action == sub_t)
        module.sub(pos);
      else if (action == randset_t)                      // actions with random values
        module.set(randpos);
      else if (action == randadd_t)
        module.add(randpos);
      else if (action == randsub_t)
        module.sub(randpos);
      else
        return false;

      return true;
    }


    /**
     * Apply quaternion calibration to given module.
     *
     * \param  module      module
     * \param  action      action
     * \param  Q           quaternion
     * \return             true if valid action; else false
     */
    static bool apply(JModule& module, const std::string& action, const JQuaternion3D& Q)
    {
      if      (action == SET_t)
        module.setQuaternion(Q);
      else if (action == ADD_t)
        module.setQuaternion(Q * module.getQuaternion());
      else if (action == SUB_t)
        module.setQuaternion(Q.getConjugate() * module.getQuaternion());
      else
        return false;

      return true;
    }
  };


  /**
   * Get modifier for given string.
   *
   * \param  id          string identifier
   * \param  modifier    modifier
   * \return             modifier
   */
  inline const JModifier& getModifier(const int id, const JModifier& modifier)
  {
    using namespace std;

    static map<int, map<string, map<size_t, JModifier> > > buffer;

    const string::size_type pos = modifier.action.find(rand_t);

    if (pos != string::npos) {

      JModifier& result = buffer[id][modifier.action][modifier.data.size()];

      if (!result.is_valid()) {      

        result.id     = id;
        result.action = modifier.action.substr(pos + rand_t.length());

        for (size_t i = 0; i != modifier.data.size(); ++i) {
          result.data.push_back(gRandom->Gaus(0.0, modifier.data[i]));
        }
      }

      return result;

    } else {

      return modifier;
    }
  }


  /**
   * Auxiliary class to apply PMT status modifications.
   */
  template<class JAddress_t>
  class JPMTModifier :
    public JAddress_t
  {
  public:
    /**
     * Default constructor.
     */
    JPMTModifier()
    {}


    /**
     * Apply modification to given PMT.
     *
     * \param  pmt         PMT
     * \return             true if valid action; else false
     */
    bool apply(JPMT& pmt) const
    {
      try {

        if        (action == set_t) {

          pmt.set(getPMTStatusBit(value));
          
        } else if (action == reset_t) {
          
          pmt.reset(getPMTStatusBit(value));

        } else {
          
          return false;
        }

        return true;
      }
      catch(const std::exception&) {}

      return false;
    }
 

    /**
     * Read PMT modifier from input.
     *
     * \param  in          input stream
     * \param  modifier    modifier
     * \return             input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JPMTModifier& modifier)
    {
      return in >> static_cast<JAddress_t&>(modifier) >> modifier.action >> 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 JPMTModifier& modifier)
    {
      out << static_cast<const JAddress_t&>(modifier);
      out << ' ';
      out << modifier.action;
      out << ' ';
      out << modifier.value;

      return out;
    }
  

    std::string action;
    std::string value;
  };


  /**
   * Range of identifiers.
   */
  struct JRange_t :
    public JRange<int> {
    /**
     * Separator between two identifier values.
     */
    static const char SEPARATOR = '-';

    /**
     * Default constructor.
     */
    JRange_t() :
      JRange(-1, -1)
    {}


    /**
     * Read range from input.
     *
     * \param  in          input stream
     * \param  range       range
     * \return             input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JRange_t& range)
    {
      if (in >> range.first) {

        range.second = range.first;

        if (in.peek() == (int) JRange_t::SEPARATOR) {

          in.get();

          in >> range.second;

        } else {

          in.clear();
        }
      }

      return in;
    }


    /**
     * Write range to output.
     *
     * \param  out         output stream
     * \param  range       range
     * \return             output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JRange_t& range)
    {
      return out << range.first << JRange_t::SEPARATOR << range.second;
    }
  };
}


/**
 * \file
 *
 * Auxiliary program to modify detector calibration.
 *
 * Syntax:
 * <pre>
 *     -M "<module identifier>     (set|add|sub|randset|randadd|randsub) x0 [x1 [x2]]"
 *     -S "<string number>         (set|add|sub|randset|randadd|randsub) x0 [x1 [x2]]"
 *     -M "<module identifier>     (rot|randrot) phi"
 *     -S "<string number>         (rot|randrot) phi"
 *     -M "<module identifier>     (mul|randmul) factor"
 *     -S "<string number>         (mul|randmul) factor"
 *     -M "<module identifier>     (via|randvia) factor"
 *     -S "<string number>         (via|randvia) factor"
 *     -M "<module identifier>     (reset)"
 *     -S "<string number>         (reset)"
 *     -M "<module identifier>     (assign)  identifier"
 *     -M "<module identifier>     (SET|ADD|SUB|) x0 [x1 x2 x3]"
 *     -S "<string number>         (SET|ADD|SUB|) x0 [x1 x2 x3]"
 *     -M "<module identifier>     (SET)"
 *     -S "<string number>         (SET)"
 *     -P "<PMT identifier>        (set|reset)   (PMT_DISABLE|HIGH_RATE_VETO_DISABLE|FIFO_FULL_DISABLE|UDP_COUNTER_DISABLE|UDP_TRAILER_DISABLE|OUT_OF_SYNC)"
 *     -p "<PMT physical address>  (set|reset)   (PMT_DISABLE|HIGH_RATE_VETO_DISABLE|FIFO_FULL_DISABLE|UDP_COUNTER_DISABLE|UDP_TRAILER_DISABLE|OUT_OF_SYNC)"
 *     -k "<string number>[-<string number>"
 *     -r "<string number>[-<string number>"
 *     -m "<module identifier>"
 *     -D "<string number> <floor>"
 *     -\@ "<key>=<value>[;<key>=<value>"
 * </pre>
 * Options <tt>-M</tt> and <tt>-S</tt> refer to a module and a string, respectively.\n
 * The values provided for a string modification coherently apply to the modules of the specified string number.
 *
 * The options <tt>randxxx</tt> correspond to a randomisation of the specified option.
 *
 * If the module identifier or string number is -1,
 * the action is applied to all modules or strings in the detector, respectively.
 *
 * For options <tt>[rand]set</tt>, <tt>[rand]add</tt> and <tt>[rand]sub</tt>,
 * the number of values apply to position or time calibration in the following way:
 *     -#  time calibration     <tt>(t = x0)</tt>
 *     -#  position calibration <tt>(x = x0, y = x1, z = 0)</tt>
 *     -#  position calibration <tt>(x = x0, y = x1, z = x2)</tt>
 *
 * For options <tt>[rand]rot</tt>,
 * the angle <tt>phi</tt> refers to an anti-clockwise rotation around the z-axis.\n
 * The rotation angle is defined in radians.
 *
 * For options <tt>[rand]mul</tt>,
 * the multiplication <tt>factor</tt> (a.k.a.\ "stretching") applies to
 * the z-coordinates of the optical modules and not to the base module (read anchor).\n
 * The factor is defined as a fraction; the actual multiplication factor is <tt>(1 + factor)</tt>. 
 *
 * For options <tt>[rand]via</tt>,
 * the <tt>factor</tt> is multiplied with the floor number and then added to the time calibration.
 *
 * For options <tt>SET ADD SUB</tt>,
 * the number of values apply to time or quaternion calibration of the module in the following way:
 *     -#  time calibration       <tt>(t = x0)</tt>
 *     -#  invalid
 *     -#  invalid
 *     -#  quaternion calibration <tt>(qa = x0, qb = x1, qc = x2, qd = x3)</tt>
 *
 * If no values are given at the option <tt>SET</tt>,
 * the time calibration of the module is set to the average time calibration of the PMTs in that module, 
 * corrected for the delay time of the pieze sensor (JDETECTOR::PIEZO_DELAYTIME_US).\n
 * For the base module, the time calibration is set to 0, 
 * corrected for the delay time of the hydrophone (JDETECTOR::HYDROPHONE_DELAYTIME_US).\n
 * The unit of time calibration is <tt>ns</tt>.
 *
 * Note that for string modifiers with option <tt>randxxx</tt>,
 * the action is coherently applied to the modules in the specified string.\n
 * Only one type of action (defined by <tt>xxx</tt> and the number of values) is then allowed per string. 
 *
 * Option <tt>-\@</tt> refers to the header information.\n
 * The list of possible keys can be obtained using JPrintDetector.cc with option <tt>-O header</tt>.
 *
 * Multiple options <tt>-M</tt>, <tt>-S</tt> or <tt>-\@</tt> will be processed in order of appearance.
 *
 * Options <tt>-k</tt> and <tt>-r</tt> can be used to keep and remove (a range of) string numbers, respectively.
 *
 * The options <tt>-m</tt> and  <tt>-D</tt> can be used to maintain a specific module (and remove all others) and 
 * to delete a floor from a string, respectively.
 *
 * Note finally that if the output file name is the same as the input file name,
 * the original file will be overwritten.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JToken<';'>  JToken_t;

  string               inputFile;
  string               outputFile;
  vector<JModifier>    mod;
  vector<JModifier>    str;
  vector<JToken_t>     hdr;
  vector< JPMTModifier<JPMTIdentifier> >      pmt;
  vector< JPMTModifier<JPMTPhysicalAddress> > alt;
  vector<JRange_t>     keep;
  vector<JRange_t>     rm;
  vector<int>          id;
  multimap<int, int>   del;
  int                  debug;
  
  try { 

    JParser<> zap("Auxiliary program to modify detector.");

    zap['a'] = make_field(inputFile);
    zap['M'] = make_field(mod)           = JPARSER::initialised();
    zap['S'] = make_field(str)           = JPARSER::initialised();
    zap['@'] = make_field(hdr)           = JPARSER::initialised();
    zap['P'] = make_field(pmt)           = JPARSER::initialised();
    zap['p'] = make_field(alt)           = JPARSER::initialised();
    zap['o'] = make_field(outputFile);
    zap['k'] = make_field(keep)          = JPARSER::initialised();
    zap['r'] = make_field(rm)            = JPARSER::initialised();
    zap['m'] = make_field(id)            = JPARSER::initialised();
    zap['D'] = make_field(del)           = JPARSER::initialised();
    zap['d'] = make_field(debug)         = 2;

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

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


  gRandom->SetSeed(0);


  if ((keep.empty() ? 0 : 1  +
       rm  .empty() ? 0 : 1  +
       id  .empty() ? 0 : 1) > 1) {
    FATAL("Use either option -k, -r or -m." << endl);
  }


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

  if (detector.setToLatestVersion()) {
    NOTICE("Set detector version to " << detector.getVersion() << endl);
  }

  if (!hdr.empty()) {

    JProperties helper = detector.getProperties();
    
    for (vector<JToken_t>::const_iterator i = hdr.begin(); i != hdr.end(); ++i) {
      
      istringstream is(*i);
      
      is >> helper;
    }
  }


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

    bool __rm__ = !keep.empty() && rm.empty();

    for (vector<JRange_t>::const_iterator i = keep.begin(); i != keep.end(); ++i) {
      if (module->getString() >= i->first && module->getString() <= i->second) {
        __rm__ = false;
      }
    }

    for (vector<JRange_t>::const_iterator i = rm.begin(); i != rm.end(); ++i) {
      if (module->getString() >= i->first && module->getString() <= i->second) {
        __rm__ = true;
      }
    }

    if (!id.empty()) {
      __rm__ = find(id.begin(), id.end(), module->getID()) == id.end();
    }

    const auto range = del.equal_range(module->getString());

    for (auto i = range.first; i != range.second; ++i) {
      if (i->second == module->getFloor()) {
        __rm__ = true;
      }
    }

    if (__rm__)
      module = detector.erase(module);
    else
      ++module;
  }


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

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

      if (module->getID() == i->id || i->id == WILDCARD ){

        DEBUG("Modifier"                                                                                                << ' '
              << "(" << FILL(4,'0') << module->getString() << "," << FILL(2,'0') << module->getFloor() << FILL() << ")" << ' '
              << setw(10) << module->getID()                                                                            << ' '
              << "action" << ' ' << i->action << JEEPZ() << i->data << endl);
        
        if (!i->apply(*module)) {
          ERROR("No valid action: " << *i << endl);
        }
      }
    }
  }


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

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

      if (module->getString() == i->id || i->id == WILDCARD) {

        const JModifier& modifier = getModifier(module->getString(), *i);

        DEBUG("Modifier"                                                                                                << ' '
              << "(" << FILL(4,'0') << module->getString() << "," << FILL(2,'0') << module->getFloor() << FILL() << ")" << ' '
              << setw(10) << module->getID()                                                                            << ' '
              << "action" << ' ' << modifier.action << JEEPZ() << modifier.data << endl);
        
        if (!modifier.apply(*module)) {
          ERROR("No valid action: " << *i << endl);
        }
      }
    }
  }


  for (vector< JPMTModifier<JPMTIdentifier> >::const_iterator i = pmt.begin(); i != pmt.end(); ++i) {
  
    for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {

      if (module->getID() == i->getModuleID() || i->getModuleID() == WILDCARD) {

        DEBUG("PMT modifier"                                                                       << ' '
              << "(" << setw(10) << module->getID() << "," << setw(2) << i->getPMTAddress() << ")" << ' '
              << "action" << ' ' << i->action << ' ' 
              << "value"  << ' ' << i->value  << endl);

        if        (i->getPMTAddress() == WILDCARD) {

          for (int pmt = 0; pmt != getNumberOfPMTs(*module); ++pmt) {
            if (!i->apply(module->getPMT(pmt))) {
              ERROR("No valid action: " << *i << endl);
            }
          }

        } else if (i->getPMTAddress() <  0                         ||
                   i->getPMTAddress() >= getNumberOfPMTs(*module)  ||
                   !i->apply(module->getPMT(i->getPMTAddress()))) {
          ERROR("No valid action: " << *i << endl);
        }
      }
    }
  }

  if (!alt.empty()) {

    if (!hasDetectorAddressMap(detector.getID())) {
      FATAL("Invalid detector identifier " << detector.getID() << endl);
    }

    const JDetectorAddressMap& demo = getDetectorAddressMap(detector.getID());

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

      const JModuleAddressMap memo = demo.get(module->getID());

      for (vector< JPMTModifier<JPMTPhysicalAddress> >::const_iterator i = alt.begin(); i != alt.end(); ++i) {

        const int tdc = memo.getAddressTranslator(*i).tdc;

        DEBUG("PMT modifier"                                                        << ' '
              << "(" << setw(10) << module->getID() << "," << setw(2) << tdc << ")" << ' '
              << "action" << ' ' << i->action << ' ' 
              << "value"  << ' ' << i->value  << endl);

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


  try {
    store(outputFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
}