//
// $Id: RecoInterfaces_t.cpp,v 1.13 2002/02/20 22:23:38 mayorov Exp $
//
// File: test_reco_interfaces.cpp
// Purpose: Tests CFT Geometry Reco Interfaces
// Created: 17-DEC-1997   Harry L. Melanson
//
// $Name:  $
//
//

// Include files

#include <iostream>
#include <vector>
#include <iomanip>
#include <stdexcept>

#include "d0om/d0_Ref.hpp"
#include "rcp/RCPManager.hpp"
#include "rcp/RCP.hpp"

#include "geometry_system/management/HandleSystem.hpp"

#include "cft_util/CFTDetector.hpp"

#include "cft_geometry/base/CftBaseGeometry.hpp"

#include "cft_geometry/base/CftGeometer.hpp"
#include "cft_geometry/channel/CftChannelGeometry.hpp"
#include "cft_geometry/channel/CftChannelDoublet.hpp"

#include "cft_geometry/base/CftDoubletSurface.hpp"
#include "cft_geometry/base/CftFiberSurface.hpp"

using namespace std;

using edm::RCPManager;
using edm::RCP;

using namespace dgs;

int main() {

  ////////////////////////////////////////////////////////////////////////
  //
  //  User interfaces to Cft geometry.
  //
  CftGeometer* cft_gm = CftGeometer::get_instance();
  CftChannelGeometry* cft_channel = CftChannelGeometry::get_instance();

  ////////////////////////////////////////////////////////////////////////
  //
  //  The following initializes the CFT geometry.  This will
  //  normally be done by the GeometryManager within the context
  //  of the Reco Framework.  It is done explicitly in this example...
  //
  cout << "Cft Geometry being initialized..." << endl;
  try {
    //
    //  Build default CftBaseGeometry using RCP
    //
    RCPManager* rcpman = RCPManager::instance();
  
    RCP rcp_base = rcpman->extract("cft_geometry","CftBaseGeometry");
    d0_Ref<CftBaseGeometry> cft_base = CftBaseGeometry::build_default(&rcp_base);
    //
    //  Convert to const geometry,
    //
    d0_Ref<const CftBaseGeometry> cft_base_const;
    cft_base_const = D0_CONST_REFCONV((CftBaseGeometry)) (cft_base);
    //
    //  and pass to Geometer.
    //
    cft_gm->refresh(cft_base_const);
    //
    //  Build default channel geometry using RCP.
    //
    RCP rcp_channel = rcpman->extract("cft_geometry","CftChannelGeometry");
    cft_channel->absRefresh(&rcp_channel);
  }
  catch(runtime_error mistake) {
    cerr << "ERROR: " << mistake.what() << endl;
    return(1);
  }

  ////////////////////////////////////////////////////////////////////////
  //
  //  The following code tests the various Reco interfaces.

  //  Get doublet indices

  CFTDoubletList dbls = CFTDetector::all_doublets();
  CFTDoubletList::const_iterator idbl;
  // 
  //  Get pointers to each CftChannelDoublet from the ChannelGeometry.
  //
  map<CFTDoublet, CftChannelDoublet*> doublets = cft_channel->get_doublets();

  cout << "Number of doublets found = " << doublets.size() << endl;

  //
  //  Exercise methods for each doublet...
  //
  for (idbl = dbls.begin(); idbl != dbls.end(); ++idbl) {

    //
    //  Get Handle to CftDoubletSurface's.  Handles are used
    //  so that underlying geometry can be updated as needed.
    //
    Handle<const CftDoubletSurface> surface = doublets[*idbl]->get_doublet_surface_handle();

    cout << "Doublet " << *idbl << std::fixed<<std::setprecision(5)<<endl;
    //
    //  First test the methods belonging to the CftBaseGeometry,
    //  using the CftDoubletSurface Handles...
    //
    cout << "  ";
    cout << "nominal radius = " << std::setw(11)<< surface->nominal_radius() << ", ";
    cout << "zmin = " << std::setw(11)<< surface->zmin() << ", ";
    cout << "zmax = " << std::setw(11)<< surface->zmax();
    cout << endl;

    cout << "  ";
    cout << "inner layer radius = " << std::setw(11)<< surface->inner_radius() << ", ";
    cout << "pitch = " << std::setw(11)<< surface->nominal_inner_pitch();
    cout << endl;

    cout << "  ";
    cout << "outer layer radius = " << std::setw(11)<< surface->outer_radius() << ", ";
    cout << "pitch = " << std::setw(11)<< surface->nominal_outer_pitch();
    cout << endl;
    //
    //  The next set of methods belong to the CftChannelGeometry.
    //
    cout << "  ";
    cout << "nfibers = " << doublets[*idbl]->nfibers() << ", ";
    cout << "dphidz  = " << std::setw(11)<< doublets[*idbl]->dphidz() << ", ";
    cout << "active radius = " << std::setw(11)<< doublets[*idbl]->get_fiber_radius();
    cout << endl;
    //
    //  Check "local to global" methods
    //

    int n = 2 * doublets[*idbl]->nfibers();  //  Total number of fibers on doublet

    //
    //  get_phi(fiber,z):
    //  Check 1st and last fiber on each singlet surface
    //  Check at z = zmin, 0, zmax
    //
    cout << "  Fiber 1";
    cout << ", phi(zmin) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_phi(1,surface->zmin());
    cout << ", phi(0.0) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_phi(1,0.0);
    cout << ", phi(zmax) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_phi(1,surface->zmax());
    cout << endl;

    cout << "  Fiber 2";
    cout << ", phi(zmin) = " <<std::setw(11)<<  
      doublets[*idbl]->get_fiber_phi(2,surface->zmin());
    cout << ", phi(0.0) = " <<std::setw(11)<<  
      doublets[*idbl]->get_fiber_phi(2,0.0);
    cout << ", phi(zmax) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_phi(2,surface->zmax());
    cout << endl;

    cout << "  Fiber n-1";
    cout << ", phi(zmin) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_phi(n-1,surface->zmin());
    cout << ", phi(0.0) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_phi(n-1,0.0);
    cout << ", phi(zmax) = " <<std::setw(11)<<  
      doublets[*idbl]->get_fiber_phi(n-1,surface->zmax());
    cout << endl;

    cout << "  Fiber n";
    cout << ", phi(zmin) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_phi(n,surface->zmin());
    cout << ", phi(0.0) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_phi(n,0.0);
    cout << ", phi(zmax) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_phi(n,surface->zmax());
    cout << endl;

    //
    //  get_r(fiber,z):
    //  Check 1st and last fiber on each singlet surface
    //  Check at z = zmin, 0, zmax
    //
    cout << "  Fiber 1";
    cout << ", r(zmin) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(1,surface->zmin());
    cout << ", r(0.0) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(1,0.0);
    cout << ", r(zmax) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(1,surface->zmax());
    cout << endl;

    cout << "  Fiber 2";
    cout << ", r(zmin) = " <<std::setw(11)<<  
      doublets[*idbl]->get_fiber_r(2,surface->zmin());
    cout << ", r(0.0) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(2,0.0);
    cout << ", r(zmax) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(2,surface->zmax());
    cout << endl;

    cout << "  Fiber n-1";
    cout << ", r(zmin) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(n-1,surface->zmin());
    cout << ", r(0.0) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(n-1,0.0);
    cout << ", r(zmax) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(n-1,surface->zmax());
    cout << endl;

    cout << "  Fiber n";
    cout << ", r(zmin) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(n,surface->zmin());
    cout << ", r(0.0) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(n,0.0);
    cout << ", r(zmax) = " << std::setw(11)<< 
      doublets[*idbl]->get_fiber_r(n,surface->zmax());
    cout << endl;

    cout << "  ";
    cout << "r(phi=0,z=zmin) = " <<std::setw(11)<< 
      surface->get_r(0.0,surface->zmin());
    cout << ", r(phi=0,z=0) = " <<std::setw(11)<< 
      surface->get_r(0.0,0.0);
    cout << ", r(phi=0,z=zmax) = " <<std::setw(11)<< 
      surface->get_r(0.0,surface->zmax());
    cout << endl;

    //
    //  get_r(phi,z):
    //  Check 1st and last fiber on each singlet surface
    //  Check at phi = (1/3)2pi, (2/3)2pi, 2pi
    //
    const double twopi = 6.283185308;
    cout << "  ";
    cout << "r(phi=(1/3)*2pi,z=0) = " <<std::setw(11)<< 
      surface->get_r(twopi/3.0,0.0);
    cout << ", r(phi=(2/3)*2pi,z=0) = " <<std::setw(11)<< 
      surface->get_r(2.0*twopi/3.0,0.0);
    cout << ", r(phi=2pi,z=0) = " <<std::setw(11)<< 
      surface->get_r(twopi,0.0);
    cout << endl;

    //
    //  get_fiber(phi,z)
    //  Check fiber 1 and 2 at z = zmin, 0, zmax
    //
    float phi1 = doublets[*idbl]->get_fiber_phi(1,0.0);
    float phi2 = doublets[*idbl]->get_fiber_phi(2,0.0);
    cout << "  Fiber(phi=" <<std::setw(11)<<  phi1 << ",z=0) = " <<std::setw(11)<< 
      doublets[*idbl]->get_fiber(phi1,0.0);
    cout << ", Fiber(phi=" <<std::setw(11)<<  phi2 << ",z=0) = " <<std::setw(11)<< 
      doublets[*idbl]->get_fiber(phi2,0.0);
    cout << ", Fiber(phi=" <<std::setw(11)<<  (phi1+phi2)/2.0 << ",z=0) = " <<std::setw(11)<< 
      doublets[*idbl]->get_fiber((phi1+phi2)/2.0,0.0);
    cout << endl;

    phi1 = doublets[*idbl]->get_fiber_phi(1,surface->zmin());
    phi2 = doublets[*idbl]->get_fiber_phi(2,surface->zmin());
    cout << "  Fiber(phi=" << std::setw(11)<< phi1 << ",z=zmin) = " <<std::setw(11)<< 
      doublets[*idbl]->get_fiber(phi1,surface->zmin());
    cout << ", Fiber(phi=" <<std::setw(11)<<  phi2 << ",z=zmin) = " <<std::setw(11)<< 
      doublets[*idbl]->get_fiber(phi2,surface->zmin());
    cout << ", Fiber(phi=" << std::setw(11)<< (phi1+phi2)/2.0 << ",z=zmin) = " <<std::setw(11)<< 
      doublets[*idbl]->get_fiber((phi1+phi2)/2.0,surface->zmin());
    cout << endl;

    phi1 = doublets[*idbl]->get_fiber_phi(1,surface->zmax());
    phi2 = doublets[*idbl]->get_fiber_phi(2,surface->zmax());
    cout << "  Fiber(phi=" <<std::setw(11)<<  phi1 << ",z=zmax) = " <<std::setw(11)<< 
      doublets[*idbl]->get_fiber(phi1,surface->zmax());
    cout << ", Fiber(phi=" <<std::setw(11)<<  phi2 << ",z=zmax) = " <<std::setw(11)<< 
      doublets[*idbl]->get_fiber(phi2,surface->zmax());
    cout << ", Fiber(phi=" << std::setw(11)<< (phi1+phi2)/2.0 << ",z=zmax) = " <<std::setw(11)<< 
      doublets[*idbl]->get_fiber((phi1+phi2)/2.0,surface->zmax());
    cout << endl;
    
    cout << "  Fiber 2370";
    cout << ", phi(0.0) = " << 
      doublets[*idbl]->get_phi(2370,0.0);
    cout << endl;


    
    phi1 = doublets[*idbl]->get_fiber_phi(1,0);
    phi2 = doublets[*idbl]->get_fiber_phi(n,0);
    int first=static_cast<int>(doublets[*idbl]->get_fiber(phi1,0.));
    int last=static_cast<int>(doublets[*idbl]->get_fiber(phi2,0.));
    if( first!=1)
      {
	cout << "wrong fiber at phi="<<std::setw(11)<< phi1<<" fiber # "<<first<<" must be 1"<<endl;
	return -1;
      }
    if( last!=n)
      {
	cout << "wrong fiber at phi="<<std::setw(11)<< phi2<<" fiber # "<<last<<" must be "<<n<<endl;
	return -1;
      }

  }

  return 0;
}