#ifndef GEOMETRY_ECALGEOMETRYLOADER_ICC
#define GEOMETRY_ECALGEOMETRYLOADER_ICC 1

#include "Geometry/EcalTestBeam/test/ee/CaloGeometryLoaderTest.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "DetectorDescription/Core/interface/DDCompactView.h"
#include "DetectorDescription/Core/interface/DDSolid.h"
#include "DetectorDescription/Core/interface/DDSpecifics.h"

#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Geometry/Plane3D.h>

#include <vector>
#include <fstream>

typedef CaloCellGeometry::CCGFloat CCGFloat;

template <class T>
const double CaloGeometryLoaderTest<T>::k_ScaleFromDDDtoGeant(0.1);

template <class T>
CaloGeometryLoaderTest<T>::CaloGeometryLoaderTest() {
  m_filter.setCriteria(DDValue("Volume",
                               "EndcapSC",
                               //				  "EESCEnv1",
                               0),
                       DDCompOp::equals);
}

template <class T>
typename CaloGeometryLoaderTest<T>::PtrType CaloGeometryLoaderTest<T>::load(const DDCompactView* cpv,
                                                                            const Alignments* alignments,
                                                                            const Alignments* globals) {
  PtrType geom = std::make_unique<T>();

  makeGeometry(cpv, dynamic_cast<T*>(geom.get()), alignments, globals);

  return geom;
}

template <class T>
void CaloGeometryLoaderTest<T>::makeGeometry(const DDCompactView* cpv,
                                             T* geom,
                                             const Alignments* alignments,
                                             const Alignments* globals) {
  geom->allocateCorners(T::k_NumberOfCellsForCorners);
  geom->allocatePar(T::k_NumberOfParametersPerShape * T::k_NumberOfShapes, T::k_NumberOfParametersPerShape);

  DDFilteredView fv(*cpv, m_filter);

  for (bool doSubDets = fv.firstChild(); doSubDets; doSubDets = fv.nextSibling()) {
    const DDSolid& solid(fv.logicalPart().solid());

    const ParmVec& parameters(solid.parameters());

    DD3Vector x, y, z;
    fv.rotation().GetComponents(x, y, z);
    const CLHEP::HepRep3x3 temp(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
    const CLHEP::HepRotation hr(temp);
    const CLHEP::Hep3Vector h3v(fv.translation().X(), fv.translation().Y(), fv.translation().Z());
    const HepGeom::Transform3D ht3d(hr,  // only scale translation
                                    k_ScaleFromDDDtoGeant * h3v);

    const unsigned int id(getDetIdForDDDNode(fv));

    myFillGeom(geom, parameters, ht3d, id);
  }
}

template <class T>
void CaloGeometryLoaderTest<T>::myFillGeom(T* geom,
                                           const ParmVec& vv,
                                           const HepGeom::Transform3D& tr,
                                           const unsigned int id) {
  static const double kRadToDeg(180 / M_PI);
  static const HepGeom::Point3D<double> PIVOT(13533.35, -000.02, 0);        //-152.23 ) ;
  static const HepGeom::Point3D<double> FLOOR(13533.35, -000.02, -152.23);  //-152.23 ) ;

  static const HepGeom::Point3D<double> IP(0, 0, 0);  //-152.23 ) ;
  HepGeom::Point3D<double> AA0;
  HepGeom::Point3D<double> BB0;
  HepGeom::Point3D<double> CC0;
  HepGeom::Point3D<double> DD0;
  HepGeom::Point3D<double> EE0;
  HepGeom::Point3D<double> FF0;
  HepGeom::Point3D<double> GG0;
  HepGeom::Point3D<double> HH0;

  std::vector<CCGFloat> pv;
  pv.reserve(vv.size());
  for (unsigned int i(0); i != vv.size(); ++i) {
    const double factor(1 == i || 2 == i || 6 == i || 10 == i ? 1 : k_ScaleFromDDDtoGeant);
    pv.push_back(factor * vv[i]);
  }

  std::vector<GlobalPoint> corners(8);

  TruncatedPyramid::createCorners(pv, tr, corners);

  if ((160 == id || 164 == id || 134 == id || 133 == id || 153 == id || 150 == id || 130 == id) &&
      (0 < corners[0].x() && 0 < corners[0].y() && 0 < corners[0].z())) {
    /*      std::cout<<"\n *************** id = "<<id
	       <<"\n "<<corners[0]
	       <<"\n "<<corners[1]
	       <<"\n "<<corners[2]
	       <<"\n "<<corners[3]
	       <<"\n "<<corners[4]
	       <<"\n "<<corners[5]
	       <<"\n "<<corners[6]
	       <<"\n "<<corners[7]
	       <<std::endl ;*/

    if (160 == id)
      AA0 = HepGeom::Point3D<double>(corners[3].x(), corners[3].y(), corners[3].z());
    if (164 == id)
      BB0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
    if (134 == id)
      CC0 = HepGeom::Point3D<double>(corners[1].x(), corners[1].y(), corners[1].z());
    if (130 == id)
      DD0 = HepGeom::Point3D<double>(corners[0].x(), corners[0].y(), corners[0].z());
    if (150 == id)
      EE0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
    if (153 == id)
      FF0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
    if (133 == id)
      GG0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
    if (130 == id)
      HH0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
  }
  static bool done(false);
  if (0.001 < AA0.mag() && 0.001 < BB0.mag() && 0.001 < CC0.mag() && 0.001 < DD0.mag() && !done) {
    done = true;

    const HepGeom::Transform3D H4ToBeamLine(PIVOT + HepGeom::Point3D<double>(0, 1, 0),
                                            PIVOT + HepGeom::Point3D<double>(0, 0, 1),
                                            PIVOT + HepGeom::Point3D<double>(1, 0, 0),
                                            HepGeom::Point3D<double>(-1, 0, 0),
                                            HepGeom::Point3D<double>(0, -1, 0),
                                            HepGeom::Point3D<double>(0, 0, 1));
    /*					  HepGeom::Point3D<double> ( 0, -1, 0),
					  HepGeom::Point3D<double> ( 1,  0, 0),
					  HepGeom::Point3D<double> ( 0,  0, 1) ) ;*/

    std::cout << "\nH4ToBeamLine"
              << "   Theta_x = " << H4ToBeamLine.getRotation().thetaX() * kRadToDeg
              << ",  Theta_y = " << H4ToBeamLine.getRotation().thetaY() * kRadToDeg
              << ",  Theta_z = " << H4ToBeamLine.getRotation().thetaZ() * kRadToDeg
              << ",  eta = " << -log(tan(H4ToBeamLine.getRotation().thetaZ() / 2.)) << std::endl;

    std::cout << "\nH4ToBeamLine"
              << "   Phi_x = " << H4ToBeamLine.getRotation().phiX() * kRadToDeg
              << ",  Phi_y = " << H4ToBeamLine.getRotation().phiY() * kRadToDeg
              << ",  Phi_z = " << H4ToBeamLine.getRotation().phiZ() * kRadToDeg << std::endl;

    std::cout << "pivot to beamline=" << H4ToBeamLine * PIVOT << std::endl;
    std::cout << "endcap to beamline=" << H4ToBeamLine * HepGeom::Point3D<double>(PIVOT.x() + 316, 0, 50) << std::endl;

    std::cout << "\nA0 = " << AA0 << "\nB0 = " << BB0 << "\nC0 = " << CC0 << "\nD0 = " << DD0 << std::endl;

    const HepGeom::Point3D<double>& A0(AA0);
    const HepGeom::Point3D<double>& B0(BB0);
    const HepGeom::Point3D<double>& C0(CC0);
    const HepGeom::Point3D<double>& D0(DD0);

    const double angle((C0 - B0).angle(A0 - B0));

    static const std::vector<HepGeom::Point3D<double> > extra = {HepGeom::Point3D<double>(13859.05, 005.79, -002.66),
                                                                 HepGeom::Point3D<double>(13860.98, -038.39, -002.11),
                                                                 HepGeom::Point3D<double>(13853.19, -039.25, -030.77),
                                                                 HepGeom::Point3D<double>(13851.23, 004.99, -031.39)};

    static const std::vector<HepGeom::Point3D<double> > sur = {
        HepGeom::Point3D<double>(13862.30, 20.68, 11.60),    HepGeom::Point3D<double>(13865.69, -52.85, 12.64),
        HepGeom::Point3D<double>(13851.53, -54.27, -44.75),  HepGeom::Point3D<double>(13848.21, 19.18, -45.89),

        HepGeom::Point3D<double>(13861.72, 28.62, 11.61),    HepGeom::Point3D<double>(13866.88, -44.81, 12.64),
        HepGeom::Point3D<double>(13852.76, -46.56, -44.73),  HepGeom::Point3D<double>(13847.67, 26.79, -45.88),
        HepGeom::Point3D<double>(13856.82, 63.41, 11.60),    HepGeom::Point3D<double>(13869.76, -9.06, 12.64),
        HepGeom::Point3D<double>(13855.90, -12.31, -44.73),  HepGeom::Point3D<double>(13843.04, 60.08, -45.89),
        HepGeom::Point3D<double>(13860.86, -36.96, 11.60),   HepGeom::Point3D<double>(13851.38, -109.96, 12.63),
        HepGeom::Point3D<double>(13837.21, -108.89, -44.73), HepGeom::Point3D<double>(13846.73, -35.98, -45.88),
        HepGeom::Point3D<double>(13854.03, -76.21, 11.60),   HepGeom::Point3D<double>(13835.82, -147.53, 12.64),
        HepGeom::Point3D<double>(13821.88, -144.76, -44.73), HepGeom::Point3D<double>(13840.13, -73.54, -45.88),
        HepGeom::Point3D<double>(13861.10, 35.48, 09.88),    HepGeom::Point3D<double>(13867.78, -37.82, 11.07),
        HepGeom::Point3D<double>(13853.42, -40.01, -46.22),  HepGeom::Point3D<double>(13846.80, 33.21, -47.52),
        HepGeom::Point3D<double>(13860.68, 35.23, 19.89),    HepGeom::Point3D<double>(13867.35, -38.08, 20.23),
        HepGeom::Point3D<double>(13854.67, -39.43, -37.48),  HepGeom::Point3D<double>(13848.07, 33.79, -37.92),
        HepGeom::Point3D<double>(13859.95, 34.97, 29.87),    HepGeom::Point3D<double>(13866.62, -38.34, 29.33),
        HepGeom::Point3D<double>(13855.61, -38.88, -28.76),  HepGeom::Point3D<double>(13849.00, 34.34, -28.32),
        HepGeom::Point3D<double>(13858.93, 34.41, 39.81),    HepGeom::Point3D<double>(13865.65, -38.90, 38.40),
        HepGeom::Point3D<double>(13856.34, -38.58, -19.97),  HepGeom::Point3D<double>(13849.70, 34.65, -18.68),
        HepGeom::Point3D<double>(13857.67, 33.69, 49.76),    HepGeom::Point3D<double>(13864.43, -39.60, 47.48),
        HepGeom::Point3D<double>(13856.81, -38.43, -11.13),  HepGeom::Point3D<double>(13850.14, 34.78, -8.95)};

    double fsumx(0);
    double fsumy(0);
    double fsumz(0);
    double fsumx2(0);
    double fsumy2(0);
    double fsumz2(0);

    double psumx(0);
    double psumy(0);
    double psumz(0);
    double psumx2(0);
    double psumy2(0);
    double psumz2(0);

    double sumx(0);
    double sumy(0);
    double sumz(0);
    double sumx2(0);
    double sumy2(0);
    double sumz2(0);
    unsigned int nsur(sur.size() / 4);

    const std::string names[10] = {"H0", "H1", "H2", "H3", "H4", "V1", "V2", "V3", "V4", "V5"};

    HepGeom::Transform3D V3Transform;

    for (unsigned int i(0); i != nsur; ++i) {
      const HepGeom::Point3D<double> SA(sur[4 * i]);
      const HepGeom::Point3D<double> SB(sur[4 * i + 1]);
      const HepGeom::Point3D<double> SC(sur[4 * i + 2]);
      const HepGeom::Point3D<double> SD(sur[4 * i + 3]);
      /*	 const HepGeom::Point3D<double>  SD ( sur[4*i] ) ;
	 const HepGeom::Point3D<double>  SA ( sur[4*i+1] ) ;
	 const HepGeom::Point3D<double>  SB ( sur[4*i+2] ) ;
	 const HepGeom::Point3D<double>  SC ( sur[4*i+3] ) ;*/

      const HepGeom::Point3D<double>& aa(SA);
      const HepGeom::Point3D<double> bb(SA + (SB - SA).unit() * (B0 - A0).mag());

      const HepGeom::Plane3D<double> plane(aa, bb, SC);

      const HepGeom::Point3D<double> cc(
          bb +
          (HepGeom::Rotate3D(-angle, plane.normal()) * HepGeom::Vector3D<double>(aa - bb)).unit() * (C0 - B0).mag());

      const HepGeom::Transform3D trform(A0, B0, C0, aa, bb, cc);

      const HepGeom::Transform3D trforminv(aa, bb, cc, A0, B0, C0);

      const HepGeom::Point3D<double> dd(trform * D0);

      const HepGeom::Transform3D BeamLineToCMS(H4ToBeamLine * aa, H4ToBeamLine * bb, H4ToBeamLine * cc, A0, B0, C0);

      if (7 == i)
        V3Transform = BeamLineToCMS;

      const CLHEP::HepRotation blcmsrot(BeamLineToCMS.getRotation());
      const HepGeom::Vector3D<double> blcmstra(BeamLineToCMS.getTranslation());

      std::cout << "\n********* Name is " << names[i];

      std::cout << "\n\n B discrepancy is " << (bb - SB) << ", " << (bb - SB).mag() << "\n\n C discrepancy is "
                << (cc - SC) << ", " << (cc - SC).mag() << "\n\n D discrepancy is " << (dd - SD) << ", "
                << (dd - SD).mag() << std::endl;
      /*
	 const double dxbb ( bb.x() - SB.x() ) ;
	 const double dybb ( bb.y() - SB.y() ) ;
	 const double dzbb ( bb.z() - SB.z() ) ;

	 const double dxcc ( cc.x() - SC.x() ) ;
	 const double dycc ( cc.y() - SC.y() ) ;
	 const double dzcc ( cc.z() - SC.z() ) ;
*/
      const double dxdd(dd.x() - SD.x());
      const double dydd(dd.y() - SD.y());
      const double dzdd(dd.z() - SD.z());

      sumx += dxdd;          //dxbb + dxcc + dxdd ;
      sumy += dydd;          //dybb + dycc + dydd ;
      sumz += dzdd;          //dzbb + dzcc + dzdd ;
      sumx2 += dxdd * dxdd;  //dxbb*dxbb + dxcc*dxcc + dxdd*dxdd ;
      sumy2 += dydd * dydd;  //dybb*dybb + dycc*dycc + dydd*dydd ;
      sumz2 += dzdd * dzdd;  //dzbb*dzbb + dzcc*dzcc + dzdd*dzdd ;
                             //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      if (0 == i) {
        const HepGeom::Point3D<double> ee(trform * EE0);
        const HepGeom::Point3D<double> ff(trform * FF0);
        const HepGeom::Point3D<double> gg(trform * GG0);
        const HepGeom::Point3D<double> hh(trform * HH0);

        std::cout << "\n\n E discrepancy is " << (ee - extra[0]) << ", " << (ee - extra[0]).mag()
                  << "\n\n F discrepancy is " << (ff - extra[1]) << ", " << (ff - extra[1]).mag()
                  << "\n\n G discrepancy is " << (gg - extra[2]) << ", " << (gg - extra[2]).mag()
                  << "\n\n H discrepancy is " << (hh - extra[3]) << ", " << (hh - extra[3]).mag() << std::endl;

        std::cout << "\n Point SA in CMS coord = " << trforminv * SA << std::endl;
      }
      const HepGeom::Point3D<double> cmsPiv(trforminv * PIVOT);
      psumx += cmsPiv.x();
      psumy += cmsPiv.y();
      psumz += cmsPiv.z();
      psumx2 += cmsPiv.x() * cmsPiv.x();
      psumy2 += cmsPiv.y() * cmsPiv.y();
      psumz2 += cmsPiv.z() * cmsPiv.z();
      const HepGeom::Point3D<double> cmsFlr(trforminv * FLOOR);
      fsumx += cmsFlr.x();
      fsumy += cmsFlr.y();
      fsumz += cmsFlr.z();
      fsumx2 += cmsFlr.x() * cmsFlr.x();
      fsumy2 += cmsFlr.y() * cmsFlr.y();
      fsumz2 += cmsFlr.z() * cmsFlr.z();

      const CLHEP::HepRotation rotInv(trforminv.getRotation());

      std::cout << "\n CMS IP in H4 coord = " << trform * IP << std::endl;
      std::cout << "\n Pivot point in CMS coord = " << cmsPiv << std::endl;
      std::cout << "\n Floor point in CMS coord = " << cmsFlr << std::endl;
      //	 std::cout<<"\n\nTransform rotation is "<<rotInv
      //		  <<"\n\n          translation="<<trforminv.getTranslation() <<std::endl ;

      /*	 std::cout<<"\n"<<names[i]
		  <<"   Theta_x == "<<(blcmsrot*HepGeom::Point3D<double> (1,0,0)).theta()*kRadToDeg
		  <<",  Theta_y == "<<(blcmsrot*HepGeom::Point3D<double> (0,1,0)).theta()*kRadToDeg
		  <<",  Theta_z == "<<(blcmsrot*HepGeom::Point3D<double> (0,0,1)).theta()*kRadToDeg<<std::endl ;

	 std::cout<<"\n"<<names[i]
		  <<"   Phi_x == "<<(blcmsrot*HepGeom::Point3D<double> (1,0,0)).phi()*kRadToDeg
		  <<",  Phi_y == "<<(blcmsrot*HepGeom::Point3D<double> (0,1,0)).phi()*kRadToDeg
		  <<",  Phi_z == "<<(blcmsrot*HepGeom::Point3D<double> (0,0,1)).phi()*kRadToDeg<<std::endl ;
*/

      std::cout << "\n"
                << names[i] << "   Theta_x = " << blcmsrot.thetaX() * kRadToDeg
                << ",  Theta_y = " << blcmsrot.thetaY() * kRadToDeg << ",  Theta_z = " << blcmsrot.thetaZ() * kRadToDeg
                << ",  eta = " << -log(tan(blcmsrot.thetaZ() / 2.)) << std::endl;

      std::cout << "\n"
                << names[i] << "   Phi_x = " << blcmsrot.phiX() * kRadToDeg
                << ",  Phi_y = " << blcmsrot.phiY() * kRadToDeg << ",  Phi_z = " << blcmsrot.phiZ() * kRadToDeg
                << std::endl;

      std::cout << "\n\n BeamLine to CMS translation" << blcmstra << std::endl;
    }
    const double ntot(1.0 * nsur);
    sumx /= ntot;
    sumy /= ntot;
    sumz /= ntot;
    sumx2 = sqrt(sumx2 / ntot - sumx * sumx);
    sumy2 = sqrt(sumy2 / ntot - sumy * sumy);
    sumz2 = sqrt(sumz2 / ntot - sumz * sumz);

    std::cout << "\n Number of survey points = " << ntot;
    std::cout << "\n====== Mean x deviation = " << sumx << " +- " << sumx2;
    std::cout << "\n====== Mean y deviation = " << sumy << " +- " << sumy2;
    std::cout << "\n====== Mean z deviation = " << sumz << " +- " << sumz2 << std::endl;

    psumx /= ntot;
    psumy /= ntot;
    psumz /= ntot;
    psumx2 = sqrt(psumx2 / ntot - psumx * psumx);
    psumy2 = sqrt(psumy2 / ntot - psumy * psumy);
    psumz2 = sqrt(psumz2 / ntot - psumz * psumz);

    std::cout << "\n====== Mean x of pivot = " << psumx << " +- " << psumx2;
    std::cout << "\n====== Mean y of pivot = " << psumy << " +- " << psumy2;
    std::cout << "\n====== Mean z of pivot = " << psumz << " +- " << psumz2 << std::endl;

    fsumx /= ntot;
    fsumy /= ntot;
    fsumz /= ntot;
    fsumx2 = sqrt(fsumx2 / ntot - fsumx * fsumx);
    fsumy2 = sqrt(fsumy2 / ntot - fsumy * fsumy);
    fsumz2 = sqrt(fsumz2 / ntot - fsumz * fsumz);

    std::cout << "\n====== Mean x of floor = " << fsumx << " +- " << fsumx2;
    std::cout << "\n====== Mean y of floor = " << fsumy << " +- " << fsumy2;
    std::cout << "\n====== Mean z of floor = " << fsumz << " +- " << fsumz2 << std::endl;

    const HepGeom::Vector3D<double> hpl[] = {HepGeom::Plane3D<double>(sur[16], sur[0], sur[8]).normal(),
                                             HepGeom::Plane3D<double>(sur[17], sur[1], sur[9]).normal(),
                                             HepGeom::Plane3D<double>(sur[18], sur[2], sur[10]).normal(),
                                             HepGeom::Plane3D<double>(sur[19], sur[3], sur[11]).normal(),
                                             HepGeom::Plane3D<double>(sur[16], sur[12], sur[8]).normal(),
                                             HepGeom::Plane3D<double>(sur[17], sur[13], sur[9]).normal(),
                                             HepGeom::Plane3D<double>(sur[18], sur[14], sur[10]).normal(),
                                             HepGeom::Plane3D<double>(sur[19], sur[15], sur[11]).normal()};

    const HepGeom::Vector3D<double> hplavr(0.125 *
                                           (hpl[0] + hpl[1] + hpl[2] + hpl[3] + hpl[4] + hpl[5] + hpl[6] + hpl[7]));

    std::cout << "\n\n+++ HNormal for A =" << hpl[0].theta() * kRadToDeg << ", " << hpl[0].phi() * kRadToDeg
              << "\n+++ HNormal for B =" << hpl[1].theta() * kRadToDeg << ", " << hpl[1].phi() * kRadToDeg
              << "\n+++ HNormal for C =" << hpl[2].theta() * kRadToDeg << ", " << hpl[2].phi() * kRadToDeg
              << "\n+++ HNormal for D =" << hpl[3].theta() * kRadToDeg << ", " << hpl[3].phi() * kRadToDeg
              << "\n+++ HNormal for A =" << hpl[4].theta() * kRadToDeg << ", " << hpl[4].phi() * kRadToDeg
              << "\n+++ HNormal for B =" << hpl[5].theta() * kRadToDeg << ", " << hpl[5].phi() * kRadToDeg
              << "\n+++ HNormal for C =" << hpl[6].theta() * kRadToDeg << ", " << hpl[6].phi() * kRadToDeg
              << "\n+++ HNormal for D =" << hpl[7].theta() * kRadToDeg << ", " << hpl[7].phi() * kRadToDeg
              << "\n+++ HAverage      =" << hplavr.theta() * kRadToDeg << ", " << hplavr.phi() * kRadToDeg << std::endl;

    const HepGeom::Vector3D<double> vpl[] = {HepGeom::Plane3D<double>(sur[20], sur[28], sur[36]).normal(),
                                             HepGeom::Plane3D<double>(sur[21], sur[29], sur[37]).normal(),
                                             HepGeom::Plane3D<double>(sur[22], sur[30], sur[38]).normal(),
                                             HepGeom::Plane3D<double>(sur[23], sur[31], sur[39]).normal(),

                                             HepGeom::Plane3D<double>(sur[20], sur[24], sur[36]).normal(),
                                             HepGeom::Plane3D<double>(sur[21], sur[25], sur[37]).normal(),
                                             HepGeom::Plane3D<double>(sur[22], sur[26], sur[38]).normal(),
                                             HepGeom::Plane3D<double>(sur[23], sur[27], sur[39]).normal(),

                                             HepGeom::Plane3D<double>(sur[20], sur[32], sur[36]).normal(),
                                             HepGeom::Plane3D<double>(sur[21], sur[33], sur[37]).normal(),
                                             HepGeom::Plane3D<double>(sur[22], sur[34], sur[38]).normal(),
                                             HepGeom::Plane3D<double>(sur[23], sur[35], sur[39]).normal(),

                                             HepGeom::Plane3D<double>(sur[20], sur[32], sur[36]).normal(),
                                             HepGeom::Plane3D<double>(sur[21], sur[33], sur[37]).normal(),
                                             HepGeom::Plane3D<double>(sur[22], sur[34], sur[38]).normal(),
                                             HepGeom::Plane3D<double>(sur[23], sur[35], sur[39]).normal()};

    const HepGeom::Vector3D<double> vplavr(0.0625 * (vpl[0] + vpl[1] + vpl[2] + vpl[3] + vpl[4] + vpl[5] + vpl[6] +
                                                     vpl[7] + vpl[8] + vpl[9] + vpl[10] + vpl[11] + vpl[12] + vpl[13] +
                                                     vpl[14] + vpl[15]));

    std::cout << "\n\n+++ VNormal for A =" << vpl[0].theta() * kRadToDeg << ", " << vpl[0].phi() * kRadToDeg
              << "\n+++ VNormal for B =" << vpl[1].theta() * kRadToDeg << ", " << vpl[1].phi() * kRadToDeg
              << "\n+++ VNormal for C =" << vpl[2].theta() * kRadToDeg << ", " << vpl[2].phi() * kRadToDeg
              << "\n+++ VNormal for D =" << vpl[3].theta() * kRadToDeg << ", " << vpl[3].phi() * kRadToDeg
              << "\n+++ VNormal for A =" << vpl[4].theta() * kRadToDeg << ", " << vpl[4].phi() * kRadToDeg
              << "\n+++ VNormal for B =" << vpl[5].theta() * kRadToDeg << ", " << vpl[5].phi() * kRadToDeg
              << "\n+++ VNormal for C =" << vpl[6].theta() * kRadToDeg << ", " << vpl[6].phi() * kRadToDeg
              << "\n+++ VNormal for D =" << vpl[7].theta() * kRadToDeg << ", " << vpl[7].phi() * kRadToDeg
              << "\n+++ VNormal for A =" << vpl[8].theta() * kRadToDeg << ", " << vpl[8].phi() * kRadToDeg
              << "\n+++ VNormal for B =" << vpl[9].theta() * kRadToDeg << ", " << vpl[9].phi() * kRadToDeg
              << "\n+++ VNormal for C =" << vpl[10].theta() * kRadToDeg << ", " << vpl[10].phi() * kRadToDeg
              << "\n+++ VNormal for D =" << vpl[11].theta() * kRadToDeg << ", " << vpl[11].phi() * kRadToDeg
              << "\n+++ VNormal for A =" << vpl[12].theta() * kRadToDeg << ", " << vpl[12].phi() * kRadToDeg
              << "\n+++ VNormal for B =" << vpl[13].theta() * kRadToDeg << ", " << vpl[13].phi() * kRadToDeg
              << "\n+++ VNormal for C =" << vpl[14].theta() * kRadToDeg << ", " << vpl[14].phi() * kRadToDeg
              << "\n+++ VNormal for D =" << vpl[15].theta() * kRadToDeg << ", " << vpl[15].phi() * kRadToDeg
              << "\n+++ VAverage      =" << vplavr.theta() * kRadToDeg << ", " << vplavr.phi() * kRadToDeg << std::endl;

    unsigned int Hinput(0);
    unsigned int Vinput(0);

    std::cout << "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"
              << "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"
              << "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"
              << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
              << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
              << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
              << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
              << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
              << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
              << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;

    bool keepGoing(true);
    while (keepGoing) {
      std::cout << "\n\n+++ Preparing to make two files for the table position you specify";

      std::cout << "\n\n**Enter desired H position (integer, 25000-300000) [0 to quit]: " << std::ends;
      std::cin >> Hinput;

      keepGoing = (0 != Hinput);

      if (keepGoing) {
        std::cout << "\n\n**Enter desired V position (integer, 5000-30000): " << std::ends;
        std::cin >> Vinput;

        std::cout << "\n\nH input value =" << Hinput << ", V input value =" << Vinput << std::endl;

        if (25000 > Hinput || 300000 < Hinput || 5000 > Vinput || 30000 < Vinput) {
          std::cout << "!!!!!!!! Invalid input(s), try again !!!!!!!" << std::endl;
        } else {
          // calculate transform parameters ================================================

          const HepGeom::Plane3D<double> plaV(sur[20], sur[28], sur[36]);

          const HepGeom::Point3D<double> projPiv(plaV.point(PIVOT));

          const double angleVPerUnit((plaV.point(sur[36]) - projPiv).angle(plaV.point(sur[20]) - projPiv) / 16000.);

          const HepGeom::Plane3D<double> plaH(sur[16], sur[0], sur[8]);

          const HepGeom::Point3D<double> projPivH(plaH.point(PIVOT));

          const double angleHPerUnit((plaH.point(sur[8]) - projPivH).angle(plaH.point(sur[16]) - projPivH) / 176458);
          /*


	       const double angleV1 ( ( plaV.point( sur[24] ) - projPiv ).angle( plaV.point( sur[20] ) - projPiv ) ) ;
	       const double angleV2 ( ( plaV.point( sur[28] ) - projPiv ).angle( plaV.point( sur[24] ) - projPiv ) ) ;
	       const double angleV3 ( ( plaV.point( sur[32] ) - projPiv ).angle( plaV.point( sur[28] ) - projPiv ) ) ;
	       const double angleV4 ( ( plaV.point( sur[36] ) - projPiv ).angle( plaV.point( sur[32] ) - projPiv ) ) ;

	       const double angleH1 ( ( plaH.point( sur[12] ) - projPivH ).angle( plaH.point( sur[16] ) - projPivH ) ) ;
	       const double angleH2 ( ( plaH.point( sur[ 0] ) - projPivH ).angle( plaH.point( sur[12] ) - projPivH ) ) ;
	       const double angleH3 ( ( plaH.point( sur[ 4] ) - projPivH ).angle( plaH.point( sur[ 0] ) - projPivH ) ) ;
	       const double angleH4 ( ( plaH.point( sur[ 8] ) - projPivH ).angle( plaH.point( sur[ 4] ) - projPivH ) ) ;

	       std::cout<<"\n V angles are: "
			<<angleV1*kRadToDeg/4000.<<", "
			<<angleV2*kRadToDeg/4000.<<", "
			<<angleV3*kRadToDeg/4000.<<", "
			<<angleV4*kRadToDeg/4000.<<std::endl ;

	       std::cout<<"\n H angles are: "
			<<angleH1*kRadToDeg/50000.<<", "
			<<angleH2*kRadToDeg/72402.<<", "
			<<angleH3*kRadToDeg/9995.<<", "
			<<angleH4*kRadToDeg/44601.<<std::endl ;

	       std::cout<<"angleHPerUnit="<<kRadToDeg*angleHPerUnit*1.e3<<std::endl ;

	       std::cout<<"angleVPerUnit="<<kRadToDeg*angleVPerUnit*1.e3<<std::endl ;
*/

          const double hAngle((1. * Hinput - 191000.) * angleHPerUnit);

          const double vAngle((1. * Vinput - 16000.) * angleVPerUnit);

          const HepGeom::Rotate3D vRotate(-vAngle, H4ToBeamLine * vpl[0]);
          const HepGeom::RotateY3D hRotate(hAngle);

          const HepGeom::Transform3D myTrf(V3Transform * vRotate * hRotate);

          const CLHEP::HepRotation myRot(myTrf.getRotation());
          const HepGeom::Vector3D<double> myTran(myTrf.getTranslation());

          const HepGeom::RotateZ3D R1(myRot.phiZ() - M_PI);
          const HepGeom::Point3D<double> xUnit(0, 1, 0);
          const HepGeom::Point3D<double> yUnit(-1, 0, 0);
          const HepGeom::Point3D<double> zUnit(0, 0, 1);

          const HepGeom::Transform3D RXRZ(HepGeom::Rotate3D(-myRot.thetaZ(), R1 * xUnit) * R1);

          const HepGeom::Point3D<double> zzUnit(RXRZ * zUnit);
          const HepGeom::Point3D<double> yyUnit(RXRZ * yUnit);
          const HepGeom::Point3D<double> xxUnit(RXRZ * xUnit);
          /*
	       std::cout<<"new z axis = theta "
			<< zzUnit.theta()*kRadToDeg 
			<< ", "
			<< zzUnit.phi()*kRadToDeg 
			<<"\nnew y axis = theta "
			<< yyUnit.theta()*kRadToDeg 
			<< ", "
			<< yyUnit.phi()*kRadToDeg 
			<<"\nnew x axis = theta "
			<< xxUnit.theta()*kRadToDeg 
			<< ", "
			<< xxUnit.phi()*kRadToDeg 
			<< std::endl ;
*/
          const HepGeom::Transform3D RPSI(HepGeom::Transform3D(myRot, HepGeom::Point3D<double>(0, 0, 0)) *
                                          (RXRZ.inverse()));

          CLHEP::Hep3Vector axis;
          double psi;
          RPSI.getRotation().getAngleAxis(psi, axis);

          const double zdiff(axis.angle(zzUnit - HepGeom::Point3D<double>(0, 0, 0)));

          if (0.95 * M_PI < fabs(zdiff))
            psi = -psi;
          /*
	       if( 0.005 < fabs( zdiff ) &&
		   0.005 < fabs( fabs( zdiff ) - M_PI ) )
	       {
		  std::cout<<"**************** BAD 3RD AXIS, zdiff="
			   <<zdiff*kRadToDeg<<" degrees"<<std::endl ;
	       }
	       std::cout<<"\n3rd rotation = "
			<< axis.angle(zzUnit-HepPoint3D(0,0,0))*kRadToDeg 
			<<", "<<angle*kRadToDeg
			<< std::endl ;

	       std::cout<<"PSI="<<(RPSI*xUnit).phi()*kRadToDeg<<std::endl;
	       std::cout<<"xUnitT="<<(RPSI*xUnit)<<std::endl;
	       std::cout<<"RPSI="<<RPSI.getRotation()<<", "<<RPSI.getTranslation()<<std::endl;

//HepRotateZ3D(-0.75*M_PI)*
	       const HepGeom::Transform3D TEMP1 ( HepGeom::RotateX3D(-myRot.thetaZ())*HepGeom::RotateZ3D(M_PI/2.-myRot.phiZ()) ) ;
	       const HepGeom::Transform3D TEMP ( HepGeom::Transform3D(myTrf.getRotation(),HepGeom::Vector3D<double> (0,0,0) )*TEMP1.inverse() ) ;
	       std::cout<<"TEMP="<<TEMP.getRotation()<<", "<<TEMP.getTranslation()<<std::endl;

	       std::cout<<"Trot="<<myTrf.getRotation()<<", "<<myTrf.getTranslation()<<std::endl;
*/
          const double TthetaX(myRot.thetaX());
          const double TthetaY(myRot.thetaY());
          const double TthetaZ(myRot.thetaZ());
          const double TphiX(myRot.phiX());
          const double TphiY(myRot.phiY());
          const double TphiZ(myRot.phiZ());
          const double Txtra(myTran.x());
          const double Tytra(myTran.y());
          const double Tztra(myTran.z());
          const double TetaZ(-log(tan(TthetaZ / 2.)));

          const double Tpsi(psi);

          // write to files ================================================================

          std::stringstream ssTpsir;
          ssTpsir << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << Tpsi;
          const std::string strTpsir(ssTpsir.str());

          std::stringstream ssTxtra;
          ssTxtra << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << Txtra << "*cm";
          const std::string strTxtra(ssTxtra.str());
          std::stringstream ssTytra;
          ssTytra << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << Tytra << "*cm";
          const std::string strTytra(ssTytra.str());
          std::stringstream ssTztra;
          ssTztra << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << Tztra << "*cm";
          const std::string strTztra(ssTztra.str());

          std::stringstream ssTetaZ;
          ssTetaZ << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TetaZ;
          const std::string strTetaZ(ssTetaZ.str());

          std::stringstream ssTthetaX;
          ssTthetaX << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TthetaX * kRadToDeg
                    << "*deg";
          const std::string strTthetaX(ssTthetaX.str());
          std::stringstream ssTphiX;
          ssTphiX << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TphiX * kRadToDeg
                  << "*deg";
          const std::string strTphiX(ssTphiX.str());
          /*	       std::stringstream ssTthetaXr ;
	       ssTthetaXr<<std::setw(7) << setiosflags( std::ios::fixed ) << std::setprecision(5)
		       << TthetaX ;
	       const std::string strTthetaXr ( ssTthetaXr.str() ) ;
	       std::stringstream ssTphiXr ;
	       ssTphiXr<<std::setw(7) << setiosflags( std::ios::fixed ) << std::setprecision(5)
		       << TphiX ;
	       const std::string strTphiXr ( ssTphiXr.str() ) ;
*/
          std::stringstream ssTthetaY;
          ssTthetaY << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TthetaY * kRadToDeg
                    << "*deg";
          const std::string strTthetaY(ssTthetaY.str());
          std::stringstream ssTphiY;
          ssTphiY << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TphiY * kRadToDeg
                  << "*deg";
          const std::string strTphiY(ssTphiY.str());
          /*	       std::stringstream ssTthetaYr ;
	       ssTthetaYr<<std::setw(7) << setiosflags( std::ios::fixed ) << std::setprecision(5)
			 << TthetaY ;
	       const std::string strTthetaYr ( ssTthetaYr.str() ) ;
	       std::stringstream ssTphiYr ;
	       ssTphiYr<<std::setw(7) << setiosflags( std::ios::fixed ) << std::setprecision(5)
		       << TphiY ;
	       const std::string strTphiYr ( ssTphiYr.str() ) ;
*/
          std::stringstream ssTthetaZ;
          ssTthetaZ << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TthetaZ * kRadToDeg
                    << "*deg";
          const std::string strTthetaZ(ssTthetaZ.str());
          std::stringstream ssTphiZ;
          ssTphiZ << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TphiZ * kRadToDeg
                  << "*deg";
          const std::string strTphiZ(ssTphiZ.str());
          std::stringstream ssTphiZr;
          ssTphiZr << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TphiZ;
          const std::string strTphiZr(ssTphiZr.str());

          std::stringstream ssTxpos;
          ssTxpos << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << psumy;
          const std::string strTxpos(ssTxpos.str());

          std::stringstream ssTypos;
          ssTypos << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << -psumx;
          const std::string strTypos(ssTypos.str());

          std::stringstream ssTzpos;
          ssTzpos << std::setw(11) << setiosflags(std::ios::fixed) << std::setprecision(2) << -26733.5 + psumz;
          const std::string strTzpos(ssTzpos.str());
          //-------------- first tbrot file ---------------------------
          std::stringstream ss;
          ss << "tbrot_H" << std::setw(6) << std::setfill('0') << Hinput << "_V" << std::setw(5) << std::setfill('0')
             << Vinput << ".xml" << std::ends;
          const std::string myname(ss.str());

          std::ofstream tbrotfile(myname.c_str());

          tbrotfile << "<?xml version=\"1.0\"?>"
                    << "<DDDefinition xmlns=\"http://www.cern.ch/cms/DDL\""
                    << " xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\""
                    << " xsi:schemaLocation=\"http://www.cern.ch/cms/DDL ../../"
                    << "DetectorDescription/Schema/DDLSchema.xsd\">"
                    << "\n\n<ConstantsSection  label=\"tbrot.xml\" eval=\"true\" >"
                    << "\n\n<Constant name=\"eta\"        value=\"" << strTetaZ
                    << "\"/>"
                    //			<<"\n\n<Constant name=\"theta\"      value=\"2*atan(exp(-[tbrot:eta]))\"/>"
                    << "\n\n<Constant name=\"theta\"      value=\"" << strTthetaZ << "\"/>"
                    << "\n\n<Constant name=\"phi\"        value=\"" << strTphiZ << "\"/>"
                    << "\n\n<Constant name=\"theta_X\"    value=\"" << strTthetaX << "\"/>"
                    << "\n\n<Constant name=\"phi_X\"      value=\"" << strTphiX << "\"/>"
                    << "\n\n<Constant name=\"theta_Y\"    value=\"" << strTthetaY << "\"/>"
                    << "\n\n<Constant name=\"phi_Y\"      value=\"" << strTphiY << "\"/>"
                    << "\n\n<Constant name=\"xtran\"        value=\"" << strTxtra << "\"/>"
                    << "\n<Constant name=\"ytran\"        value=\"" << strTytra << "\"/>"
                    << "\n<Constant name=\"ztran\"        value=\"" << strTztra << "\"/>"
                    << "\n\n</ConstantsSection>"
                    << "\n\n</DDDefinition>" << std::endl;

          tbrotfile.close();

          //--------------- now beam paramter file ----------------------

          std::stringstream ss2;
          ss2 << "ee_beam_H" << std::setw(6) << std::setfill('0') << Hinput << "_V" << std::setw(5) << std::setfill('0')
              << Vinput << "_cff.py" << std::ends;
          const std::string myname2(ss2.str());

          std::ofstream beamfile(myname2.c_str());

          beamfile << "import FWCore.ParameterSet.Config as cms"
                   << "\n\ncommon_beam_direction_parameters = cms.PSet("
                   << "\n    MinEta = cms.double(" << strTetaZ << "),"
                   << "\n    MaxEta = cms.double(" << strTetaZ << "),"
                   << "\n    MinPhi = cms.double(" << strTphiZr << "),"
                   << "\n    MaxPhi = cms.double(" << strTphiZr << "),"
                   << "\n    Psi    = cms.double(" << strTpsir << "),"
                   << "\n    BeamMeanX = cms.double(" << strTxpos << "),"
                   << "\n    BeamMeanY = cms.double(" << strTypos << "),"
                   << "\n    BeamPosition = cms.double(" << strTzpos << ")\n    )" << std::endl;

          beamfile.close();
        }
      }
    }
  }
}

template <class T>
unsigned int CaloGeometryLoaderTest<T>::getDetIdForDDDNode(const DDFilteredView& fv) {
  // perform some consistency checks
  // get the parents and grandparents of this node

  const DDGeoHistory& parents(fv.geoHistory());
  const DDGeoHistory::size_type psize(parents.size());

  EcalBaseNumber baseNumber;
  baseNumber.setSize(psize);

  for (unsigned int i = 1; i <= psize; ++i) {
    baseNumber.addLevel(parents[psize - i].logicalPart().name().name(), parents[psize - i].copyno());
  }

  return baseNumber.getCopyNumber(0);

  //  return m_scheme.getUnitID( baseNumber );
}

#endif