#include <ICOMETLog.hxx>
#include <IOADatabase.hxx>
#include <IGeomIdManager.hxx>
#include <ICOMETInput.hxx>
#include "IBeamlineCoordinates.hxx"
#include "IBeamlineParameterExtractor.hxx"
#include "IBeamlineCoordinatesPainter.hxx"
#include <HEPUnits.hxx>

#include <TRandom3.h>
#include <TPolyMarker3D.h>
#include <TPointSet3D.h>
#include <TVector3.h>
#include <TGeoNode.h>
#include <TGeoVolume.h>
#include <TGeoShape.h>
#include <TGeoBBox.h>
#include <TGeoManager.h>

#include <cstdlib>
#include <ios>
#include <sstream>
#include <string>
#include <map>
#include <utility>
#include <vector>
#include <fstream>

bool operator<(const TVector3& left, const TVector3& right);

namespace COMET {
namespace GeomTools {

    class PointSelector{
        public:
        PointSelector(int totalPoints): fTotalPoints(totalPoints), 
                                        fPointIndex(0){}
        virtual ~PointSelector(){}

        virtual bool GetNextPoint(TVector3& setPoint, bool verb=true) {
            if (fTotalPoints >= 10 && verb){
                int tenPerCent=fTotalPoints/10;
                // Print progress if needed
                if (fPointIndex%tenPerCent==0)  
                    std::cout<<" "<<fPointIndex/tenPerCent*10
                                  <<"% complete"<<std::endl;
            }
            // Get the next point
            if (fPointIndex < fTotalPoints){
                setPoint = GeneratePoint();
                ++fPointIndex;
                return true;
            } else return false;
        };

        virtual TVector3 GeneratePoint()=0;

        protected:
        int fTotalPoints;
        int fPointIndex;
    };


    struct PointDetails{
        std::string material;
        std::string full_path;
        bool inOverlap;
        PointDetails():material("unknown"),full_path("unknown"),inOverlap(false){}
        PointDetails(TGeoNode* volume){
            material=volume->GetMedium()->GetMaterial()->GetName();
            full_path=volume->GetVolume()->GetTitle();
            inOverlap=volume->IsOverlapping();
        }
        bool operator!=(const PointDetails& rhs)const{
            return  material!=rhs.material 
                || full_path!=rhs.full_path
                || inOverlap!=rhs.inOverlap;
        }
        std::string toString()const{
            std::stringstream stream;
            stream<<"(Path="<< full_path <<", ";
            stream<<"Material="<< material <<", ";
            stream<<"Overlap="<< std::boolalpha << inOverlap <<")";
            return stream.str();
        }
    };

    struct PointPair{
        enum Severity_t{
            kNoDiff=0,
            kPathName=1<<0,
            kMaterial=1<<1,
            kOverlap=1<<2,
            kMaxSeverity=(1<<3)-1
        };
        PointDetails first;
        PointDetails second;
        PointPair():first(),second(){}
        PointPair(const PointDetails& one, const PointDetails& two):first(one),second(two){}
        std::string toString()const{
            std::stringstream stream;
            if(first.material!=second.material) {
                stream<<"Material: "<<first.material <<" --> "<<second.material<<" ";
            }
            if(first.full_path!=second.full_path) {
                stream<<"Path: "<<first.full_path <<" --> "<<second.full_path<<" ";
            }
            if(first.inOverlap!=second.inOverlap) {
                stream<<"Overlap: "<<std::boolalpha<<first.inOverlap <<" --> "<<second.inOverlap;
            }
            return stream.str();
        }
        int GetDifferenceSeverity()const{
            int output=kNoDiff;
            if(first.material!=second.material) output|=kMaterial;
            if(first.full_path!=second.full_path) output |=kPathName;
            if(first.inOverlap!=second.inOverlap) output |=kOverlap;
            return output;
        }
        std::string SeverityString()const{ return SeverityString(GetDifferenceSeverity()); } 
        static std::string SeverityString(int severity){
            std::string output;
            if( (severity&kMaterial)) output+="Material";
            if( (severity&kPathName)) output+="Path";
            if( (severity&kOverlap )) output+="Overlap";
            return output;
        }
    };

    typedef std::map<TVector3,PointPair> PointDetailsList;

    struct PointDiffLists: private std::vector<PointDetailsList>{
        PointDiffLists():std::vector<PointDetailsList>(PointPair::kMaxSeverity){};
        void AddDiff(const TVector3& place, const PointPair& diff){
            int severity=diff.GetDifferenceSeverity();
            at(severity)[place]=diff;
        }
        const PointDetailsList& GetList(int severity)const{return at(severity);}
        size_t size()const{
            size_t sz=0;
            for(const_iterator i_list=this->begin(); i_list!=this->end(); ++i_list)
                    sz+=i_list->size();
            return sz;
        }
        using std::vector<PointDetailsList>::begin;
        using std::vector<PointDetailsList>::end;
        using std::vector<PointDetailsList>::iterator;
        using std::vector<PointDetailsList>::const_iterator;
        int SeverityColour(int severity)const{
            switch(severity){
                case PointPair::kNoDiff:   return kBlack;
                case PointPair::kPathName: return kMagenta;
                case PointPair::kMaterial: return kOrange;
                case PointPair::kMaterial|PointPair::kPathName: return kRed;
                case PointPair::kOverlap:  return kGreen;
                case PointPair::kOverlap |PointPair::kMaterial: return kCyan;
                case PointPair::kOverlap |PointPair::kPathName: return kCyan;
                case PointPair::kOverlap |PointPair::kMaterial|PointPair::kPathName: return kBlue;
            }
            return kOrange;
        }
    };

    typedef std::vector<std::string> SubStrings;

    class RandomDaughterPoints: public PointSelector{
        public:
            RandomDaughterPoints(int totalPoints, TGeoManager* geom_one):
                PointSelector(totalPoints),
                fGeoManager(geom_one),
                fNode(NULL),
                fVolume(NULL) {
                    fNode=geom_one->GetCurrentNode();
                    fVolume=geom_one->GetCurrentVolume();

                    const TGeoShape *shape = fVolume->GetShape();
                    TGeoBBox *box = (TGeoBBox *)shape;
                    fDx = box->GetDX();
                    fDy = box->GetDY();
                    fDz = box->GetDZ();
                    fOx = (box->GetOrigin())[0];
                    fOy = (box->GetOrigin())[1];
                    fOz = (box->GetOrigin())[2];
                }
            virtual TVector3 GeneratePoint(){
                double xyz[3];
                bool point_found=false;
                do{
                    xyz[0] = fOx-fDx+2*fDx*gRandom->Rndm();
                    xyz[1] = fOy-fDy+2*fDy*gRandom->Rndm();
                    xyz[2] = fOz-fDz+2*fDz*gRandom->Rndm();
                    fGeoManager->PushPath();
                    fGeoManager->SetCurrentPoint(xyz);
                    TGeoNode* node = fGeoManager->FindNode();
                    fGeoManager->PopPath();
                    if (!node) continue;
                    if(node==fNode) continue;
                    point_found=true;
                }while(!point_found);
                double master[3];
                fGeoManager->LocalToMaster(xyz,master);

                return TVector3(master[0],master[1],master[2]); 
            }
        private:
            TGeoManager* fGeoManager;
            TGeoNode* fNode;
            TGeoVolume* fVolume;
            double fDx, fDy, fDz, fOx, fOy, fOz;
    };

    class PointsAlongBeamLine: public PointSelector{
        public:
            PointsAlongBeamLine(int totalPoints, TGeoManager* geom_one):
                PointSelector(totalPoints),
                fGeoManager(geom_one){
                    COMET::IBeamlineParameterExtractor::Load(fBeamline,
                                                             fGeoManager);
                    fBeamline.Initialise();
                    fBeamline.Print();
                }
            virtual TVector3 GeneratePoint(){
                // Get the length point here
                double aFract  = fPointIndex/(double)fTotalPoints;
                double aLength = aFract * fBeamline.GetTotalLength();
                // Return the global coordinate
                return fBeamline.ToGlobalCoords(aLength, 0, 0) * (1/unit::cm);
            }
        private:
            TGeoManager* fGeoManager;
            COMET::IBeamlineCoordinates fBeamline;
    };

    class PointsFromFile: public PointSelector{
        /// A container of points to read file into
        typedef std::vector<TVector3> PointsArray;
        public:
            /// Constructor from a file name
            PointsFromFile(int totalPoints, std::string& fileName):
                PointSelector(totalPoints),
                fFileName(fileName),
                fPointsArray(){
                    // Open the input file
                    std::ifstream input(fFileName.c_str());
                    // Cache the line data
                    std::string line_data;
                    // Cache the line data values
                    double xPos, yPos, zPos;
                    // Stream the line data
                    std::stringstream line_stream(line_data);
                    while (std::getline(input, line_data)){
                        // Set the string stream
                        line_stream.str(line_data);
                        if (line_stream >> xPos >> yPos >> zPos){
                            fPointsArray.push_back(TVector3(xPos,yPos,zPos)*(1/unit::cm));
                        }
                        // Clear the string stream
                        line_stream.str("");
                        line_stream.clear();
                    }
                    // Set the total number of points to the number of points in 
                    // the file if there are less points in the file than in the
                    // provided max, totalPoints
                    fTotalPoints = std::min(fTotalPoints, 
                                            (int)fPointsArray.size());
                    // Set the first point to the beginning
                    fThisPoint = fPointsArray.begin();
                }
            /// Read the point from the vector
            virtual TVector3 GeneratePoint(){
                TVector3 returnVal(*fThisPoint);
                fThisPoint++; 
                return returnVal;
            }
        private:
            std::string fFileName;
            PointsArray fPointsArray;
            PointsArray::iterator fThisPoint;
    };

    void SplitString(const std::string& inputString, SubStrings& outStrings, char delim=' ');
    TGeoManager* FetchGeometry(COMET::ICOMETInput*,std::string,std::string test_volume="");
    PointSelector* BuildPointSelector(const std::string&, int, TGeoManager*);
    bool ScanFirstGeometry(TGeoManager*,PointSelector*,PointDetailsList&);
    bool ScanSecondGeometry(TGeoManager*,PointDetailsList&,PointDiffLists&);
    void PrintDiff(const PointDiffLists&,bool printEverything);
    bool PrintDetailList(PointDetailsList& outList, 
                         bool printFirst=true, 
                         bool printSecond=true);
    bool SaveEverything(const std::string& output_file, const PointDiffLists& pointDiffs, const PointDetailsList& pointDetails);
    void SavePoints(const PointDetailsList& list, int style,int colour,std::string name="");
    TGeoNode* GetNode(TGeoManager* geometry, const TVector3& position);
}
}