#include <iostream>
#include <cmath>
#include <algorithm>

#include <TGeoManager.h>
#include <TGeoNode.h>
#include <TGeoArb8.h>
#include <TVector3.h>
#include <TRandom.h>

#include <HEPUnits.hxx>
#include <IGeometryId.hxx>
#include <COMETGeomId.hxx>
#include <COMETGeomIdDef.hxx>

#include "IGeomVisitor.hxx"
#include "IP0DGeom.hxx"

#include "IOADatabase.hxx"
#include "IGeomIdManager.hxx"

COMET::IP0DGeom::IP0DGeom() 
    : IGeomBase(), fNeighborsXZ(NULL), fNeighborsYZ(NULL) {}

COMET::IP0DGeom::~IP0DGeom() {
    this->Clear();
}

namespace COMET {
/// For information see the discription of the base class 'IGeomVisitor'
    class IP0DGeomVisitor: public COMET::IGeomVisitor {
    public:
        explicit IP0DGeomVisitor(COMET::IP0DGeom* geom) {
            fP0DGeom = geom;
        }
        bool VisitNode(int d, const std::string& name, const TGeoNode* node);
        std::vector<COMET::IP0DGeom::Bar*> fXBars;
        std::vector<COMET::IP0DGeom::Bar*> fYBars;
    private:
        COMET::IP0DGeom* fP0DGeom;
    };
};

bool COMET::IP0DGeomVisitor::VisitNode(int, 
                                    const std::string& name, 
                                    const TGeoNode*) {
    // If the node is in the wrong sub-detector, skip the entire heirarchy.
    if (name.find("/TPC") != std::string::npos) return false;
    if (name.find("/FGD") != std::string::npos) return false;
    if (name.find("/DsECal") != std::string::npos) return false;

    // If the node isn't deep enough to be in the P0D, skip this node.
    if (name.find("/P0D_") == std::string::npos) return true;

    // Check if we have found P0D top volume
    if( name.find("/P0D_" ) != std::string::npos 
        && fP0DGeom->fModules.size() < 1) {      
      COMET::IGeomModuleBase* p0d = new IGeomModuleBase(name);
      
      p0d->SetX0(25*unit::cm);  // THIS APPROXIMATES WATER TARGET.
      fP0DGeom->fModules.push_back(p0d);
    }
    COMET::IGeometryId bgeom;
    COMET::IOADatabase::Get().GeomId().FindGeometryId(bgeom);

    if (name.find("/P0Dule_") != std::string::npos
        && name.find("/Bar_") != std::string::npos
        && name.find("/X_") != std::string::npos) {
        // This node is an X scintillator bar.
        TVector3 coord = LocalToMaster(TVector3(0,0,0));
        fP0DGeom->fActiveZ.push_back(coord.Z());
        fP0DGeom->fActiveXZ.push_back(coord.Z());
        COMET::IP0DGeom::Bar *bar
            = new COMET::IP0DGeom::Bar(bgeom, 0,
                                   coord.X(), coord.Z());
        fXBars.push_back(bar);
        fP0DGeom->fActiveMin.SetX(std::min(coord.X(),
                                           fP0DGeom->fActiveMin.X()));
        fP0DGeom->fActiveMax.SetX(std::max(coord.X(),
                                           fP0DGeom->fActiveMax.X()));
        fP0DGeom->fActiveMin.SetZ(std::min(coord.Z(),
                                           fP0DGeom->fActiveMin.Z()));
        fP0DGeom->fActiveMax.SetZ(std::max(coord.Z(),
                                           fP0DGeom->fActiveMax.Z()));
        if (name.find("/USECal_") != std::string::npos
            || name.find("/USECalSP0Dule_") != std::string::npos
            || name.find("/USPreCal_") != std::string::npos) {
            fP0DGeom->fUSECalMin.SetX(std::min(coord.X(),
                                               fP0DGeom->fUSECalMin.X()));
            fP0DGeom->fUSECalMax.SetX(std::max(coord.X(),
                                               fP0DGeom->fUSECalMax.X()));
            fP0DGeom->fUSECalMin.SetZ(std::min(coord.Z(),
                                               fP0DGeom->fUSECalMin.Z()));
            fP0DGeom->fUSECalMax.SetZ(std::max(coord.Z(),
                                               fP0DGeom->fUSECalMax.Z()));
        }
        if (name.find("/CECal_") != std::string::npos
            || name.find("/CECalSP0Dule_") != std::string::npos
            || name.find("/CPreCal_") != std::string::npos) {
            fP0DGeom->fCECalMin.SetX(std::min(coord.X(),
                                              fP0DGeom->fCECalMin.X()));
            fP0DGeom->fCECalMax.SetX(std::max(coord.X(),
                                              fP0DGeom->fCECalMax.X()));
            fP0DGeom->fCECalMin.SetZ(std::min(coord.Z(),
                                              fP0DGeom->fCECalMin.Z()));
            fP0DGeom->fCECalMax.SetZ(std::max(coord.Z(),
                                              fP0DGeom->fCECalMax.Z()));
        }
        if (name.find("/CTarget_") != std::string::npos
            || name.find("/USTarget_") != std::string::npos
            || name.find("/WaterSP0Dule_") != std::string::npos
            || name.find("/TargetSP0Dule_") != std::string::npos) {
            fP0DGeom->fWaterMin.SetX(std::min(coord.X(),
                                              fP0DGeom->fWaterMin.X()));
            fP0DGeom->fWaterMax.SetX(std::max(coord.X(),
                                              fP0DGeom->fWaterMax.X()));
            fP0DGeom->fWaterMin.SetZ(std::min(coord.Z(),
                                              fP0DGeom->fWaterMin.Z()));
            fP0DGeom->fWaterMax.SetZ(std::max(coord.Z(),
                                              fP0DGeom->fWaterMax.Z()));
        }
        if (name.find("/CarbonSP0Dule_") != std::string::npos) {
            // Deprecated from a previous design.
            fP0DGeom->fCarbonMin.SetX(std::min(coord.X(),
                                               fP0DGeom->fCarbonMin.X()));
            fP0DGeom->fCarbonMax.SetX(std::max(coord.X(),
                                               fP0DGeom->fCarbonMax.X()));
            fP0DGeom->fCarbonMin.SetZ(std::min(coord.Z(),
                                               fP0DGeom->fCarbonMin.Z()));
            fP0DGeom->fCarbonMax.SetZ(std::max(coord.Z(),
                                               fP0DGeom->fCarbonMax.Z()));
        }
        // Don't look any deeper into the hierarchy.
        return false;
    }
    else if (name.find("/P0Dule_") != std::string::npos
             && name.find("/Bar_") != std::string::npos
             && name.find("/Y_") != std::string::npos) {
        // This node is a Y scintillator bar.
        TVector3 coord = LocalToMaster(TVector3(0,0,0));
        fP0DGeom->fActiveZ.push_back(coord.Z());
        fP0DGeom->fActiveYZ.push_back(coord.Z());
        COMET::IP0DGeom::Bar *bar
            = new COMET::IP0DGeom::Bar(bgeom, 1,
                                   coord.Y(), coord.Z());
        fYBars.push_back(bar);
        fP0DGeom->fActiveMin.SetY(std::min(coord.Y(),
                                           fP0DGeom->fActiveMin.Y()));
        fP0DGeom->fActiveMax.SetY(std::max(coord.Y(),
                                           fP0DGeom->fActiveMax.Y()));
        fP0DGeom->fActiveMin.SetZ(std::min(coord.Z(),
                                           fP0DGeom->fActiveMin.Z()));
        fP0DGeom->fActiveMax.SetZ(std::max(coord.Z(),
                                           fP0DGeom->fActiveMax.Z()));
        if (name.find("/USECal_") != std::string::npos
            || name.find("/USECalSP0Dule") != std::string::npos
            || name.find("/USPreCal_") != std::string::npos) {
            fP0DGeom->fUSECalMin.SetY(std::min(coord.Y(),
                                               fP0DGeom->fUSECalMin.Y()));
            fP0DGeom->fUSECalMax.SetY(std::max(coord.Y(),
                                               fP0DGeom->fUSECalMax.Y()));
            fP0DGeom->fUSECalMin.SetZ(std::min(coord.Z(),
                                               fP0DGeom->fUSECalMin.Z()));
            fP0DGeom->fUSECalMax.SetZ(std::max(coord.Z(),
                                               fP0DGeom->fUSECalMax.Z()));
        }
        if (name.find("/CECal_") != std::string::npos
            || name.find("/CECalSP0Dule") != std::string::npos
            || name.find("/CPreCal_") != std::string::npos) {
            fP0DGeom->fCECalMin.SetY(std::min(coord.Y(),
                                              fP0DGeom->fCECalMin.Y()));
            fP0DGeom->fCECalMax.SetY(std::max(coord.Y(),
                                              fP0DGeom->fCECalMax.Y()));
            fP0DGeom->fCECalMin.SetZ(std::min(coord.Z(),
                                              fP0DGeom->fCECalMin.Z()));
            fP0DGeom->fCECalMax.SetZ(std::max(coord.Z(),
                                              fP0DGeom->fCECalMax.Z()));
        }
        if (name.find("/USTarget_") != std::string::npos
            || name.find("/CTarget_") != std::string::npos
            || name.find("/WaterSP0Dule_") != std::string::npos
            || name.find("/TargetSP0Dule_") != std::string::npos) {
            fP0DGeom->fWaterMin.SetY(std::min(coord.Y(),
                                              fP0DGeom->fWaterMin.Y()));
            fP0DGeom->fWaterMax.SetY(std::max(coord.Y(),
                                              fP0DGeom->fWaterMax.Y()));
            fP0DGeom->fWaterMin.SetZ(std::min(coord.Z(),
                                              fP0DGeom->fWaterMin.Z()));
            fP0DGeom->fWaterMax.SetZ(std::max(coord.Z(),
                                              fP0DGeom->fWaterMax.Z()));
        }
        if (name.find("/CarbonSP0Dule_") != std::string::npos) {
            // Deprecated from a previous design.
            fP0DGeom->fCarbonMin.SetY(std::min(coord.Y(),
                                               fP0DGeom->fCarbonMin.Y()));
            fP0DGeom->fCarbonMax.SetY(std::max(coord.Y(),
                                               fP0DGeom->fCarbonMax.Y()));
            fP0DGeom->fCarbonMin.SetZ(std::min(coord.Z(),
                                               fP0DGeom->fCarbonMin.Z()));
            fP0DGeom->fCarbonMax.SetZ(std::max(coord.Z(),
                                               fP0DGeom->fCarbonMax.Z()));
        }
        // Don't look any deeper into the hierarchy.
        return false;
    }
    else if ((name.find("/USTarget_") != std::string::npos
              || name.find("/CTarget_") != std::string::npos
              || name.find("/WaterSP0Dule_") != std::string::npos
              || name.find("/TargetSP0Dule_") != std::string::npos)
             && name.find("/Target_") != std::string::npos
             && name.find("/Water_") != std::string::npos) {
        // This node is water target.
        TVector3 coord = LocalToMaster(TVector3(0,0,0));
        fP0DGeom->fWaterZ.push_back(coord.Z());
    }

    return true;
}

void COMET::IP0DGeom::CompressPlanes(std::vector<double>& planes) {
    std::vector<double>::iterator r = planes.begin();
    for (std::vector<double>::iterator p = r; 
         p != planes.end();) {
        double sum = 0.0;
        double count = 0.0;
        std::vector<double>::iterator q = p;
        for(q = p;
            q != planes.end() && std::abs(*p - *q) < 1*unit::cm;
            ++q) {
            sum += *q;
            count += 1.0;
        }
        *(r++) = sum/count;
        p = q;
    }
    planes.erase(r,planes.end());
}

void COMET::IP0DGeom::Fill() {
    Clear();
    fBarMap.clear();
    fActiveMin = TVector3(unit::km, unit::km, unit::km);
    fActiveMax = TVector3(-unit::km, -unit::km, -unit::km);
    fUSECalMin = TVector3(unit::km, unit::km, unit::km);
    fUSECalMax = TVector3(-unit::km, -unit::km, -unit::km);
    fCECalMin = TVector3(unit::km, unit::km, unit::km);
    fCECalMax = TVector3(-unit::km, -unit::km, -unit::km);
    fWaterMin = TVector3(unit::km, unit::km, unit::km);
    fWaterMax = TVector3(-unit::km, -unit::km, -unit::km);
    fCarbonMin = TVector3(unit::km, unit::km, unit::km);
    fCarbonMax = TVector3(-unit::km, -unit::km, -unit::km);
    IP0DGeomVisitor visitor(this);
    visitor.VisitGeometry();
    fNeighborsXZ = new COMET::IDelaunay2D;
    fNeighborsYZ = new COMET::IDelaunay2D;
    // Shuffle the bar order since the Delaunay triangulation assumes "random"
    // points.  Adding the bars in order is a worst case scenerio for speed.
    // This cuts the time spent on these operations by about three.
    std::random_shuffle(visitor.fXBars.begin(),visitor.fXBars.end());
    for (std::vector<COMET::IP0DGeom::Bar*>::iterator b = visitor.fXBars.begin();
         b != visitor.fXBars.end();
         ++b) {
        fNeighborsXZ->AddPoint(*b);
    }
    fNeighborsXZ->Close();
    std::random_shuffle(visitor.fYBars.begin(),visitor.fYBars.end());
    for (std::vector<COMET::IP0DGeom::Bar*>::iterator b = visitor.fYBars.begin();
         b != visitor.fYBars.end();
         ++b) {
        fNeighborsYZ->AddPoint(*b);
    }
    fNeighborsYZ->Close();
    fXBars.clear();
    fYBars.clear();
    fBarMap.clear();
    for (COMET::TMeshPointSet::iterator p = fNeighborsXZ->BeginPoints();
         p != fNeighborsXZ->EndPoints();
         ++p) {
        COMET::IP0DGeom::Bar* bar = dynamic_cast<COMET::IP0DGeom::Bar*>(*p);
        for (COMET::IMeshPoint::EdgeIterator e = bar->BeginEdges();
             e != bar->EndEdges();
             ++e) {
            COMET::IP0DGeom::Bar* s 
                = dynamic_cast<COMET::IP0DGeom::Bar*>((*e)->GetOther(*p));
            bar->GetNeighbors().push_back(s->GetGeomId());
        }
        fBarMap[bar->GetGeomId()] = bar;
        fXBars.push_back(bar->GetGeomId());
    }
    for (COMET::TMeshPointSet::iterator p = fNeighborsYZ->BeginPoints();
         p != fNeighborsYZ->EndPoints();
         ++p) {
        COMET::IP0DGeom::Bar* bar = dynamic_cast<COMET::IP0DGeom::Bar*>(*p);
        for (COMET::IMeshPoint::EdgeIterator e = bar->BeginEdges();
             e != bar->EndEdges();
             ++e) {
            COMET::IP0DGeom::Bar* s 
                = dynamic_cast<COMET::IP0DGeom::Bar*>((*e)->GetOther(*p));
            bar->GetNeighbors().push_back(s->GetGeomId());
        }
        fBarMap[bar->GetGeomId()] = bar;
        fYBars.push_back(bar->GetGeomId());
    }
    sort(fActiveZ.begin(), fActiveZ.end());
    CompressPlanes(fActiveZ);
    sort(fActiveXZ.begin(), fActiveXZ.end());
    CompressPlanes(fActiveXZ);
    sort(fActiveYZ.begin(), fActiveYZ.end());
    CompressPlanes(fActiveYZ);
    sort(fWaterZ.begin(), fWaterZ.end());
    CompressPlanes(fWaterZ);
    // Find the pointing direction for each bar in the ROOT Geometry.
    gGeoManager->PushPath();
    for (std::map<COMET::IGeometryId,COMET::IP0DGeom::Bar*>::iterator b = fBarMap.begin();
         b != fBarMap.end();
         ++b) {
        //gGeoManager->CdNode(b->first);
	COMET::IOADatabase::Get().GeomId().CdId(b->first);
        TGeoVolume* volume = gGeoManager->GetCurrentVolume();
        TGeoTrap* shape = dynamic_cast<TGeoTrap*>(volume->GetShape()); 
        if (!shape) {
            std::cout << "*** Node " << b->first << "isn't a triangular bar"
                      << std::endl;
            continue;
        }
        COMET::IP0DGeom::Bar* bar = b->second;
        double local[3] = {0,0,0};
        double master[3] = {0,0,0};
        if (shape->GetTl1() < 0.5*shape->GetBl1()) {
            local[1] = 1.0; // point is on y axis.
            gGeoManager->LocalToMasterVect(local,master);
        }
        else if (shape->GetBl2() < 0.5*shape->GetBl1()) {
            local[2] = 1.0; // point is on z axis.
            gGeoManager->LocalToMasterVect(local,master);
        }
        else {
            // Possible square bar: Later software depends on triangles.
            std::cout << "*** Node " << b->first << "is invalid shape"
                      << std::endl;
        }
        bar->GetPointing().SetXYZ(master[0],master[1],master[2]);
    }
    gGeoManager->PopPath();

    
}

void COMET::IP0DGeom::Clear() {
    if (fNeighborsXZ) delete fNeighborsXZ;
    if (fNeighborsYZ) delete fNeighborsYZ;
    fNeighborsXZ = NULL;
    fNeighborsYZ = NULL;
    fXBars.clear();
    fYBars.clear();
    fActiveZ.clear();
    fActiveXZ.clear();
    fActiveYZ.clear();
    fWaterZ.clear();
    fActiveMin = TVector3(0,0,0);
    fActiveMax = TVector3(0,0,0);
    fUSECalMin = TVector3(0,0,0);
    fUSECalMax = TVector3(0,0,0);
    fCECalMin = TVector3(0,0,0);
    fCECalMax = TVector3(0,0,0);
    fWaterMin = TVector3(0,0,0);
    fWaterMax = TVector3(0,0,0);
    
    std::vector<COMET::IGeomModuleBase*>::iterator modules_it = fModules.begin();
	std::vector<COMET::IGeomModuleBase*>::iterator modules_end = fModules.end();
	for( ; modules_it != modules_end; ++modules_it)
	{
		if(*modules_it) delete (*modules_it);
	}
	fModules.clear();
}

const COMET::IP0DGeom::Bar& COMET::IP0DGeom::GetBar(COMET::IGeometryId geomId) const {
    std::map<COMET::IGeometryId,COMET::IP0DGeom::Bar*>::const_iterator s 
        = fBarMap.find(geomId);
    if (s == fBarMap.end()) throw COMET::ENoSuchBar();
    return *(s->second);
};

const COMET::IP0DGeom::Bar& COMET::IP0DGeom::GetBar(const COMET::IHit& hit) const {
    return GetBar(hit.GetGeomId());
};

const COMET::IP0DGeom::BarMap& COMET::IP0DGeom::GetBars() const {
    return fBarMap;
}

int COMET::IP0DGeom::FindZPlane(std::vector<double> planes, 
                             double z, double toler) const {
    double minDelta = 10*unit::meter;
    int bestPlane = -1; 
    for(unsigned int i = 0; i < planes.size(); ++i) {
        double surfaceZ = planes[i];
        double deltaZ = z - surfaceZ;
        if (deltaZ<0) deltaZ = - deltaZ;
        if (deltaZ<minDelta) {
            minDelta = deltaZ;
            bestPlane = i;
        }
    }
    if (minDelta>toler) return -1;
    return bestPlane;
}

int COMET::IP0DGeom::ActivePlane(double z) const {
    return FindZPlane(fActiveZ,z,1*unit::cm);
}

double COMET::IP0DGeom::ActivePlaneZ(int index) const {
    if (index < 0) return -100*unit::meter;
    if ((unsigned) index >= fActiveZ.size()) return 100*unit::meter;
    return fActiveZ[index];
}

int COMET::IP0DGeom::ActivePlaneCount(void) const {
    return fActiveZ.size();
}

int COMET::IP0DGeom::ActiveXPlane(double z) const {
    return FindZPlane(fActiveXZ,z,1*unit::cm);
}

double COMET::IP0DGeom::ActiveXPlaneZ(int index) const {
    if (index < 0) return -100*unit::meter;
    if ((unsigned) index >= fActiveXZ.size()) return 100*unit::meter;
    return fActiveXZ[index];
}

int COMET::IP0DGeom::ActiveXPlaneCount(void) const {
    return fActiveXZ.size();
}

int COMET::IP0DGeom::ActiveYPlane(double z) const {
    return FindZPlane(fActiveYZ,z,1*unit::cm);
}

double COMET::IP0DGeom::ActiveYPlaneZ(int index) const {
    if (index < 0) return -100*unit::meter;
    if ((unsigned) index >= fActiveYZ.size()) return 100*unit::meter;
    return fActiveYZ[index];
}

int COMET::IP0DGeom::ActiveYPlaneCount(void) const {
    return fActiveYZ.size();
}

double COMET::IP0DGeom::WaterDistance(double checkZ) const {
    double minDist = 100*unit::meter;
    for (std::vector<double>::const_iterator z = fWaterZ.begin();
         z != fWaterZ.end();
         ++z) {
        double deltaZ = (*z) - checkZ;
        if (std::abs(deltaZ) < std::abs(minDist)) minDist = deltaZ;
    }
    return minDist;
}

COMET::IP0DGeom::Bar::Bar(COMET::IGeometryId geomId, 
                       int orientation, double x, double y)
    : COMET::IScintBarGeom(geomId,orientation,1), fX(x), fY(y) {}

COMET::IP0DGeom::Bar::~Bar() {}

int COMET::IP0DGeom::Bar::GetSuperP0Dule() const {
    int p0dule = GetP0Dule();
    if (p0dule<7) return 0;     // Upstream ECal.
    if (p0dule<20) return 1;    // Upstream Water Target.
    if (p0dule<33) return 2;    // Upstream Water Target.
    return 3;
}

int COMET::IP0DGeom::Bar::GetP0Dule() const {
    return GetGeomId().GetField(COMET::GeomId::Def::P0D::Bar::kP0DuleMSB,
                                COMET::GeomId::Def::P0D::Bar::kP0DuleLSB);
}

int COMET::IP0DGeom::Bar::GetLayer() const {
    return GetGeomId().GetField(COMET::GeomId::Def::P0D::Bar::kLayerMSB,
                                COMET::GeomId::Def::P0D::Bar::kLayerLSB);
}

int COMET::IP0DGeom::Bar::GetNumber() const {
    return GetGeomId().GetField(COMET::GeomId::Def::P0D::Bar::kBarMSB,
                                COMET::GeomId::Def::P0D::Bar::kBarLSB);
}

double COMET::IP0DGeom::EdgeDistance(double Z) const {
    double zDist = std::min(Z - fActiveMin.Z(), fActiveMax.Z() - Z);
    return zDist;
}

double COMET::IP0DGeom::EdgeDistance(double X, double Y, double Z) const {
    double xDist = std::min(X - fActiveMin.X(), fActiveMax.X() - X);
    double yDist = std::min(Y - fActiveMin.Y(), fActiveMax.Y() - Y);
    double zDist = 100*unit::m;
    if (Z > -99*unit::m) zDist = EdgeDistance(Z);
    return std::min(xDist,std::min(yDist,zDist));
}

double COMET::IP0DGeom::WaterFiducialDistance(double Z) const {
    double zDist = std::min(Z - fWaterMin.Z(), fWaterMax.Z() - Z);
    return zDist;
}

double COMET::IP0DGeom::WaterFiducialDistance(double X, double Y, double Z,
                                           double xyFiducial,
                                           double zFiducial) const {
    double xDist = std::min(X - fWaterMin.X(), fWaterMax.X() - X)-xyFiducial;
    double yDist = std::min(Y - fWaterMin.Y(), fWaterMax.Y() - Y)-xyFiducial;
    double zDist = 100*unit::m;
    if (Z > -99*unit::m) zDist = WaterFiducialDistance(Z) - zFiducial;
    return std::min(xDist,std::min(yDist,zDist));
}

double COMET::IP0DGeom::TargetDepth() const {
    double depth = 0;
    for (int t = 0; t < 25; ++t) {
        depth += TargetDepth(COMET::GeomId::P0D::Target(t));
    }
    return depth/25.0;
}

double COMET::IP0DGeom::TargetDepth(COMET::IGeometryId target, int side) const {
    double depth = -99*unit::m;
    gGeoManager->PushPath();
    if (COMET::IOADatabase::Get().GeomId().CdId(target)) {
        TGeoVolume *targetVolume = gGeoManager->GetCurrentVolume();
        depth = 0.0;
        if (side >= 0) {
            TGeoNode* node = targetVolume->GetNode("LeftTarget_0");
            TGeoVolume* volume = node->GetVolume();
            TGeoBBox* box = dynamic_cast<TGeoBBox*>(volume->GetShape());
            depth += 2*box->GetDY();
        }
        if (side <= 0) {
            TGeoNode* node = targetVolume->GetNode("RightTarget_0");
            TGeoVolume* volume = node->GetVolume();
            TGeoBBox* box = dynamic_cast<TGeoBBox*>(volume->GetShape());
            depth += 2*box->GetDY();
        }
        if (side == 0) {
            depth /= 2;
        }
    }
    gGeoManager->PopPath();
    return depth;
}

void COMET::IP0DGeom::ls(const char*) {

    COMETLog( "P0D minimum active edges : "
              << fActiveMin.X() 
              << " " << fActiveMin.Y() 
              << " " << fActiveMin.Z());

    COMETLog( "P0D maximum active edges : "
              << fActiveMax.X() 
              << " " << fActiveMax.Y() 
              << " " << fActiveMax.Z());

    COMETLog( "P0D active volume center : "
              << (fActiveMax.X()+fActiveMin.X())/2.0
              << " " << (fActiveMax.Y()+fActiveMin.Y())/2.0
              << " " << (fActiveMax.Z()+fActiveMin.Z())/2.0);

}