//
// $Id: check_ideal.cpp,v 1.5 2000/03/06 15:45:30 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 "d0om/d0_Ref.hpp"
#include "rcp/RCPManager.hpp"
#include "rcp/RCP.hpp"

#include "geometry_system/management/HandleSystem.hpp"

#include "cft_geometry/base/CftBaseGeometry.hpp"

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

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

#include "sftalign/FlatFile.h"

using edm::RCPManager;
using edm::RCP;
using std::vector;
using std::runtime_error;

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;
    RCP rcp_base = rcpman.create("CftBaseGeometry.rcp");
    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.create("CftChannelGeometry.rcp");
    cft_channel->refresh(&rcp_channel);
  }
  catch(runtime_error mistake) {
    cerr << "ERROR: " << mistake.what() << endl;
    exit(1);
  }

  ////////////////////////////////////////////////////////////////////////
  //
  //  Initialize ideal geometry

  FlatFile CFTFile;
  CFTFile.Open("geometry.dat");
  CFTIdealGeometry* FileGeometry = CFTFile.FileGeometry();

  ////////////////////////////////////////////////////////////////////////
  //
  //  The following code tests the various Reco interfaces.
 
  // 
  //  Get pointers to each CftDoublet from the ChannelGeometry.
  //
  vector<CftDoublet*> doublets = cft_channel->get_doublets();
  cout << "Number of doublets REAL  = " << doublets.size() << endl;

  //
  //  And pointers to ideal doublets.
  //
  vector<const CFTIdealGeomDoublet*> ideal;
  for (int i = 1; i <= CFTFile.NumberOfCFTDoublets(); i++)
    ideal.push_back(FileGeometry->DoubletPtr(i));
  cout << "Number of doublets IDEAL = " << CFTFile.NumberOfCFTDoublets() << endl;

  //
  //  Now get Handles to CftDoubletSurface's.  Handles are used
  //  so that underlying geometry can be updated as needed.
  //
  vector< Handle<const CftDoubletSurface> > surfaces;
  for (int i = 0; i < doublets.size(); i++)
    surfaces.push_back(doublets[i]->get_doublet_surface_handle());


  //
  //  Exercise methods for each doublet...
  //
  for (int i = 0; i < doublets.size(); i++) {
    cout << "Doublet " << i+1 << endl;
    //
    //  First test the methods belonging to the CftBaseGeometry,
    //  using the CftDoubletSurface Handles...
    //
    cout << " nominal radius: DIFF = " << 
      surfaces[i]->nominal_radius() - ideal[i]->nominal_radius();
    cout << " REAL  = " << 
      surfaces[i]->nominal_radius();
    cout << " IDEAL = " << 
      ideal[i]->nominal_radius();
    cout << endl;

    cout << " zmin: DIFF = " << 
      surfaces[i]->zmin() - ideal[i]->zmin();
    cout << " REAL  = " << 
      surfaces[i]->zmin();
    cout << " IDEAL = " << 
      ideal[i]->zmin();
    cout << endl;

    cout << " zmax: DIFF = " << 
      surfaces[i]->zmax() - ideal[i]->zmax();
    cout << " REAL  = " << 
      surfaces[i]->zmax();
    cout << " IDEAL = " << 
      ideal[i]->zmax();
    cout << endl;

    cout << " inner radius: DIFF = " << 
      surfaces[i]->inner_radius() - ideal[i]->inner_radius();
    cout << " REAL  = " << 
      surfaces[i]->inner_radius();
    cout << " IDEAL = " << 
      ideal[i]->inner_radius();
    cout << endl;

    cout << " outer radius: DIFF = " << 
      surfaces[i]->outer_radius() - ideal[i]->outer_radius();
    cout << " REAL  = " << 
      surfaces[i]->outer_radius();
    cout << " IDEAL = " << 
      ideal[i]->outer_radius();
    cout << endl;

    cout << " nominal pitch: DIFF = " << 
      doublets[i]->get_fiber_pitch() - 
      ideal[i]->get_fiber_pitch();
    cout << " REAL = " << 
      doublets[i]->get_fiber_pitch();
    cout << " IDEAL  = " << 
      ideal[i]->get_fiber_pitch();
    cout << endl;

    //
    //  The next set of methods belong to the CftChannelGeometry.
    //
    cout << " nfibers: DIFF = " << 
      doublets[i]->nfibers() - ideal[i]->nfibers();
    cout << " REAL  = " << 
      doublets[i]->nfibers();
    cout << " IDEAL = " << 
      ideal[i]->nfibers();
    cout << endl;

    cout << " dphidz: DIFF = " << 
      doublets[i]->dphidz() - ideal[i]->dphidz();
    cout << " REAL  = " << doublets[i]->dphidz();
    cout << " IDEAL = " << ideal[i]->dphidz();
    cout << endl;

    cout << " active radius: DIFF = " << 
      doublets[i]->get_fiber_radius() - ideal[i]->get_fiber_radius();
    cout << " REAL  = " << 
      doublets[i]->get_fiber_radius();
    cout << " IDEAL = " << 
      ideal[i]->get_fiber_radius();
    cout << endl;
    //
    //  Check "local to global" methods
    //

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

    //
    //  get_phi(fiber,z):
    //  Check 1st and 2nd fiber
    //  Check at z = zmin, 0, zmax
    //
    cout << " Fiber 1, phi(zmin) DIFF = " << 
      doublets[i]->get_phi(1,surfaces[i]->zmin()) - 
      ideal[i]->get_phi(1,ideal[i]->zmin());
    cout << " REAL = " << 
      doublets[i]->get_phi(1,surfaces[i]->zmin());
    cout << " IDEAL = " << 
      ideal[i]->get_phi(1,ideal[i]->zmin());
    cout << endl;

    cout << " Fiber 1, phi(0.0) DIFF = " << 
      doublets[i]->get_phi(1,0.0) - 
      ideal[i]->get_phi(1,0.0);
    cout << " REAL = " << 
      doublets[i]->get_phi(1,0.0);
    cout << " IDEAL = " << 
      ideal[i]->get_phi(1,0.0);
    cout << endl;

    cout << " Fiber 1, phi(zmax) DIFF = " << 
      doublets[i]->get_phi(1,surfaces[i]->zmax()) - 
      ideal[i]->get_phi(1,ideal[i]->zmax());
    cout << " REAL = " << 
      doublets[i]->get_phi(1,surfaces[i]->zmax());
    cout << " IDEAL = " << 
      ideal[i]->get_phi(1,ideal[i]->zmax());
    cout << endl;

    cout << " Fiber 2, phi(zmin) DIFF = " << 
      doublets[i]->get_phi(2,surfaces[i]->zmin()) - 
      ideal[i]->get_phi(2,ideal[i]->zmin());
    cout << " REAL = " << 
      doublets[i]->get_phi(2,surfaces[i]->zmin());
    cout << " IDEAL = " << 
      ideal[i]->get_phi(2,ideal[i]->zmin());
    cout << endl;

    cout << " Fiber 2, phi(0.0) DIFF = " << 
      doublets[i]->get_phi(2,0.0) - 
      ideal[i]->get_phi(2,0.0);
    cout << " REAL = " << 
      doublets[i]->get_phi(2,0.0);
    cout << " IDEAL = " << 
      ideal[i]->get_phi(2,0.0);
    cout << endl;

    cout << " Fiber 2, phi(zmax) DIFF = " << 
      doublets[i]->get_phi(2,surfaces[i]->zmax()) - 
      ideal[i]->get_phi(2,ideal[i]->zmax());
    cout << " REAL = " << 
      doublets[i]->get_phi(2,surfaces[i]->zmax());
    cout << " IDEAL = " << 
      ideal[i]->get_phi(2,ideal[i]->zmax());
    cout << endl;

    //
    //  get_r(fiber,z):
    //  Check 1st and 2nd fiber
    //  Check at z = zmin, 0, zmax
    //
    cout << " Fiber 1, r(zmin) DIFF = " << 
      doublets[i]->get_r(1,surfaces[i]->zmin()) - 
      ideal[i]->get_r(1,ideal[i]->zmin());
    cout << " REAL = " << 
      doublets[i]->get_r(1,surfaces[i]->zmin());
    cout << " IDEAL = " << 
      ideal[i]->get_r(1,ideal[i]->zmin());
    cout << endl;

    cout << " Fiber 1, r(z=0.0) DIFF = " << 
      doublets[i]->get_r(1,0.0) - 
      ideal[i]->get_r(1,0.0);
    cout << " REAL = " << 
      doublets[i]->get_r(1,0.0);
    cout << " IDEAL = " << 
      ideal[i]->get_r(1,0.0);
    cout << endl;

    cout << " Fiber 1, r(zmax) DIFF = " << 
      doublets[i]->get_r(1,surfaces[i]->zmax()) - 
      ideal[i]->get_r(1,ideal[i]->zmax());
    cout << " REAL = " << 
      doublets[i]->get_r(1,surfaces[i]->zmax());
    cout << " IDEAL = " << 
      ideal[i]->get_r(1,ideal[i]->zmax());
    cout << endl;

    cout << " Fiber 2, r(zmin) DIFF = " << 
      doublets[i]->get_r(2,surfaces[i]->zmin()) - 
      ideal[i]->get_r(2,ideal[i]->zmin());
    cout << " REAL = " << 
      doublets[i]->get_r(2,surfaces[i]->zmin());
    cout << " IDEAL = " << 
      ideal[i]->get_r(2,ideal[i]->zmin());
    cout << endl;

    cout << " Fiber 2, r(z=0.0) DIFF = " << 
      doublets[i]->get_r(2,0.0) - 
      ideal[i]->get_r(2,0.0);
    cout << " REAL = " << 
      doublets[i]->get_r(2,0.0);
    cout << " IDEAL = " << 
      ideal[i]->get_r(2,0.0);
    cout << endl;

    cout << " Fiber 2, r(zmax) DIFF = " << 
      doublets[i]->get_r(2,surfaces[i]->zmax()) - 
      ideal[i]->get_r(2,ideal[i]->zmax());
    cout << " REAL = " << 
      doublets[i]->get_r(2,surfaces[i]->zmax());
    cout << " IDEAL = " << 
      ideal[i]->get_r(2,ideal[i]->zmax());
    cout << endl;

     cout << " r(phi=0,z=zmin): DIFF = " <<
       surfaces[i]->get_r(0.0,surfaces[i]->zmin()) -
       ideal[i]->get_r(0.0,ideal[i]->zmin());
     cout << " REAL = " <<
       surfaces[i]->get_r(0.0,surfaces[i]->zmin());
     cout << " IDEAL = " <<
       ideal[i]->get_r(0.0,ideal[i]->zmin());
     cout << endl;

     cout << " r(phi=0,z=0.0): DIFF = " <<
       surfaces[i]->get_r(0.0,0.0) -
       ideal[i]->get_r(0.0,0.0);
     cout << " REAL = " <<
       surfaces[i]->get_r(0.0,0.0);
     cout << " IDEAL = " <<
       ideal[i]->get_r(0.0,0.0);
     cout << endl;

     cout << " r(phi=0,z=zmax): DIFF = " <<
       surfaces[i]->get_r(0.0,surfaces[i]->zmax()) -
       ideal[i]->get_r(0.0,ideal[i]->zmax());
     cout << " REAL = " <<
       surfaces[i]->get_r(0.0,surfaces[i]->zmax());
     cout << " IDEAL = " <<
       ideal[i]->get_r(0.0,ideal[i]->zmax());
     cout << endl;

     //
     //  get_fiber(phi,z)
     //  Check fiber 1 and 2 at z = zmin, 0, zmax
     //
     float phi1 = doublets[i]->get_phi(1,0.0); 
     float phi2 = doublets[i]->get_phi(2,0.0);

     cout << " Fiber(phi=" << phi1 << ",z=0) DIFF = " <<
       doublets[i]->get_fiber(phi1,0.0) - 
       ideal[i]->get_fiber(phi1,0.0);
     cout << " REAL = " <<
       doublets[i]->get_fiber(phi1,0.0);
     cout << " IDEAL = " <<
       ideal[i]->get_fiber(phi1,0.0);
     cout << endl;

     cout << " Fiber(phi=" << phi2 << ",z=0) DIFF = " <<
       doublets[i]->get_fiber(phi2,0.0) - 
       ideal[i]->get_fiber(phi2,0.0);
     cout << " REAL = " <<
       doublets[i]->get_fiber(phi2,0.0);
     cout << " IDEAL = " <<
       ideal[i]->get_fiber(phi2,0.0);
     cout << endl;

     cout << " Fiber(phi=" << (phi1+phi2)/2.0 << ",z=0) DIFF = " <<
       doublets[i]->get_fiber((phi1+phi2)/2.0,0.0) - 
       ideal[i]->get_fiber((phi1+phi2)/2.0,0.0);
     cout << " REAL = " <<
       doublets[i]->get_fiber((phi1+phi2)/2.0,0.0);
     cout << " IDEAL = " <<
       ideal[i]->get_fiber((phi1+phi2)/2.0,0.0);
     cout << endl;

  }
}