#ifndef __JDETECTOR__JMONTECARLODETECTOR__
#define __JDETECTOR__JMONTECARLODETECTOR__

#include <string>
#include <istream>
#include <vector>
#include <cmath>
#include <algorithm>

#include "JMath/JConstants.hh"
#include "JLang/JObjectID.hh"
#include "JDetector/JDetector.hh"
#include "JGeometry3D/JEulerAngle3D.hh"
#include "JGeometry3D/JEulerMatrix3D.hh"
#include "JGeometry3D/JAngle3D.hh"


/**
 * \author mdejong
 */

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

namespace JDETECTOR {

  using JMATH::PI;
  using JLANG::JObjectID;
  using JLANG::getUndefinedObjectID;
  using JGEOMETRY3D::JEulerAngle3D;
  using JGEOMETRY3D::JEulerMatrix3D;
  using JGEOMETRY3D::JAngle3D;
  using JGEOMETRY3D::JDirection3D;
  using JGEOMETRY3D::JAxis3D;


  static const int LED_BEACON_PMT_TYPE = 2;          //!< PMT type of LED beacon


  /**
   * Monte Carlo detector (i.e.\ \".det\" file).
   */
  class JMonteCarloDetector : 
    public JDetector
  {
  private:
    /**
     * Auxiliary class for OM 
     */
    class OM : 
      public JObjectID 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      OM(std::istream& in)
      {
	in >> static_cast<JObjectID&>(*this) >> pmtType >> serialNumber >> address;
      }

      int pmtType;        //!< PMT type
      int serialNumber;   //!< PMT serial number
      int address;        //!< address
    };


    /**
     * Auxiliary class for OM cluster parameters
     */
    class OM_cluster : 
      public JObjectID, 
      public std::vector<int> 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      OM_cluster(std::istream& in)
      {
	int n;
	  
	in >> static_cast<JObjectID&>(*this) >> string_id >> z >> n;
	  
	this->resize(n);

	iterator p = this->begin();

	for ( ; n != 0; --n, ++p)
	  in >> *p;
      }


      int    string_id;   //!< string identifier
      double z;           //!< z position
    };


    /**
     * Auxiliary class for OM cluster data
     */
    class OM_cluster_data : 
      public JObjectID
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      OM_cluster_data(std::istream& in)
      {
	in >> static_cast<JObjectID&>(*this) >> x >> y >> z >> theta >> phi >> psi;
      }


      double x;           //!< x position
      double y;           //!< y position
      double z;           //!< z position
      double theta;       //!< Euler angle
      double phi;         //!< Euler angle
      double psi;         //!< Euler angle
    };

    /**
     * Auxiliary class for OM position
     */
    class OM_position : 
      public JObjectID 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      OM_position(std::istream& in)
      {
	in >> static_cast<JObjectID&>(*this) >> x >> y >> z >> theta >> phi;
      }


      double x;           //!< x position
      double y;           //!< y position
      double z;           //!< z position
      double theta;       //!< zenit angle of orientation
      double phi;         //!< azimuth angle of orientation
    };


    /**
     * Auxiliary class for string parameters
     */
    class String : 
      public JObjectID 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      String(std::istream& in)
      {
	in >> static_cast<JObjectID&>(*this) >> x >> y >> z >> tOffset >> twist >> twistRate;
      }


      double x;           //!< x position
      double y;           //!< y position
      double z;           //!< z position
      double tOffset;     //!< time offset
      double twist;       //!< twist offset
      double twistRate;   //!< twist rate
    };


    /**
     * Auxiliary class for LCM logic parameters
     */
    class LCM_logic : 
      public JObjectID, 
      public std::vector<int> 
    {
    public:
      /**
       * Constructor.
       *
       * \param  in        input stream
       */
      LCM_logic(std::istream& in)
      {
	int n;
	  
	in >> static_cast<JObjectID&>(*this) >> string_id >> lcm_id >> n;

	this->resize(n);

	iterator p = this->begin();
	  
	for ( ; n != 0; --n, ++p)
	  in >> *p;
      }

      int string_id;      //!< string identifier
      int lcm_id;         //!< LCM identifier
    };


    /**
     * Auxiliary class for LCM logic parameters
     */
    class LCM_reverse_logic : 
      public JObjectID 
    {
    public:
      /**
       * Default constructor.
       */
      LCM_reverse_logic() :
	JObjectID(),
	lcm_id(-1),
	pmt_id(-1)
      {}


      /**
       * Constructor.
       *
       * \param  __id      identifier
       * \param  __lcm_id  LCM identifier
       * \param  __pmt_id  PMT identifier
       */
      LCM_reverse_logic(const int __id, 
			const int __lcm_id, 
			const int __pmt_id) :
	JObjectID(__id),
	lcm_id   (__lcm_id),
	pmt_id   (__pmt_id)
      {}


      int lcm_id;         //!< LCM identifier
      int pmt_id;         //!< PMT identifier
    };


  public:
    /**
     * Default constructor.
     */
    JMonteCarloDetector() :
      JDetector(),
      use_logic(false)
    {}


    /**
     * Constructor.
     *
     * \param  useLogic    apply logic data in file
     */
    JMonteCarloDetector(const bool useLogic) :
      JDetector(),
      use_logic(useLogic)
    {}



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

      detector.clear();


      vector<OM>              OMaddress;
      vector<OM_cluster>      OMcluster;
      vector<OM_position>     OMposition;
      vector<String>          DetectorString;
      vector<OM_cluster_data> OMclusterdata;
      vector<LCM_logic>       LCMlogic;


      for (string key, buffer; in >> key; ) {

	if      (key == "OM:")

	  OMaddress.push_back(OM(in));

	else if (key == "OM_cluster:")

	  OMcluster.push_back(OM_cluster(in));

	else if (key == "OM_position:")

	  OMposition.push_back(OM_position(in));

	else if (key == "string:")

	  DetectorString.push_back(String(in));

	else if (key == "OM_cluster_data:")

	  OMclusterdata.push_back(OM_cluster_data(in));

	else if (key == "LCM_logic:")

	  LCMlogic.push_back(LCM_logic(in));

	else

	  getline(in, buffer);
      }


      // invert LCM_Logic data for fast lookup

      vector<LCM_reverse_logic> LCMreverselogic;

      for (vector<LCM_logic>::const_iterator i  = LCMlogic.begin(); i != LCMlogic.end(); ++i) {
	
	const int lcm_id = i->getID();

	for (LCM_logic::const_iterator j = i->begin(); j != i->end(); ++j) {

	  const int id     = *j;
	  const int pmt_id = distance(i->begin(),j) + 1;

	  LCMreverselogic.push_back(LCM_reverse_logic(id, lcm_id, pmt_id));
	}
      }


      // sort data according to identifier

      sort(OMaddress      .begin(), OMaddress      .end());
      sort(OMcluster      .begin(), OMcluster      .end());
      sort(OMposition     .begin(), OMposition     .end());
      sort(DetectorString .begin(), DetectorString .end());
      sort(LCMreverselogic.begin(), LCMreverselogic.end());


      vector<OM_cluster>     ::const_iterator   OMclusterIterator;
      vector<String>         ::const_iterator   StringIterator;
      vector<int>            ::const_iterator   addressIterator;
      vector<OM>             ::const_iterator   OMIterator;
      vector<OM_position>    ::const_iterator   OMpositionIterator;
      vector<OM_cluster_data>::const_iterator   OMclusterdataIterator;


      for (OMclusterIterator  = OMcluster.begin(); 
	   OMclusterIterator != OMcluster.end(); 
	   OMclusterIterator++) {
	
	StringIterator = find(DetectorString.begin(),
			      DetectorString.end(),
			      OMclusterIterator->string_id);
	
	if (StringIterator != DetectorString.end()) {
	  
	  const JLocation location(StringIterator->getID(), -1);

	  for (addressIterator  = OMclusterIterator->begin(); 
	       addressIterator != OMclusterIterator->end(); 
	       addressIterator++) {
	    
	    OMIterator = find(OMaddress.begin(),
			      OMaddress.end(),
			      *addressIterator);
	    
	    if (OMIterator != OMaddress.end()) {
	      
	      OMpositionIterator = find(OMposition.begin(),
					OMposition.end(),
					OMIterator->address);
	      
	      if (OMpositionIterator != OMposition.end()) {
		
		// rotate string
		
		double ct = cos(StringIterator->twist * PI / 180.0);
		double st = sin(StringIterator->twist * PI / 180.0);
		
		double x  =  ct * OMpositionIterator->x - st * OMpositionIterator->y;
		double y  =  st * OMpositionIterator->x + ct * OMpositionIterator->y;

		JPosition3D  pos(StringIterator->x + x,
				 StringIterator->y + y,
				 StringIterator->z + OMpositionIterator->z + OMclusterIterator->z);
		
		JDirection3D dir(JAngle3D(OMpositionIterator->theta, 
					  StringIterator->twist * PI / 180.0 + OMpositionIterator->phi));
		
		int lcm_id  = OMclusterIterator->getID();
		int pmt_id  = OMIterator->address;
		
		if (detector.use_logic) {
		  
		  vector<LCM_reverse_logic>::const_iterator p = find(LCMreverselogic.begin(), 
								     LCMreverselogic.end(), 
								     OMIterator->getID());
		  
		  if (p != LCMreverselogic.end()) {
		    lcm_id = p->lcm_id;
		    pmt_id = p->pmt_id;
		  }
		}
		
		detector.put(lcm_id, pmt_id - 1, location, JPMT(OMIterator->getID(), JAxis3D(pos, dir)));
	      }
	    }
	  }
	}
      }


      // 'event' data

      for (OMclusterdataIterator  = OMclusterdata.begin(); 
	   OMclusterdataIterator != OMclusterdata.end(); 
	   OMclusterdataIterator++) {
	  
	const JPosition3D    P(OMclusterdataIterator->x,
			       OMclusterdataIterator->y,
			       OMclusterdataIterator->z);
	
	const JEulerMatrix3D R(JEulerAngle3D(OMclusterdataIterator->theta,
					     OMclusterdataIterator->phi,
					     OMclusterdataIterator->psi));
	
	OMclusterIterator = find(OMcluster.begin(),
				 OMcluster.end(),
				 OMclusterdataIterator->getID());

	if (OMclusterIterator != OMcluster.end()) {

	  const JLocation location(OMclusterIterator->string_id, -1);

	  for (addressIterator  = OMclusterIterator->begin(); 
	       addressIterator != OMclusterIterator->end(); 
	       addressIterator++) {

	    OMIterator = find(OMaddress.begin(),
			      OMaddress.end(),
			      *addressIterator);
	      
	    if (OMIterator != OMaddress.end()) {
		
	      OMpositionIterator = find(OMposition.begin(),
					OMposition.end(),
					OMIterator->address);
		
	      if (OMpositionIterator != OMposition.end()) {
		  
		JPosition3D  pos(OMpositionIterator->x,
				 OMpositionIterator->y,
				 OMpositionIterator->z);
		  
		JDirection3D dir(JAngle3D(OMpositionIterator->theta, 
					  OMpositionIterator->phi));
		  		  
		pos.transform(R);
		dir.transform(R);
		  
		pos += P;
		  
		int lcm_id  = OMclusterIterator->getID();
		int pmt_id  = OMIterator->address;

		if (detector.use_logic) {

		  vector<LCM_reverse_logic>::const_iterator p = find(LCMreverselogic.begin(), 
								     LCMreverselogic.end(), 
								     OMIterator->getID());
		  
		  if (p != LCMreverselogic.end()) {
		    lcm_id = p->lcm_id;
		    pmt_id = p->pmt_id;
		  }
		}

		detector.put(lcm_id, pmt_id - 1, location, JPMT(OMIterator->getID(), JAxis3D(pos, dir)));
	      }
	    }
	  }
	}
      }

      // set floor identifier

      sort(detector.begin(), detector.end(), JMonteCarloDetector::compare);

      int string = -1;
      int floor  = -1;
      
      for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
	
	if (module->getString() != string) {
	  string = module->getString();
	  floor  = 0;
	}

	module->setLocation(JLocation(string,++floor));
      }
      
      return in;
    }


    /**
     * Set usage of logic.
     *
     * \param  useLogic    apply logic data in file
     */
    void setUseLogic(const bool useLogic) 
    {
      use_logic = useLogic;
    }


  protected:
    /**
     * Set PMT.
     * 
     * \param  moduleID    module identifier
     * \param  pmtAddress  PMT address
     * \param  location    module location
     * \param  pmt         JPMT
     * \return             this JDetector
     */
    JDetector& put(const int        moduleID,
		   const int        pmtAddress,
		   const JLocation& location,
		   const JPMT&      pmt)
    {
      JDetector::iterator p = std::lower_bound(this->begin(), this->end(), moduleID);

      if (p == this->end() || p->getID() != moduleID) {
        
        JModule module(moduleID, location);
        
        module.resize(pmtAddress + 1);

        module[pmtAddress] = pmt;

	module.setPosition(pmt.getPosition());
        
        this->insert(p, module);
        
      } else {
        
        // average position
        
	int n = 0;

	for (JModule::const_iterator i = p->begin(); i != p->end(); ++i) {
	  if (i->getID() != getUndefinedObjectID().getID()) {
	    ++n;
	  }
	}

	if (pmtAddress + 1 > (int) p->size()) {
	  p->resize(pmtAddress + 1);
	}

        JPosition3D& P = static_cast<JPosition3D&>(*p);

	P.mul(n);

	if (p->at(pmtAddress).getID() != getUndefinedObjectID().getID())
	  P.sub(p->at(pmtAddress).getPosition());
	else
	  ++n;

	P.add(pmt.getPosition());
	P.div(n);

	p->at(pmtAddress) = pmt;
      }

      return *this;
    }


    /**
     * Binary search method.
     *
     * \param  __begin     begin of data
     * \param  __end       end   of data
     * \param  id          target value
     * \return             pointer to corresponding element 
     */
    template<class T>
    static T find(T __begin, T __end, const int id) 
    {
      T i = std::lower_bound(__begin, __end, id);

      if (i != __end && *i != id)
	return __end;
      else
	return i;
    }


    /**
     * Module comparator.\n
     * The comparison is based on:
     *    -# string identifier;
     *    -# z-position.
     *
     * \param  first       first  module
     * \param  second      second module
     * \return             true if first module before second; else false
     */
    static bool compare(const JModule& first, const JModule& second)
    {
      if (first.getString() == second.getString())
	return first.getZ()      < second.getZ();
      else
	return first.getString() < second.getString();
    }

  private:
    bool use_logic;
  };
}

#endif