JEditPMTParametersMap.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/JDetectorAddressMapToolkit.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;

  /**
   * Wild card for physical address of PMT.
   */
  static const JPMTPhysicalAddress WILDCARD('*', -1);
  

  /**
   * Auxiliary class to apply modifications to PMT parameters.
   */
  template<class JAddress_t>
  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>
 *     -A "<PMT physical address> (set|add|sub|mul|div) <key> <value>"
 * </pre>
 * The PMT physical address corresponds to the data structure JDETECTOR::JPMTPhysicalAddress.
 * The key corresponds to one of the data members of the JDETECTOR::JPMTParameters data structure.
 *
 * 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.
 *
 * Multiple options <tt>-A</tt> will be processed in order of appearance.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JModifier<JPMTPhysicalAddress> JModifier_t;  
  typedef JRange<double>                 JRange_t;

  string               detectorFile;
  int                  detectorID;
  JPMTParametersMap    parameters;
  vector<JModifier_t>  modifier;
  JRange_t             T_ns;
  string               outputFile;
  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['A'] = make_field(modifier,         "PMT parameter modifier.")                             = JPARSER::initialised();
    zap['T'] = make_field(T_ns,             "time-over-threshold range.")                          = JRange_t();
    zap['o'] = make_field(outputFile,       "output file.")                                        = "pmt.txt";
    zap['d'] = make_field(debug)        = 3;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << 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()) {

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

  if (!modifier.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_t>::const_iterator i = modifier.begin(); i != modifier.end(); ++i) {
        
        if (i->address == address || i->address == WILDCARD) {

          NOTICE("Modifying parameters for PMT " << ps->first << endl);

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


  if (T_ns != JRange_t()) {

    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),
                                                          NPE)
                       /
                       cpu.getIntegralOfChargeProbability(cpu.getNPE(T_ns.getLowerLimit(), NPE),
                                                          cpu.getNPE(T_ns.getUpperLimit(), NPE),
                                                          NPE));
    }
  }
  
  
  if (outputFile != "") {

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

    ofstream out(outputFile.c_str());

    out << parameters << endl;

    out.close();
  }
}