//==========================================================================
//  AIDA Detector description implementation
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author     : A. Sailer
//
//==========================================================================
#ifndef DDCAD_UTILITIES_H
#define DDCAD_UTILITIES_H

#include <vector>
#include <cmath>
#include <limits>

#include <TGeoTessellated.h>
#include <TGeoVector3.h>

/// Namespace for the AIDA detector description toolkit
namespace dd4hep {

  /// Namespace for implementation details of the AIDA detector description toolkit
  namespace cad  {

#if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
    inline std::string streamFacet(TGeoFacet const& facet,
                                   TGeoTessellated const& shape) {
      using ::operator<<;
      std::stringstream str;
      str << "{";
      for (int i = 0; i < facet.GetNvert(); ++i) {
        str << shape.GetVertex(facet[i]);
        if (i != facet.GetNvert() - 1)
          str << ", ";
      }
      str << "}";
      return str.str();
    }
#else
    inline std::string streamFacet(TGeoFacet const& facet,
                                   TGeoTessellated const& /* shape */) {
      using ::operator<<;
      std::stringstream str;
      str << facet;
      return str.str();
    }
#endif
    
    inline std::string streamVertices(ROOT::Geom::Vertex_t const& v1,
					    ROOT::Geom::Vertex_t const& v2,
					    ROOT::Geom::Vertex_t const& v3) {
      using ::operator<<;
      std::stringstream str;
      str << "{" << v1 << ", " << v2 <<  ", " << v3 << "}";
      return str.str();
    }
    // Determine if the facet is degenerated by calculating its determinant
    inline bool facetIsDegenerated(std::vector<ROOT::Geom::Vertex_t> const& vertices){
      using ROOT::Geom::Vertex_t;
      // Compute normal using non-zero segments
      constexpr double kTolerance = 1.e-20;
      bool degenerated = true;
      Vertex_t normal;
      int fNvert = int(vertices.size());
      for (int i = 0; i < fNvert - 1; ++i) {
        Vertex_t e1 = vertices[i + 1] - vertices[i];
        if (e1.Mag2() < kTolerance)
          continue;
        for (int j = i + 1; j < fNvert; ++j) {
          Vertex_t e2 = vertices[(j + 1) % fNvert] - vertices[j];
          if (e2.Mag2() < kTolerance)
            continue;
          normal = Vertex_t::Cross(e1, e2);
          // e1 and e2 may be colinear
          if (normal.Mag2() < kTolerance)
            continue;
          normal.Normalize();
          degenerated = false;
          break;
        }
        if (!degenerated)
          break;
      }
      return degenerated;
    }
  }
}
#endif
