#include <cmath>
#include <algorithm>
#include <memory>
#include <iostream>
#include <tut.h>

#include <TH2F.h>
#include <TFile.h>
#include <TRandom.h>
#include <TGeoManager.h>

#include <IGeomIdManager.hxx>
#include <IOADatabase.hxx>
#include <oaEvent_version.h>

#include "IGeomInfo.hxx"

// Namespace surrounding all tut testing code
namespace tut {
    struct baseGeomInfo {
        baseGeomInfo() {
            // Run before each test.
            COMET::IOADatabase::Get().RegisterGeometryLookup(NULL);
            COMET::IOADatabase::Get().RegisterAlignmentLookup(NULL);
            if (!gGeoManager) {
                COMET::IOADatabase::Get().GeomId().ReadGeometry(
                    COMET::ISHAHashValue());
            }
            else {
                COMET::IOADatabase::Get().GeomId().ResetGeometry();
            }
            COMET::IOADatabase::Get().ClearGeometryCallbacks();
            COMET::IOADatabase::Get().SetGeometryOverride(
                COMET::IOADatabase::Get().GeomId().GetHash());
        }
        ~baseGeomInfo() {
            // Run after each test.
            COMET::IOADatabase::Get().RegisterGeometryLookup(NULL);
            COMET::IOADatabase::Get().RegisterAlignmentLookup(NULL);
            COMET::IOADatabase::Get().ClearGeometryCallbacks();
        }
    };

    // Declare the test
    typedef test_group<baseGeomInfo>::object testGeomInfo;
    test_group<baseGeomInfo> groupGeomInfo("GeomInfo");

    // Test the default constructor and destructor. 
    template<> template<>
    void testGeomInfo::test<1> () {
        COMET::IGeomInfo::Get();
    }

    /// Check the node id and the node position match.
    template<> template<>
    void testGeomInfo::test<2> () {
        int geomIdCount = 0;
        for (int i=0; i<100000; ++i) {
            TVector3 p0(gRandom->Uniform(-100*unit::cm,100*unit::cm),
                         gRandom->Uniform(-100*unit::cm,100*unit::cm),
                         gRandom->Uniform(-300*unit::cm,300*unit::cm));
	    COMET::IGeometryId n1 
                = COMET::IGeomInfo::Get().GeomId(p0.X(), p0.Y(), p0.Z());
            if (n1 == COMET::IGeometryId()) continue;
            ++geomIdCount;
            TVector3 p1 = COMET::IGeomInfo::Get().NodePosition(n1);
	    COMET::IGeometryId n2
                = COMET::IGeomInfo::Get().GeomId(p1.X(), p1.Y(), p1.Z());
            TVector3 p2 = COMET::IGeomInfo::Get().NodePosition(n2);
            if (n1 != n2) {
                std::string name1 = COMET::IGeomInfo::Get().NodeName(n1);
                std::string name2 = COMET::IGeomInfo::Get().NodeName(n2);
                if (name1.find(name2) != std::string::npos) continue;
                if (name2.find(name1) != std::string::npos) continue;
            }
            ensure_equals("Node position matches node number",n1,n2);
        }
        ensure_greaterthan("Geometry objects found", geomIdCount, 0);
    }

    /// Check that scintillator bars are found.
    template<> template<>
    void testGeomInfo::test<3> () {
        ensure_greaterthan("P0D bars exist",
                           COMET::IGeomInfo::P0D().GetBars().size(),
                           (size_t) 0);
	/*
        ensure_greaterthan("FGD bars exist",
                           COMET::IGeomInfo::FGD().GetBars().size(),
                           (size_t) 0);
        ensure_greaterthan("ECAL Modules exist",
                           COMET::IGeomInfo::ECAL().GetECALModules().size(),
                           (size_t) 0);
        for (std::vector<COMET::IECALModule*>::const_iterator m 
                 = COMET::IGeomInfo::ECAL().GetECALModules().begin();
             m != COMET::IGeomInfo::ECAL().GetECALModules().end();
             ++m) {
            ensure_greaterthan("ECAL module bars exist",
                               (*m)->GetBarMap().size(),
                               (size_t) 0);
        }
	*/
        ensure_greaterthan("SMRD Modules exist",
                           COMET::IGeomInfo::SMRD().GetBars().size(),
                           (size_t) 0);

        /// Check that a P0D bar can be found.
        ensure_greaterthan("P0D X Bars exist", 
                           COMET::IGeomInfo::P0D().GetXBars().size(), 
                           (size_t) 0);
        COMET::IGeometryId p0dNode =  COMET::IGeomInfo::P0D().GetXBars().front();
        const COMET::IScintBarGeom& p0dBar 
            = COMET::IGeomInfo::Get().GetScintBar(p0dNode);
        ensure_equals("P0D bar nodes match", p0dNode, p0dBar.GetGeomId());
	/*
        /// Check that an FGD bar can be found.
        COMET::IGeometryId fgdNode = COMET::IGeomInfo::FGD().GetBars().begin()->first;
        const COMET::IScintBarGeom& fgdBar
            = COMET::IGeomInfo::Get().GetScintBar(fgdNode);
        ensure_equals("FGD bar nodes match", fgdNode, fgdBar.GetGeomId());

        for (std::vector<COMET::IECALModule*>::const_iterator m 
                 = COMET::IGeomInfo::ECAL().GetECALModules().begin();
             m != COMET::IGeomInfo::ECAL().GetECALModules().end();
             ++m) {
            COMET::IGeometryId ecalNode = (*m)->GetBarMap().begin()->first;
            const COMET::IScintBarGeom& ecalBar
                = COMET::IGeomInfo::Get().GetScintBar(ecalNode);
            ensure_equals("ECAL bar nodes match", 
                          ecalNode, ecalBar.GetGeomId());
        }
	*/
        /// Check that an SMRD bar can be found.
        COMET::IGeometryId smrdNode
            = COMET::IGeomInfo::SMRD().GetBars().begin()->first;
        const COMET::IScintBarGeom& smrdBar
            = COMET::IGeomInfo::Get().GetScintBar(smrdNode);
        ensure_equals("SMRD bar nodes match", smrdNode, smrdBar.GetGeomId());

        /// Check that the exception is thrown if an incorrect node id is
        /// passed.  Node 2 is probably the hall, but is definately not a
        /// scintillating bar.
        try {
            COMET::IGeomInfo::Get().GetScintBar(COMET::IGeometryId(2));
            ensure("Exception should be thrown",false);
        }
        catch (COMET::EVolumeNotFound& e) { }
    }

    /// Check the coordinate system translation methods.
    template<> template<>
    void testGeomInfo::test<4> () {
      /*
        // Find a bar that is part of the front FGD face.
        ensure_greaterthan("Enough FGD bars",
                           COMET::IGeomInfo::FGD().GetBars().size(),
                           (size_t) 200);

        // Get an FGD bar that should be in about the center of the first
        // plane.  This uses an FGD bar since they are rectangular boxes.  The
        // P0D uses triangular bars and complicates the tests.
        COMET::IFGDGeom::BarMap::const_iterator bar 
            = COMET::IGeomInfo::FGD().GetBars().begin();
        for (int i=0; i<100; ++i) ++bar;
        COMET::IGeometryId nodeId = bar->first;

        // Make sure that it is in the list of scintillating bars.  This will
        // thrown an COMET::EVolumeNotFound exception if there is a problem and
        // it will be caught by tut.
        COMET::IGeomInfo::Get().GetScintBar(nodeId);

        // Get the size of the bar in the master coordinate system.
        TVector3 size = COMET::IGeomInfo::Get().NodeSize(nodeId);
        ensure_distance("Bar width", size.X(), 4.8*unit::mm, 0.5*unit::mm);
        ensure_distance("Bar height", size.Y(), 4.8*unit::mm, 0.5*unit::mm);
        ensure_greaterthan("Bar length", size.Z(), 900*unit::mm);

        // Get extent of the bar in the master coordinate system
        TVector3 extent = COMET::IGeomInfo::Get().NodeExtent(nodeId);
        ensure_distance("Bar X extent", extent.X(), size.X(), 0.001);
        ensure_distance("Bar Y extent", extent.Y(), size.Z(), 0.001);
        ensure_distance("Bar Z extent", extent.Z(), size.Y(), 0.001);

        // Check that random points within "extent" of the node position are
        // assigned to this node.
        TVector3 position = COMET::IGeomInfo::Get().NodePosition(nodeId);
        for (int i=0; i<1000; ++i) {
            // Generate a trial position.
            TVector3 offset(gRandom->Uniform(-0.99*extent.X(), 
                                             0.99*extent.X()),
                            gRandom->Uniform(-0.99*extent.Y(),
                                             0.99*extent.Y()),
                            gRandom->Uniform(-0.99*extent.Z(),
                                             0.99*extent.Z()));
            TVector3 trial = position + offset;
            COMET::IGeometryId trialNode = COMET::IGeomInfo::Get().GeomId(trial.X(), 
                                                        trial.Y(), 
                                                        trial.Z());
            ensure_equals("Trial point is inside bar", trialNode, nodeId);
            TVector3 local 
                = COMET::IGeomInfo::Get().MasterToLocal(nodeId,
                                                     trial.X(),
                                                     trial.Y(),
                                                     trial.Z());
            ensure_lessthan("Local X coordinate inside bar",
                            std::abs(local.X()), size.X());
            ensure_lessthan("Local Y coordinate inside bar",
                            std::abs(local.Y()), size.Y());
            ensure_lessthan("Local Z coordinate inside bar",
                            std::abs(local.Z()), size.Z());
            
            ensure_distance("Offset X and Local X match",
                            offset.X(), local.X(), 0.001);
            ensure_distance("Offset Y and Local Z match",
                            offset.Y(), local.Z(), 0.001);
            ensure_distance("Offset -Z and Local Y match",
                            -offset.Z(), local.Y(), 0.001);

        }

        // Check that random points not within "extent" of the node position
        // are not assigned to this node.
        for (int i=0; i<1000; ++i) {
            // Generate a trial position.
            TVector3 offset(gRandom->Uniform(-30*extent.X(), 30*extent.X()),
                            gRandom->Uniform(-30*extent.Y(), 30*extent.Y()),
                            gRandom->Uniform(-30*extent.Z(), 30*extent.Z()));
            if (std::abs(offset.X()) < extent.X()
                && std::abs(offset.Y()) < extent.Y()
                && std::abs(offset.Z()) < extent.Z()) continue;
            TVector3 trial = position + offset;
            COMET::IGeometryId trialNode = COMET::IGeomInfo::Get().GeomId(trial.X(), 
                                                        trial.Y(), 
                                                        trial.Z());
            ensure("Trial point is outside bar", trialNode != nodeId);
            TVector3 local 
                = COMET::IGeomInfo::Get().MasterToLocal(nodeId,
                                                     trial.X(),
                                                     trial.Y(),
                                                     trial.Z());
            ensure("Local coordinate outside bar",
                   (size.X() < std::abs(local.X()))
                   || (size.Y() < std::abs(local.Y()))
                   || (size.Z() < std::abs(local.Z())));

            ensure_distance("Offset X and Local X match",
                            offset.X(), local.X(), 0.001);
            ensure_distance("Offset Y and Local Z match",
                            offset.Y(), local.Z(), 0.001);
            ensure_distance("Offset -Z and Local Y match",
                            -offset.Z(), local.Y(), 0.001);

        }
      */
    }

};