JDetectorToolkit.hh

Go to the documentation of this file.

#ifndef __JDETECTOR__JDETECTORTOOLKIT__
#define __JDETECTOR__JDETECTORTOOLKIT__

#include <limits>
#include <istream>
#include <ostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <set>
#include <algorithm>

#include "JDetector/JDetector.hh"
#include "JDetector/JMonteCarloDetector.hh"
#include "JDetector/JTimeRange.hh"
#include "JDetector/JModuleAddressMap.hh"
#include "JDetector/JLocation.hh"
#include "JDetector/JModuleGeometry.hh"
#include "JDetector/JStringCounter.hh"
#include "JDetector/JPMTIdentifier.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JTools/JConstants.hh"
#include "JIO/JFileStreamIO.hh"
#include "JLang/JString.hh"
#include "JLang/JException.hh"
#include "JLang/gzstream.h"
#include "JFit/JPoint3DEstimator.hh"
#include "JFit/JPoint3D.hh"
#include "Jeep/JPrint.hh"


/**
 * \author mdejong
 */

namespace JDETECTOR {}
namespace JPP { using namespace JDETECTOR; }

/**
 * Auxiliary classes and methods for detector calibration and simulation.
 */
namespace JDETECTOR {

  using JGEOMETRY3D::JVector3D;
  using JGEOMETRY3D::JVersor3D;
  using JTOOLS::getInverseSpeedOfLight;
  using JTOOLS::getIndexOfRefraction;
  using JTOOLS::getTanThetaC;
  using JTOOLS::getSinThetaC;
  using JTOOLS::PI;
  using JLANG::JString;
  using JLANG::JException;
  using JLANG::JFileOpenException;
  using JLANG::JProtocolException;
  using JGEOMETRY3D::JCylinder3D;
  using JFIT::JPoint3D;



  /**
   * File name extensions.
   */
  static const char* const GENDET_DETECTOR_FILE_FORMAT = ".det";      //!< file format used by gendet
  static const char* const BINARY_DETECTOR_FILE_FORMAT = ".dat";      //!< JIO binary file format
  static const char* const KM3NET_DETECTOR_FILE_FORMAT = ".detx";     //!< KM3NeT standard ASCII format
  static const char* const ZIPPED_DETECTOR_FILE_FORMAT = ".gz";       //!< zipped KM3NeT standard ASCII format
  static const char* const GDML_DETECTOR_FILE_FORMAT   = ".gdml";     //!< KM3Sim input format
  static const char* const GDML_SCHEMA                 = getenv("GDML_SCHEMA_DIR");  //!< directory necessary for correct GDML header output
  static const char* const CAN_MARGIN                  = getenv("CAN_MARGIN");       //!< extension of the detector size to comply with the can definition
  static const char* const G4GDML_DEFAULT_SCHEMALOCATION = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd";


  /**
   * Get maximal distance between modules in detector.
   *
   * \param  detector          detector
   * \return                   maximal distance [m]
   */
  inline double getMaximalDistance(const JDetector& detector)
  {
    double dmax = 0.0;
    
    for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
      for (JDetector::const_iterator j = detector.begin(); j != i; ++j) {
        if (i->getDistance(*j) > dmax) {
          dmax = i->getDistance(*j);
        }
      }
    }
    
    return dmax;
  }


  /**
   * Get rotation over X axis in Geant4 coordinate system
   *
   * \param  dir            direction
   * \return                X-rotation [deg]
   */
  inline double GetXrotationG4(const JVersor3D dir)
  {
    const double phi = atan2(dir.getDY(), dir.getDZ())*(180.0/PI);

    if (phi < 0.0){
      return phi + 360.0;
    }
    else{
      return phi;       
    }
  }


  /**
   * Get rotation over Y axis in Geant4 coordinate system
   *
   * \param  dir            direction
   * \return                Y-rotation [deg]
   */
  inline double GetYrotationG4(const JVersor3D dir)
  {
    return asin(-dir.getDX())*(180.0/PI);
  }


  inline void read_gdml(std::istream&, const JDetector&)
  {}


  /**
   * Writes KM3Sim GDML input file from detector
   *
   * \param  out           output stream
   * \param  detector      detector
   */
  inline void write_gdml(std::ostream& out, const JDetector& detector)
  {
      

    const JCylinder3D cylinder(detector.begin(), detector.end());
    double can_margin;
    if(CAN_MARGIN) can_margin = atof(CAN_MARGIN);
    else{
         std::cout<<"CAN_MARGIN not defined! Setting standard CAN_MARGIN = 350 m"<<std::endl;
         can_margin = 5*70; //this is given in meters like all the dimensions in the GDML file (look at the unit field of every position and solid definition)
    }
    // Detector height is bascially just a distance between the highest and the lowest DOM
    //double Detector_Height = cylinder.getZmax() - cylinder.getZmin();
    // Distance from the center of the lowest DOM to the seabed
    //double Seabed_Distance = cylinder.getZmin();
                
    // To translate Z coordinate from detx format which corresponds to the distance from the seabed one has to perform a following calculation
    // Dom_Z_GDML = Dom_Z_Detx - Distance to the seabed from the lowest DOM - Detector_height/2

    //double DistanceOfSeabedFromLowestOM = Seabed_Distance;
    // For now, the size of the WorldBox is fixed to 2200 meters
    //const double WorldBoxHeight = 2200.;

    const double WorldCylinderHeight = 2*(cylinder.getZmax() - cylinder.getZmin() + can_margin + 50); //
    const double WorldRadius = cylinder.getRadius() + can_margin + 50;

    std::cout<<WorldCylinderHeight/2<<std::endl;
    std::cout<<WorldRadius<<std::endl;

    // Below formula comes directly from the note by Apostolos, can be used to reproduce the first GDML format of km3sim
    //double Crust_Z_size = WorldBoxHeight/2-Detector_Height/2-DistanceOfSeabedFromLowestOM; Formula for obsolete coordinate covention, where Z_seabed != 0
    //double Crust_Z_position = -WorldBoxHeight/4-Detector_Height/4-DistanceOfSeabedFromLowestOM/2; Formula for obsolete coordinate covention, where Z_seabed != 0

    //use this if you want to use WorldBox as the World Volume
    //const double Crust_Z_size = WorldBoxHeight/2; 
    //const double Crust_Z_position= -WorldBoxHeight/4;

    const double Crust_Z_size = WorldCylinderHeight/2;
    const double Crust_Z_position= -WorldCylinderHeight/4;                      

    out<<"<?xml version=\"1.0\" encoding=\"UTF-8\" ?>\n<gdml xmlns:gdml=\"http://cern.ch/2001/Schemas/GDML\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xsi:noNamespaceSchemaLocation=\"";
    if(!GDML_SCHEMA){
      std::cout<<"GDML_SCHEMA_DIR NOT DEFINED! Setting default path"<<std::endl;
      //out<<"http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd\">\n\n\n";
      out<<G4GDML_DEFAULT_SCHEMALOCATION<<"\">\n\n\n";
    }
    else{
      out<<GDML_SCHEMA<<"gdml.xsd\">\n\n\n";
    }
    out<<"<define>\n";
    out<<"<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";

    int PMTs_NO = 0;
    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
                        
      const JFIT::JEstimator<JPoint3D> center(module->begin(), module->end());

      //out<<"<position name=\"PosString"<<module->getString()<<"_Dom"<<module->getID()<<"\" unit=\"m\" x=\""<<module->getX()<<"\" y=\""<<module->getY()<<"\" z=\""<<module->getZ() - Seabed_Distance-Detector_Height/2<<"\"/>\n"; Formula for obsolete coordinate covention, where Z_seabed != 0
      out<<"<position name=\"PosString"<<module->getString()<<"_Dom"<<module->getID()<<"\" unit=\"m\" x=\""<<module->getX()<<"\" y=\""<<module->getY()<<"\" z=\""<<module->getZ()<<"\"/>\n";
                                        
      for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
        //JVector3D vec;
        //static_cast<JVector3D> (vec) = module->sub(*pmt);
                                
        const JVector3D vec = static_cast<JVector3D>(*pmt).sub(center);
        out<<"<position name=\"CathodPosition"<<pmt->getID()<<"_"<<module->getID()<<"\" unit=\"m\" x=\""<<vec.getX()<<"\" y=\""<<vec.getY()<<"\" z=\""<<vec.getZ()<<"\"/>\n";
        out<<"<rotation name=\"CathodRotation"<<pmt->getID()<<"_"<<module->getID()<<"\" unit=\"deg\" x=\""<<GetXrotationG4(*pmt)<<"\" y=\""<<GetYrotationG4(*pmt)<<"\" z=\"0.000000\"/>\n";
        out<<"<constant name=\"CathodID_"<<PMTs_NO<<"\" value=\""<<pmt->getID()<<"\"/>\n";
        PMTs_NO++;
      }
                        
    }

    out<<"<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
    out<<"\n\n\n";
    out<<"<!-- end of DU position definitions -->\n<position name=\"CrustPosition\"   unit=\"m\" x=\"0\" y=\"0\" z=\""<<Crust_Z_position<<"\"/>\n";
                
    out<<"</define>\n";
    out<<"<materials>\n";
    out<<"</materials>\n";
    out<<"<solids>\n";
    //out<<"    <box name=\"WorldBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"2200\"/>\n";
    out<<"      <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\""<<Crust_Z_size<<"\"/>\n";
    out<<"      <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
    out<<"      <sphere name=\"OMSphere\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"21.6\" startphi=\"0.0\" deltaphi=\"360.0\" starttheta=\"0.0\" deltatheta=\"180.0\"/>\n";
    out<<"      <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
    out<<"      <tube name=\"WorldTube\" lunit=\"m\" aunit=\"deg\" rmin=\"0.0\" rmax=\""<<WorldRadius<<"\" z=\""<<WorldCylinderHeight<<"\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";    
    out<<"</solids>\n\n\n";
                
    out<<"<structure>\n";
    out<<"      <volume name=\"CathodVolume\">\n";
    out<<"              <materialref ref=\"Cathod\"/>\n";
    out<<"              <solidref ref=\"CathodTube\"/>\n";
    out<<"      </volume>\n";

    out<<"<!-- OMVolume(s) construction -->\n";

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

      out<<"    <volume name=\"OMVolume"<<module->getID()<<"\">\n";
      out<<"            <materialref ref=\"Water\"/>\n";
      out<<"            <solidref ref=\"OMSphere\"/>\n";

      for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
        out<<"          <physvol>\n";
        out<<"                  <volumeref ref=\"CathodVolume\"/>\n";
        out<<"                  <positionref ref=\"CathodPosition"<<pmt->getID()<<"_"<<module->getID()<<"\"/>\n";
        out<<"                  <rotationref ref=\"CathodRotation"<<pmt->getID()<<"_"<<module->getID()<<"\"/>\n";
        out<<"          </physvol>\n";
      }

      out<<"    </volume>\n";
    }

    out<<"<!-- StoreyVolume(s) construction -->\n";

    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
      out<<"    <volume name=\"StoreyVolume"<<module->getID()<<"\">\n";
      out<<"            <materialref ref=\"Water\"/>\n";
      out<<"            <solidref ref=\"StoreyBox\"/>\n";
      out<<"            <physvol>\n";
      out<<"                    <volumeref ref=\"OMVolume"<<module->getID()<<"\"/>\n";
      out<<"                    <positionref ref=\"OMShift\"/>\n";
      //out<<"                  <positionref ref=\"zero\"/>\n";
      out<<"                    <rotationref ref=\"identity\"/>\n";
      out<<"            </physvol>\n";
      out<<"    </volume>\n";
    }

    out<<"<!-- Crust Volume construction-->\n";
    out<<"<volume name=\"CrustVolume\">\n";
    out<<"      <materialref ref=\"Crust\"/>\n";
    out<<"      <solidref ref=\"CrustBox\"/>\n";
    out<<"</volume>\n";

    out<<"<!-- World Volume construction-->\n";
    out<<"<volume name=\"WorldVolume\">\n";
    out<<"      <materialref ref=\"Water\"/>\n";
    //out<<"    <solidref ref=\"WorldBox\"/>\n";
    out<<"      <solidref ref=\"WorldTube\"/>\n";

    out<<"      <physvol>\n";
    out<<"              <volumeref ref=\"CrustVolume\"/>\n";
    out<<"              <positionref ref=\"CrustPosition\"/>\n";
    out<<"              <rotationref ref=\"identity\"/>\n";
    out<<"      </physvol>\n";

    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
      out<<"    <physvol>\n";
      out<<"            <volumeref ref=\"StoreyVolume"<<module->getID()<<"\"/>\n";
      out<<"            <positionref ref=\"PosString"<< module->getString() <<"_Dom"<< module->getID() <<"\"/>\n";
      //out<<"          <rotation x="{}" y="{}" z="0.0"/>\n'.format(Dom_Phi[i], Dom_Theta[i]))
      out<<"                    <rotationref ref=\"identity\"/>\n";
      out<<"    </physvol>\n";  
    }

    out<<"</volume>\n";

    out<<"</structure>\n";
    out<<"<setup name=\"Default\" version=\"1.0\">\n";
    out<<"<world ref=\"WorldVolume\"/>\n";
    out<<"</setup>\n";
    out<<"</gdml>\n";
  }


  /**
   * Get maximal time between modules in detector following causality.
   *
   * \param  detector          detector
   * \return                   maximal time [ns]
   */
  inline double getMaximalTime(const JDetector& detector)
  {
    return getMaximalDistance(detector) * getIndexOfRefraction() * getInverseSpeedOfLight();
  }


  /**
   * Get maximal time between modules in detector following causality.
   * The road width corresponds to the maximal distance traveled by the light.
   *
   * \param  detector          detector
   * \param  roadWidth_m       road width [m]
   * \return                   maximal time [ns]
   */
  inline double getMaximalTime(const JDetector& detector, const double roadWidth_m)
  {
    const double Dmax_m = getMaximalDistance(detector);

    return (sqrt((Dmax_m + roadWidth_m*getSinThetaC()) * 
                 (Dmax_m - roadWidth_m*getSinThetaC()))  +  
            roadWidth_m * getSinThetaC() * getTanThetaC()) * getInverseSpeedOfLight();
  }


  /**
   * Get de-calibrated time range.
   *
   * The de-calibrated time range is corrected for minimal and maximal time offset of PMTs in given module.
   *
   * \param  timeRange         time range [ns]
   * \param  module            module
   * \return                   time range [ns]
   */
  inline JTimeRange getTimeRange(const JTimeRange& timeRange, const JModule& module)
  {
    if (!module.empty()) {

      JTimeRange time_range(JTimeRange::DEFAULT_RANGE);
      
      for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
        
        const JCalibration& calibration = pmt->getCalibration();
        
        time_range.include(putTime(timeRange.getLowerLimit(), calibration));
        time_range.include(putTime(timeRange.getUpperLimit(), calibration));
      }
      
      return time_range;

    } else {
      
      return timeRange;
    }
  }

  
  /**
   * Get number of PMTs.
   *
   * \param  module            module
   * \return                   number of PMTs
   */
  inline int getNumberOfPMTs(const JModule& module)
  {
    return module.size();
  }

  
  /**
   * Get number of PMTs.
   *
   * \param  detector          detector
   * \return                   number of PMTs
   */
  inline int getNumberOfPMTs(const JDetector& detector)
  {
    int number_of_pmts = 0;
    
    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
      number_of_pmts += getNumberOfPMTs(*module);
    }

    return number_of_pmts;
  }

  
  /**
   * Get list of strings  IDs.
   *
   * \param  detector          detector
   * \return                   list of string IDs
   */
  inline std::set<int> getStringIDs(const JDetector& detector)
  {
    std::set<int> buffer;
    
    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
      buffer.insert(module->getString());
    }

    return buffer;
  }

  
  /**
   * Get number of floors.
   *
   * \param  detector          detector
   * \return                   number of floors
   */
  inline int getNumberOfFloors(const JDetector& detector)
  {
    std::set<int> buffer;
    
    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
      buffer.insert(module->getFloor());
    }

    return buffer.size();
  }

  
  /**
   * Get number of modules.
   *
   * \param  detector          detector
   * \return                   number of modules
   */
  inline int getNumberOfModules(const JDetector& detector)
  {
    std::set<int> buffer;
    
    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
      buffer.insert(module->getID());
    }

    return buffer.size();
  }
  
  
  /**
   * Load detector from input file.
   *
   * Supported file formats:
   *   - ASCII   file, extension JDETECTOR::GENDET_DETECTOR_FILE_FORMAT, gendet          format; 
   *   - binary  file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp    internal format;
   *   - ASCII   file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, KM3NeT standard format; 
   *   - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, KM3NeT standard format; 
   *
   * \param  file_name         file name
   * \param  detector          detector
   */
  inline void load(const JString& file_name, JDetector& detector)
  {
    using namespace std;
    using namespace JIO;
    
    if      (file_name.endsWith(GENDET_DETECTOR_FILE_FORMAT)) {

      JMonteCarloDetector buffer(true);

      ifstream in(file_name.c_str());

      if (!in) {
        THROW(JFileOpenException, "File not opened: " << file_name);
      }
      
      in >> buffer;

      in.close();

      detector.swap(buffer);

    } else if (file_name.endsWith(BINARY_DETECTOR_FILE_FORMAT)) {

      JFileStreamReader in(file_name.c_str());

      if (!in) {
        THROW(JFileOpenException, "File not opened: " << file_name);
      }
      
      detector.read(in);

      in.close();
 
    } else if (file_name.endsWith(KM3NET_DETECTOR_FILE_FORMAT)) {
      
      ifstream in(file_name.c_str());

      if (!in) {
        THROW(JFileOpenException, "File not opened: " << file_name);
      }
      
      in >> detector;

      in.close();

    } else if (file_name.endsWith(ZIPPED_DETECTOR_FILE_FORMAT)) {
      
      igzstream in(file_name.c_str());

      if (!in) {
        THROW(JFileOpenException, "File not opened: " << file_name);
      }
      
      in >> detector;

      in.close();

    } else {

      THROW(JProtocolException, "Protocol not defined: " << file_name);
    }
  }
  
  
  /**
   * Store detector to output file.
   *
   * Supported file formats:
   *   - binary  file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
   *   - ASCII   file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, KM3NeT standard format; 
   *   - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, KM3NeT standard format; 
   *   - gdml    file, extension JDETECTOR::GDML_DETECTOR_FILE_FORMAT,   KM3Sim input format;
   *
   * \param  file_name         file name
   * \param  detector          detector
   */
  inline void store(const JString& file_name, const JDetector& detector)
  {
    using namespace std;
    using namespace JIO;
    
    if        (file_name.endsWith(BINARY_DETECTOR_FILE_FORMAT)) {

      JFileStreamWriter out(file_name.c_str());

      detector.write(out);

      out.close();
 
    } else if (file_name.endsWith(KM3NET_DETECTOR_FILE_FORMAT)) {
      
      std::ofstream out(file_name.c_str());

      out << detector;

      out.close();

    } else if (file_name.endsWith(ZIPPED_DETECTOR_FILE_FORMAT)) {
      
      ogzstream out(file_name.c_str());

      out << detector;

      out.close();

    } else if (file_name.endsWith(GDML_DETECTOR_FILE_FORMAT)) {
      
      std::ofstream out(file_name.c_str());

      write_gdml(out, detector);

      out.close();

    }  else {

      THROW(JProtocolException, "Protocol not defined: " << file_name);
    }
  }


  /**
   * Get module according module address map.
   *
   * \param  memo              module address map
   * \param  id                module identifier
   * \param  location          module location
   * \return                   module
   */
  inline const JModule& getModule(const JModuleAddressMap& memo, 
                                  const int                id,
                                  const JLocation&         location = JLocation())
  {
    static JModule module;


    module.setID(id);

    module.setLocation(location);

    module.resize(memo.size());
    
    if (memo.has( 0)) { module[memo[ 0].tdc] = JPMT( 1, JAxis3D(JVector3D(+0.000, +0.000, -0.200), JVersor3D(+0.000, +0.000, -1.000))); }    

    if (memo.has( 1)) { module[memo[ 1].tdc] = JPMT( 2, JAxis3D(JVector3D(+0.000, +0.105, -0.170), JVersor3D(+0.000, +0.527, -0.850))); }
    if (memo.has( 2)) { module[memo[ 2].tdc] = JPMT( 3, JAxis3D(JVector3D(+0.091, +0.053, -0.170), JVersor3D(+0.456, +0.263, -0.850))); }
    if (memo.has( 3)) { module[memo[ 3].tdc] = JPMT( 4, JAxis3D(JVector3D(+0.091, -0.053, -0.170), JVersor3D(+0.456, -0.263, -0.850))); }
    if (memo.has( 4)) { module[memo[ 4].tdc] = JPMT( 5, JAxis3D(JVector3D(+0.000, -0.105, -0.170), JVersor3D(+0.000, -0.527, -0.850))); }
    if (memo.has( 5)) { module[memo[ 5].tdc] = JPMT( 6, JAxis3D(JVector3D(-0.091, -0.053, -0.170), JVersor3D(-0.456, -0.263, -0.850))); }
    if (memo.has( 6)) { module[memo[ 6].tdc] = JPMT( 7, JAxis3D(JVector3D(-0.091, +0.053, -0.170), JVersor3D(-0.456, +0.263, -0.850))); }

    if (memo.has( 7)) { module[memo[ 7].tdc] = JPMT( 8, JAxis3D(JVector3D(+0.083, +0.144, -0.111), JVersor3D(+0.416, +0.720, -0.555))); }
    if (memo.has( 8)) { module[memo[ 8].tdc] = JPMT( 9, JAxis3D(JVector3D(+0.166, +0.000, -0.111), JVersor3D(+0.832, +0.000, -0.555))); }
    if (memo.has( 9)) { module[memo[ 9].tdc] = JPMT(10, JAxis3D(JVector3D(+0.083, -0.144, -0.111), JVersor3D(+0.416, -0.720, -0.555))); }
    if (memo.has(10)) { module[memo[10].tdc] = JPMT(11, JAxis3D(JVector3D(-0.083, -0.144, -0.111), JVersor3D(-0.416, -0.720, -0.555))); }
    if (memo.has(11)) { module[memo[11].tdc] = JPMT(12, JAxis3D(JVector3D(-0.166, +0.000, -0.111), JVersor3D(-0.832, +0.000, -0.555))); }
    if (memo.has(12)) { module[memo[12].tdc] = JPMT(13, JAxis3D(JVector3D(-0.083, +0.144, -0.111), JVersor3D(-0.416, +0.720, -0.555))); }

    if (memo.has(13)) { module[memo[13].tdc] = JPMT(14, JAxis3D(JVector3D(+0.000, +0.191, -0.059), JVersor3D(+0.000, +0.955, -0.295))); }
    if (memo.has(14)) { module[memo[14].tdc] = JPMT(15, JAxis3D(JVector3D(+0.165, +0.096, -0.059), JVersor3D(+0.827, +0.478, -0.295))); }
    if (memo.has(15)) { module[memo[15].tdc] = JPMT(16, JAxis3D(JVector3D(+0.165, -0.096, -0.059), JVersor3D(+0.827, -0.478, -0.295))); }
    if (memo.has(16)) { module[memo[16].tdc] = JPMT(17, JAxis3D(JVector3D(+0.000, -0.191, -0.059), JVersor3D(+0.000, -0.955, -0.295))); }
    if (memo.has(17)) { module[memo[17].tdc] = JPMT(18, JAxis3D(JVector3D(-0.165, -0.096, -0.059), JVersor3D(-0.827, -0.478, -0.295))); }
    if (memo.has(18)) { module[memo[18].tdc] = JPMT(19, JAxis3D(JVector3D(-0.165, +0.096, -0.059), JVersor3D(-0.827, +0.478, -0.295))); }

    if (memo.has(19)) { module[memo[19].tdc] = JPMT(20, JAxis3D(JVector3D(+0.096, +0.165, +0.059), JVersor3D(+0.478, +0.827, +0.295))); }
    if (memo.has(20)) { module[memo[20].tdc] = JPMT(21, JAxis3D(JVector3D(+0.191, +0.000, +0.059), JVersor3D(+0.955, +0.000, +0.295))); }
    if (memo.has(21)) { module[memo[21].tdc] = JPMT(22, JAxis3D(JVector3D(+0.096, -0.165, +0.059), JVersor3D(+0.478, -0.827, +0.295))); }
    if (memo.has(22)) { module[memo[22].tdc] = JPMT(23, JAxis3D(JVector3D(-0.096, -0.165, +0.059), JVersor3D(-0.478, -0.827, +0.295))); }
    if (memo.has(23)) { module[memo[23].tdc] = JPMT(24, JAxis3D(JVector3D(-0.191, +0.000, +0.059), JVersor3D(-0.955, +0.000, +0.295))); }
    if (memo.has(24)) { module[memo[24].tdc] = JPMT(25, JAxis3D(JVector3D(-0.096, +0.165, +0.059), JVersor3D(-0.478, +0.827, +0.295))); }

    if (memo.has(25)) { module[memo[25].tdc] = JPMT(26, JAxis3D(JVector3D(+0.000, +0.166, +0.111), JVersor3D(+0.000, +0.832, +0.555))); }
    if (memo.has(26)) { module[memo[26].tdc] = JPMT(27, JAxis3D(JVector3D(+0.144, +0.083, +0.111), JVersor3D(+0.720, +0.416, +0.555))); }
    if (memo.has(27)) { module[memo[27].tdc] = JPMT(28, JAxis3D(JVector3D(+0.144, -0.083, +0.111), JVersor3D(+0.720, -0.416, +0.555))); }
    if (memo.has(28)) { module[memo[28].tdc] = JPMT(29, JAxis3D(JVector3D(+0.000, -0.166, +0.111), JVersor3D(+0.000, -0.832, +0.555))); }
    if (memo.has(29)) { module[memo[29].tdc] = JPMT(30, JAxis3D(JVector3D(-0.144, -0.083, +0.111), JVersor3D(-0.720, -0.416, +0.555))); }
    if (memo.has(30)) { module[memo[30].tdc] = JPMT(31, JAxis3D(JVector3D(-0.144, +0.083, +0.111), JVersor3D(-0.720, +0.416, +0.555))); }

    module.compile();

    return module;
  }


  /**
   * Get module corresponding to Antares storey.
   *
   * \param  id                module identifier
   * \param  location          module location
   * \return                   module
   */
  inline const JModule& getModule(const int        id,
                                  const JLocation& location = JLocation())
  {
    using JTOOLS::PI;

    static JModule module;

    module.setID(id);

    module.setLocation(location);

    if (module.empty()) {

      module.resize(3);
      
      const double R  = 0.5;           // [m]
      
      const double st = sin(0.75*PI);
      const double ct = cos(0.75*PI);
      
      for (int i = 0; i != 3; ++i) {
        
        const double phi = (2.0*PI*i) / 3.0;
        const double cp  = cos(phi);
        const double sp  = sin(phi);
        
        module[i] = JPMT(i, JAxis3D(JVector3D(R*st*cp, R*st*sp, R*ct), JVersor3D(st*cp, st*sp, ct)));
      }
    }

    return module;
  }


  /**
   * Auxiliary class to get rotation matrix between two optical modules.
   */
  struct JRotation :
    public JRotation3D
  {

    static const size_t NUMBER_OF_DIMENSIONS = 3;    //!< Number of dimensions
    

    /**
     * Get rotation matrix to go from first to second module.
     *
     * \param  first       first  module
     * \param  second      second module
     * \return             rotation matrix
     */ 
    const JRotation3D& operator()(const JModule& first, const JModule& second) 
    {
      this->setIdentity();

      if (first.size() == second.size()) {

        const size_t N = first.size();

        if (N >= NUMBER_OF_DIMENSIONS) {

          in .resize(N);
          out.resize(N);
        
          for (size_t i = 0; i != N; ++i) {
            in [i] = first .getPMT(i).getDirection();
            out[i] = second.getPMT(i).getDirection();
          }

          for (size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
            if (!orthonormalise(i)) {
              THROW(JException, "Failure to orthonormalise direction " << i);
            }
          }

          this->a00 = out[0].getX() * in[0].getX()  +  out[1].getX() * in[1].getX()  +  out[2].getX() * in[2].getX();
          this->a01 = out[0].getX() * in[0].getY()  +  out[1].getX() * in[1].getY()  +  out[2].getX() * in[2].getY();
          this->a02 = out[0].getX() * in[0].getZ()  +  out[1].getX() * in[1].getZ()  +  out[2].getX() * in[2].getZ();
          
          this->a10 = out[0].getY() * in[0].getX()  +  out[1].getY() * in[1].getX()  +  out[2].getY() * in[2].getX();
          this->a11 = out[0].getY() * in[0].getY()  +  out[1].getY() * in[1].getY()  +  out[2].getY() * in[2].getY();
          this->a12 = out[0].getY() * in[0].getZ()  +  out[1].getY() * in[1].getZ()  +  out[2].getY() * in[2].getZ();
          
          this->a20 = out[0].getZ() * in[0].getX()  +  out[1].getZ() * in[1].getX()  +  out[2].getZ() * in[2].getX();
          this->a21 = out[0].getZ() * in[0].getY()  +  out[1].getZ() * in[1].getY()  +  out[2].getZ() * in[2].getY();
          this->a22 = out[0].getZ() * in[0].getZ()  +  out[1].getZ() * in[1].getZ()  +  out[2].getZ() * in[2].getZ();

        } else {

          THROW(JException, "Module size" << N << " < " << NUMBER_OF_DIMENSIONS);
        }

      } else {

        THROW(JException, "Module size" << first.size() << " != " << second.size());
      }

      return *this;
    }

  private:
    /**
     * Put normalised primary direction at specified index and orthoganilise others.\n
     * This procedure follows Gram-Schmidt process.
     *
     * \param  index       index 
     * \return             true if primary direction exists; else false
     */
    bool orthonormalise(const size_t index)
    {
      using namespace std;

      size_t pos = index;

      for (size_t i = index + 1; i != in.size(); ++i) {
        if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
          pos = i;
        }
      }

      const double u = in[pos].getLength();

      if (u > 0.0) {

        in [pos] /= u;
        out[pos] /= u;

        if (pos != index) {
          swap(in [pos], in [index]);
          swap(out[pos], out[index]);
        }

        for (size_t i = index + 1; i != in.size(); ++i) {

          const double dot = in[index].getDot(in[i]);

          in [i] -= dot * in [index];
          out[i] -= dot * out[index];
        }

        return true;

      } else {

        return false;
      }
    }


    std::vector<JVector3D> in;
    std::vector<JVector3D> out;
  };

  
  /**
   * Function object to get rotation matrix to go from first to second module.
   */
  static JRotation getRotation;

  
  /**
   * Get position to go from first to second module.
   *
   * \param  first       first  module
   * \param  second      second module
   * \return             position
   */ 
  inline JPosition3D getPosition(const JModule& first, const JModule& second)
  {
    return second.getPosition() - first.getPosition();
  }

  
  /**
   * Get calibration to go from first to second calibration.
   *
   * \param  first       first  calibration
   * \param  second      second calibration
   * \return             calibration
   */ 
  inline JCalibration getCalibration(const JCalibration& first, const JCalibration& second)
  {
    return JCalibration(second.getT0() - first.getT0());
  }
}

#endif