#include <iostream>
#include <vector>
#include <string>
#include <cstring>
#include <cmath>

#include <TGeoManager.h>
#include <TGeoNode.h>
#include <TGeoVolume.h>
#include <TGeoBBox.h>

#include <IOADatabase.hxx>
#include <IHit.hxx>
#include <HEPUnits.hxx>
#include <HEPConstants.hxx>
#include <oaEvent_version.h>

#include "IGeomInfo.hxx"

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

// Declare the callback, but make sure it's in the private namespace.
void COMET::IGeomInfo::IGeometryChange::Callback(const COMET::ICOMETEvent* const) {
    if (!COMET::IGeomInfo::fGeomInfo) return;
    delete COMET::IGeomInfo::fGeomInfo;
    COMET::IGeomInfo::fGeomInfo = new COMET::IGeomInfo;
}
namespace {COMET::IGeomInfo::IGeometryChange geomInfoGeometryChange;};

// The static member pointer to the singleton.
COMET::IGeomInfo* COMET::IGeomInfo::fGeomInfo = NULL;

COMET::IGeomInfo::IGeomInfo()
    : fGeoManager(NULL),
      fP0DGeom(NULL),
      /*
      fFGDGeom(NULL),
      fECALGeom(NULL),
      */
      fSMRDGeom(NULL),
      //fTPCGeom(NULL),
      fIngridGeom(NULL),
      fScintDetGeom(NULL),
      fCDCGeom(NULL) {
    FillInformation();
    COMET::IOADatabase::Get().RegisterGeometryCallback(&geomInfoGeometryChange);
}

COMET::IGeomInfo::~IGeomInfo() {
    if (fP0DGeom) delete fP0DGeom;
    /*
    if (fFGDGeom) delete fFGDGeom;
    if (fECALGeom) delete fECALGeom;
    */
    if (fSMRDGeom) delete fSMRDGeom;
    //if (fTPCGeom) delete fTPCGeom;
    if (fIngridGeom) delete fIngridGeom;
    if (fScintDetGeom) delete fScintDetGeom;
    if (fCDCGeom) delete fCDCGeom;
    COMET::IOADatabase::Get().RemoveGeometryCallback(&geomInfoGeometryChange);
    fGeomInfo = NULL;
}

COMET::IGeomInfo& COMET::IGeomInfo::Get(void){
    if (!fGeomInfo 
        || fGeomInfo->fGeoManager != COMET::IOADatabase::Get().Geometry()) {
        if (fGeomInfo) delete fGeomInfo;
        COMETInfo("Create a new Geometry Info object.");
        fGeomInfo = new COMET::IGeomInfo;
    }
    return *fGeomInfo;
}

COMET::IP0DGeom& COMET::IGeomInfo::P0D(void){
    COMET::IGeomInfo::Get();
    if (!fGeomInfo->fP0DGeom) {
        fGeomInfo->fP0DGeom = new COMET::IP0DGeom();
        fGeomInfo->fP0DGeom->Fill();
    }
    return *(fGeomInfo->fP0DGeom);
}
/*
COMET::IFGDGeom& COMET::IGeomInfo::FGD(void){
    COMET::IGeomInfo::Get();
    if (!fGeomInfo->fFGDGeom) {
        fGeomInfo->fFGDGeom = new COMET::IFGDGeom();
        fGeomInfo->fFGDGeom->Fill();
    }
    return *(fGeomInfo->fFGDGeom);
}

COMET::IECALGeom& COMET::IGeomInfo::ECAL(void){
    COMET::IGeomInfo::Get();
    if (!fGeomInfo->fECALGeom) {
        fGeomInfo->fECALGeom = new COMET::IECALGeom();
        fGeomInfo->fECALGeom->Fill();
    }
    return *(fGeomInfo->fECALGeom);
}
*/
COMET::ISMRDGeom& COMET::IGeomInfo::SMRD(void){
    COMET::IGeomInfo::Get();
    if (!fGeomInfo->fSMRDGeom) {
        fGeomInfo->fSMRDGeom = new COMET::ISMRDGeom("/Magnet","/MRD","/Bar");
        fGeomInfo->fSMRDGeom->Fill();
    }
    return *(fGeomInfo->fSMRDGeom);
}


COMET::IScintDetGeom& COMET::IGeomInfo::ScintDet(void){
    COMET::IGeomInfo::Get();
    if (!fGeomInfo->fScintDetGeom) {
        fGeomInfo->fScintDetGeom = new COMET::IScintDetGeom("/Magnet","/MRD","/Bar");
        fGeomInfo->fScintDetGeom->Fill();
    }
    return *(fGeomInfo->fScintDetGeom);
}

/*
COMET::ITPCGeom& COMET::IGeomInfo::TPC(void){
    COMET::IGeomInfo::Get();
    if (!fGeomInfo->fTPCGeom) {
        fGeomInfo->fTPCGeom = new COMET::ITPCGeom();
        fGeomInfo->fTPCGeom->Fill();
    }
    return *(fGeomInfo->fTPCGeom);
}
*/

COMET::IIngridGeom& COMET::IGeomInfo::Ingrid(void){
    COMET::IGeomInfo::Get();
    if (!fGeomInfo->fIngridGeom) {
        fGeomInfo->fIngridGeom = new COMET::IIngridGeom();
        fGeomInfo->fIngridGeom->Fill();
    }
    return *(fGeomInfo->fIngridGeom);
}


COMET::ICDCGeom& COMET::IGeomInfo::CDC(void){
    COMET::IGeomInfo::Get();
    if (!fGeomInfo->fCDCGeom) {
        fGeomInfo->fCDCGeom = new COMET::ICDCGeom();
        fGeomInfo->fCDCGeom->Fill();
    }
    return *(fGeomInfo->fCDCGeom);
}


void COMET::IGeomInfo::Clear() {
    if (fP0DGeom) fP0DGeom->Clear();
    /*
    if (fFGDGeom) fFGDGeom->Clear();
    if (fECALGeom) fECALGeom->Clear();
    */
    if (fSMRDGeom) fSMRDGeom->Clear();
    //if (fTPCGeom)  fTPCGeom->Clear();
    if (fIngridGeom) fIngridGeom->Clear();
    if (fScintDetGeom) fScintDetGeom->Clear();
    if (fCDCGeom) fCDCGeom->Clear();
}

void COMET::IGeomInfo::FillInformation(void) {
    fGeoManager = COMET::IOADatabase::Get().Geometry();
    Clear();
    if (fP0DGeom) fP0DGeom->Fill();
    /*
    if (fFGDGeom) fFGDGeom->Fill();
    if (fECALGeom) fECALGeom->Fill();
    */
    if (fSMRDGeom) fSMRDGeom->Fill();
    //if (fTPCGeom)  fTPCGeom->Fill();
    if (fIngridGeom) fIngridGeom->Fill();
    if (fScintDetGeom) fScintDetGeom->Fill();
    if (fCDCGeom) fCDCGeom->Fill();
}

COMET::IGeometryId COMET::IGeomInfo::GeomId(double x, double y, double z) {
    COMET::IGeometryId geomId;
    COMET::IOADatabase::Get().GeomId().GetGeometryId(x, y, z, geomId);
    return geomId;
}

TVector3 COMET::IGeomInfo::NodeExtent(const COMET::IHit& hit) {
    return NodeExtent(hit.GetGeomId());
}

TVector3 COMET::IGeomInfo::NodeExtent(COMET::IGeometryId geomId) {
    gGeoManager->PushPath();
    COMET::IOADatabase::Get().GeomId().CdId(geomId);
    TGeoNode* node = gGeoManager->GetCurrentNode();
    TGeoBBox *shape
        =  dynamic_cast<TGeoBBox*>(node->GetVolume()->GetShape());
    double spread[3] = {0,0,0};
    const double localX[3] = {shape->GetDX(),0,0};
    const double localY[3] = {0,shape->GetDY(),0};
    const double localZ[3] = {0,0,shape->GetDZ()};
    double master[3] = {0,0,0};
    gGeoManager->LocalToMasterVect(localX,master);
    for (int i=0; i<3; ++i)
        spread[i] = std::max(spread[i],std::abs(master[i]));
    gGeoManager->LocalToMasterVect(localY,master);
    for (int i=0; i<3; ++i)
        spread[i] = std::max(spread[i],std::abs(master[i]));
    gGeoManager->LocalToMasterVect(localZ,master);
    for (int i=0; i<3; ++i)
        spread[i] = std::max(spread[i],std::abs(master[i]));
    gGeoManager->PopPath();
    return TVector3(spread);
}

TVector3 COMET::IGeomInfo::NodeSize(COMET::IGeometryId geomId) {
    gGeoManager->PushPath();
    COMET::IOADatabase::Get().GeomId().CdId(geomId);
    TGeoNode* node = gGeoManager->GetCurrentNode();
    TGeoBBox *shape =  dynamic_cast<TGeoBBox*>(node->GetVolume()->GetShape());
    gGeoManager->PopPath();
    return TVector3(shape->GetDX(), shape->GetDY(), shape->GetDZ());
}

std::string COMET::IGeomInfo::NodeName(COMET::IGeometryId geomId) { 
    return geomId.GetName(); 
}

TVector3 COMET::IGeomInfo::NodePosition(COMET::IGeometryId geomId) {
    return geomId.GetPosition();
}

TVector3 COMET::IGeomInfo::MasterToLocal(COMET::IGeometryId geomId,
                                      double x, double y, double z) {
    gGeoManager->PushPath();
    COMET::IOADatabase::Get().GeomId().CdId(geomId);
    const double master[3] = {x,y,z};
    double local[3];
    gGeoManager->MasterToLocal(master,local);
    gGeoManager->PopPath();
    return TVector3(local);
}

TVector3 COMET::IGeomInfo::MasterToLocalVect(COMET::IGeometryId geomId,
                                          double x, double y, double z) {
    gGeoManager->PushPath();
    COMET::IOADatabase::Get().GeomId().CdId(geomId);
    const double master[3] = {x,y,z};
    double local[3];
    gGeoManager->MasterToLocalVect(master,local);
    gGeoManager->PopPath();
    return TVector3(local);
}

TVector3 COMET::IGeomInfo::LocalToMaster(COMET::IGeometryId geomId,
                                      double x, double y, double z) {
    gGeoManager->PushPath();
    COMET::IOADatabase::Get().GeomId().CdId(geomId);
    const double local[3] = {x,y,z};
    double master[3];
    gGeoManager->LocalToMaster(local,master);
    gGeoManager->PopPath();
    return TVector3(master);
}

TVector3 COMET::IGeomInfo::LocalToMasterVect(COMET::IGeometryId geomId,
                                          double x, double y, double z) {
    gGeoManager->PushPath();
    COMET::IOADatabase::Get().GeomId().CdId(geomId);
    const double local[3] = {x,y,z};
    double master[3];
    gGeoManager->LocalToMasterVect(local,master);
    gGeoManager->PopPath();
    return TVector3(master);
}

const COMET::IScintBarGeom& 
COMET::IGeomInfo::GetScintBar(COMET::IGeometryId geomId) const {
  TString pathname(geomId.GetName());
  if (pathname.Contains("/P0D_") &&
      pathname.Contains("/P0Dule") &&
      pathname.Contains("/Bar")){
      return COMET::IGeomInfo::P0D().GetBar(geomId);
  }
  /*
  else if (pathname.Contains("/Tracker") &&
      pathname.Contains("/Scint") &&
      pathname.Contains("/Bar")){
      return COMET::IGeomInfo::FGD().GetBar(geomId);
  }
  else if (pathname.Contains("ECal_") &&
	   pathname.Contains("/Scint") &&
	   pathname.Contains("/Bar")){
      return COMET::IGeomInfo::ECAL().GetBar(geomId);
  }
  */
  else if (pathname.Contains("Clam_") &&
	   pathname.Contains("/MRD") &&
	   pathname.Contains("/Bar_")){
      return COMET::IGeomInfo::SMRD().GetBar(geomId);
  }
  else if (pathname.Contains("/ingrid") &&
            ( ( ( ( pathname.Contains("/standard") && pathname.Contains("/tracker")  ) ||
                  ( pathname.Contains("/proton")   && pathname.Contains("/ptracker") )   ) &&
                ( pathname.Contains("/scinti_v") || pathname.Contains("/scinti_h") ) ) ||
              ( pathname.Contains("/veto_h") && pathname.Contains("/vscinti_h") ) ||
              ( pathname.Contains("/veto_v") && pathname.Contains("/vscinti_v") )
            )
          ){
      return COMET::IGeomInfo::Ingrid().GetBar(geomId);
  }
  else {
      COMETSevere("Volume not found " << pathname << " (" << geomId.AsInt() << ")");
      throw COMET::EVolumeNotFound();
  }

  static const COMET::IScintBarGeom missing(IGeometryId(0),0,0);
  return missing;
}

double COMET::IGeomInfo::X0(TGeoMaterial * mat, int method)
{
  double x0 = 0;
  switch ( method ) {
    case 2: // simplified PDG method (PDG 2006 sec 27.4)
    {
      double Z = mat->GetZ();
      if ( Z > 0. ) {
        x0 = 716.4*unit::g/unit::cm2*mat->GetA()/(unit::g/unit::mole)/
             (Z*(Z+1)*log(287./sqrt(Z)));
      }
      x0 /= mat->GetDensity();
    }
    break;
    case 1: // full PDG method (PDG 2006 sec 27.4)
    {
      double Lrad       = 0.;
      double Lrad_prime = 0.;
      double Z          = mat->GetZ();
      if      ( Z == 1. ) {
        Lrad       = 5.31;
        Lrad_prime = 6.144;
      }
      else if ( Z == 2. ) {
        Lrad       = 4.79;
        Lrad_prime = 5.621;
      }
      else if ( Z == 3. ) {
        Lrad       = 4.74;
        Lrad_prime = 5.805;
      }
      else if ( Z == 4. ) {
        Lrad       = 4.71;
        Lrad_prime = 5.924;
      }
      else {
        Lrad       = log(184.15*pow(Z,-1./3.));
        Lrad_prime = log(1194. *pow(Z,-2./3.));
      }
      double a2    = pow(unit::fine_structure_const*Z,2.);
      double f_Z   = a2*(1./(1.+a2)+0.20206-0.0369*a2+0.0083*a2*a2+0.0082*a2*a2*a2);
      double x0inv = (1./(716.408*unit::g/unit::cm2))/
      (mat->GetA()/(unit::g/unit::mole))*
      (Z*Z*(Lrad-f_Z)+Z*Lrad_prime);
      if ( x0inv != 0. ) x0 = 1./x0inv;
      x0 /= mat->GetDensity();
    }
    break;
    default:
    case 0:
    {
      x0 = mat->GetRadLen()/1000.*unit::cm;
    }
    break;
  }
  return x0;
}

double COMET::IGeomInfo::X0(TGeoMaterial * mat)
{
  double x0 = 0;
  double Z = mat->GetZ();
  if ( Z > 0. ) {
    x0 = 716.4*unit::g/unit::cm2*mat->GetA()/(unit::g/unit::mole)/
      (Z*(Z+1)*log(287./sqrt(Z)));
  }
  x0 /= mat->GetDensity();
  return x0;
}