// -*- C++ -*-
//
// Package:     DRGeomProcessor
// Module:      DRGeomProcessor
// 
// Description: Example of a Processor using DRLayer (to be named  ADRStore) 
//
// Implementation:
//     <Notes on implementation>
//
// Author:      Inga Karliner
// Created:     Sun Jun 21 18:55:00 EDT 1998
// $Id$
//
// Revision history
//
// $Id$
//

#include "Experiment/Experiment.h"

// system include files
#include <iomanip.h> 
#include <fstream.h>
#if defined(AMBIGUOUS_STRING_FUNCTIONS_BUG)
#include <string>
#endif             
#if defined(STL_TEMPLATE_DEFAULT_PARAMS_FIRST_BUG)
#endif /* STL_TEMPLATE_DEFAULT_PARAMS_FIRST_BUG */

// user include files
#include "DRGeomProcessor/DRGeomProcessor.h"
#include "Experiment/report.h"


#include "FrameAccess/FAItem.h"
#include "FrameAccess/FATable.h"
#include "DataHandler/Frame.h"
#include "FrameAccess/extract.h"         
#include "DataHandler/Record.h"

#include "ADRGeom/ADRSenseWireStore.h"
#include "ActiveElements/AEWireInfo.h"
#include "ActiveElements/AEWireLayerCylindrical.h"
//#include "ActiveElements/AEWireLayerRealistic.h"
#include "ADRGeom/ADRWireAddress.h"

// STL classes

//
// constants, enums and typedefs
//
static const char* const kFacilityString = "Processor.DRGeomProcessor" ;

//
// static data member definitions
//

//
// constructors and destructor
//
DRGeomProcessor::DRGeomProcessor( void )               // anal1
   : Processor( "DRGeomProcessor" )
{
   report( INFO, kFacilityString ) << "here in ctor()" << endl;

   // ---- bind a method to a stream -----
   // These lines ARE VERY IMPORTANT! If you don't bind the 
   // code you've just written (the "action") to a stream, 
   // your code won't get executed!

   //   bind( &DRGeomProcessor::event,    Stream::kEvent );
   //   bind( &DRGeomProcessor::beginRun, Stream::kBeginRun );
   bind( &DRGeomProcessor::baseGeometry,Stream::kBaseGeometry); 
   //   bind( &DRGeomProcessor::drAlignment, Stream::kDRAlignment );
   //   bind( &DRGeomProcessor::endRun,   Stream::kEndRun );

   // do anything here that needs to be done at creation time
   // (e.g. allocate resources etc.)

}

// DRGeomProcessor::DRGeomProcessor( const DRGeomProcessor& iDRGeomProcessor )
// {
//    *this = iDRGeomProcessor;
// }

DRGeomProcessor::~DRGeomProcessor()                    // anal5
{
   report( INFO, kFacilityString ) << "here in dtor()" << endl;
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}

//
// assignment operators
//
// const DRGeomProcessor& DRGeomProcessor::operator=( const DRGeomProcessor& iDRGeomProcessor )
// {
//   if( this != &iDRGeomProcessor ) {
//      // do actual copying here, plus:
//      // "SuperClass"::operator=( iDRGeomProcessor );
//   }
//
//   return *this;
// }

//
// member functions
//

// ------------ methods for beginning/end "Interactive" ------------
// --------------------------- init method -------------------------
void
DRGeomProcessor::init( void )          // anal1 "Interactive"
{
   report( INFO, kFacilityString ) << "here in init()" << endl;

   // do any initialization here based on Parameter Input by User
   // (e.g. run expensive algorithms that are based on parameters
   //  specified by user at run-time)

}

// -------------------- terminate method ----------------------------
void
DRGeomProcessor::terminate( void )     // anal5 "Interactive"
{
   report( INFO, kFacilityString ) << "here in terminate()" << endl;

   // do anything here BEFORE New Parameter Change
   // (e.g. write out result based on parameters from user-input)
 
}

// ---------------- standard place to book histograms ---------------
void
DRGeomProcessor::hist_book( TBHistoManager& )
{
   report( INFO, kFacilityString ) << "here in hist_book()" << endl;

   // book your histograms here

}

// --------------------- methods bound to streams -------------------
//ActionBase::ActionResult
//DRGeomProcessor::event( Frame& )          // anal3 equiv.
//{
//   report( INFO, kFacilityString ) << "here in event()" << endl;
//
//   return ActionBase::kPassed;
//}

//ActionBase::ActionResult
//DRGeomProcessor::beginRun( Frame& )       // anal2 equiv.
//{
//   report( INFO, kFacilityString ) << "here in beginRun()" << endl;
//
//   return ActionBase::kPassed;
//}


/*
ActionBase::ActionResult
DRGeomProcessor::endRun( Frame& )         // anal4 equiv.
{
   report( INFO, kFacilityString ) << "here in endRun()" << endl;

   return ActionBase::kPassed;
}
*/

ActionBase::ActionResult
DRGeomProcessor::drAlignment( Frame& iFrame )       // anal2 equiv.
{
//  bind to DRAlignment
  report( INFO, kFacilityString ) << "here in DRAlign" << endl;
  //   FAItem< ADRSenseWireStore > iADRSenseWireStore ;
  //   extract( iFrame.record( Stream::kDRAlignment ), iADRSenseWireStore );
  return ActionBase::kPassed;
}

ActionBase::ActionResult
DRGeomProcessor::baseGeometry( Frame& iFrame )       // anal2 equiv.
{ 
  //  bind to baseGeometry,till Mike implements drAlignment
  report( INFO, kFacilityString ) << "here in baseGeometry" << endl;
  //extract ADRSenseWireStore
  FAItem< ADRSenseWireStore > iADRSenseWireStore ;
  extract( iFrame.record( Stream::kBaseGeometry ), iADRSenseWireStore );
  if( iADRSenseWireStore.valid()  )
    {  
       // EXAMPLES OF ACTIVE DR_GEOMETRY ACCESS FUNCTIONS

       //----------------------------+
       //Global Active DR properties |
       //----------------------------+

  ofstream outf("ADrGeomProcessor.dump");

  cout << endl<<"*** ADRGeom Processor ***** " 
       << endl<<"============================"<< endl; 
  cout << "#DR Layers        "<< iADRSenseWireStore->numberOfLayers()<<endl;
  cout << "#DR Layers Axial  "<< iADRSenseWireStore->numberOfAxialLayers()<< endl;
  cout << "#DR Layers Stereo "<< iADRSenseWireStore->numberOfStereoLayers()<< endl;
  cout << "#DR Layers U      "<< iADRSenseWireStore->numberOfULayers()<< endl;
  cout << "#DR Layers V      "<< iADRSenseWireStore->numberOfVLayers()<< endl;
  cout << "#DR Wires Total   "<< iADRSenseWireStore->numberOfWiresTotal()<<endl;

  outf << endl<<"*** ADRGeom Processor ***** "
       << endl<<"============================"<< endl; 
  outf << "#DR Layers        "<< iADRSenseWireStore->numberOfLayers()<<endl;
  outf << "#DR Layers Axial  "<< iADRSenseWireStore->numberOfAxialLayers()<< endl;
  outf << "#DR Layers Stereo "<< iADRSenseWireStore->numberOfStereoLayers()<< endl;
  outf << "#DR Layers U      "<< iADRSenseWireStore->numberOfULayers()<< endl;
  outf << "#DR Layers V      "<< iADRSenseWireStore->numberOfVLayers()<< endl;
  outf << "#DR Wires Total   "<< iADRSenseWireStore->numberOfWiresTotal()<<endl;

       //       unsigned int itLayer  = 1;
       //       unsigned int itLayer  = 16;
       //                     unsigned int itLayer  = 17;
  unsigned int itLayer  = 30;
  unsigned int itLayer2  = 47;
  unsigned int iwire ;
  const unsigned int* alayer = iADRSenseWireStore->aLayers();
  const unsigned int* ulayer = iADRSenseWireStore->uLayers();
  const unsigned int* vlayer = iADRSenseWireStore->vLayers();
  
  int ika = iADRSenseWireStore->numberOfAxialLayers() -1;
  int iku = iADRSenseWireStore->numberOfULayers()-1; 
  int ikv = iADRSenseWireStore->numberOfVLayers()-1;
  cout << endl<< "First/Last Axial layers:"
       << setw(3) << alayer[0] <<setw(3) << alayer[ika]<<endl
       <<"           U     layers:"
       <<setw(3) << ulayer[0] <<setw(3) << ulayer[iku]
       << endl 
       <<"           V     layers:"<<setw(3) << vlayer[0] <<setw(3) << vlayer[ikv]
       <<endl;
  cout <<endl
       <<"+----------------------------------------+"<<endl
       <<"| CellID  <------> Layer, (Wire in Layer)|"<<endl
       <<"| Trackers Beware and NOTE Range:        |"<<endl
       <<"|  CellID     =1:9796                    |"<<endl
       <<"|  Layer      =1:  47                    |"<<endl
       <<"|  WireInLayer=0:(NWiresInLayer-1)       |"<<endl
       <<"+----------------------------------------+"<<endl;

  outf << endl<< "First/Last Axial layers:"
       << setw(3) << alayer[0] <<setw(3) << alayer[ika]<<endl
       <<"           U     layers:"<<setw(3) 
       << ulayer[0] <<setw(3) << ulayer[iku] << endl 
       <<"           V     layers:"<<setw(3) 
       << vlayer[0] <<setw(3) << vlayer[ikv] <<endl;

  outf <<endl
       <<"+----------------------------------------+"<<endl
       <<"| CellID  <------> Layer, (Wire in Layer)|"<<endl
       <<"| Trackers Beware and NOTE Range:        |"<<endl
       <<"|  CellID     =1:9796                    |"<<endl
       <<"|  Layer      =1:  47                    |"<<endl
       <<"|  WireInLayer=0:(NWiresInLayer-1)       |"<<endl
       <<"+----------------------------------------+"<<endl;  
       //-------------------------------------------------------------+
       //CellID-> Layer, Wire in Layer                                |      
       //    Trackers Beware NOTE Range:                              |
       //    CellID=1:9796  Layer=1:47 WireInLayer=0:(NWiresInLayer-1)|
       //-------------------------------------------------------------+
  UInt32 CellID = 1 ;
  //       UInt32 CellID = 9796 ;
  ADRSenseWireStore::LayerWire lw=iADRSenseWireStore->layerWire(CellID);
  unsigned short layerCell=lw.first;
  unsigned short wireCell =lw.second;
  cout<<endl<<"CellID "<<CellID<<"->Layer, Wire:"
      <<setw(3)<< layerCell<<setw(4)<<wireCell<<endl;
  outf<<endl<<"CellID "<<CellID<<"->Layer, Wire:"
      <<setw(3)<< layerCell<<setw(4)<<wireCell<<endl;       
       //-------------------------------------------------------------+
       //Layer, Wire in Layer -> CellID                               |      
       //    Trackers Beware NOTE Range:                              |
       //    CellID=1:9796  Layer=1:47  WireInLayer=0:NWiresInLayer-1 |
       //-------------------------------------------------------------+
  layerCell = 18;
  wireCell  = 22;
  cout<<endl
      <<"Layer,Wire="<< setw(3)<<layerCell<<setw(3)<<wireCell
      <<"->CellID="  <<iADRSenseWireStore->cellID(layerCell,wireCell)<<endl;
  outf<<endl
      <<"Layer,Wire="<< setw(3)<<layerCell<<setw(3)<<wireCell
      <<"->CellID="  <<iADRSenseWireStore->cellID(layerCell,wireCell)<<endl; 
  
       //------------------------------+
       // Layer Properties Access      |
       // access AEWireLayerCylindrical| 
       // via ADRSenseWireStore        |
       //------------------------------+
// get layer # 
  const AEWireLayerCylindrical& layer2 = iADRSenseWireStore->layer(itLayer2);    
  int layernumber =  iADRSenseWireStore->layerNumber(layer2);
  cout<<"*************************"<<endl;
  cout<<"Layer" <<setw(3)<<itLayer2<<" is ="
      <<setw(3)<<layernumber<<endl;
  cout<<"*************************"<<endl;
  
  const AEWireLayerCylindrical& layer = iADRSenseWireStore->layer(itLayer);
//**       const AEWireLayerRealistic& layer 
//**                = iADRSenseWireStore->layer(itLayer);
   unsigned int numberOfWiresLyr 
     = iADRSenseWireStore->numberOfWiresLyr(itLayer);
   double radiusLayerEndCap      = layer.radiusEndCap();
   double phi0West               = layer.phiOfFirstWireInVolume2();
   double phi0East               = layer.phiOfFirstWireInVolume1();
   double radiusLayer            = iADRSenseWireStore->radius(itLayer);
   double phi0                   = iADRSenseWireStore->phi0(itLayer);
   double phiCellWidth           = iADRSenseWireStore->cellPhiWidth(itLayer);
   double tanStereoAngle         = iADRSenseWireStore->tanStereoAngle(itLayer);
   double zMin                   = iADRSenseWireStore->zMin(itLayer);
   double zMax                   = iADRSenseWireStore->zMax(itLayer);
   DABoolean isAxial             = iADRSenseWireStore->isAxial(itLayer);
   DABoolean isStereo            = iADRSenseWireStore->isStereo(itLayer);
   DABoolean isStereoU           = iADRSenseWireStore->isStereoU(itLayer);
   DABoolean isStereoV           = iADRSenseWireStore->isStereoV(itLayer);
/*
  unsigned int numberOfWiresLyr = layer.numberOfWires();
  double radiusLayerEndCap      = layer.radiusEndCap();
  double phi0West               = layer.phiOfFirstWireInVolume2();
  double phi0East               = layer.phiOfFirstWireInVolume1();
  double radiusLayer            = layer.radius();
  double phi0                   = layer.phi0();
  double phiCellWidth           = layer.cellPhiWidth();
  double tanStereoAngle         = layer.tanStereoAngle();
  double zMin                   = layer.zMin();
  double zMax                   = layer.zMax();
  DABoolean isAxial             = layer.isAxial();
  DABoolean isStereo            = layer.isStereo();
  DABoolean isStereoU           = layer.isStereoU();
  DABoolean isStereoV           = layer.isStereoV();
  */
//get radius at the layer's sense wire at z=0 and one end, compare with input
    DABoolean successR;
    double rMin =  iADRSenseWireStore->radiusAtz(successR,itLayer,zMin);
    //    double rMin = layer.radiusAtz(successR,zMin);
    if(!successR) {cout<<"? rMin="<<rMin<<endl;} 
    double rMax =  iADRSenseWireStore->radiusAtz(successR,itLayer,zMax);
    //    double rMax = layer.radiusAtz(successR, zMax);
    if(!successR) {cout<<"? rMax="<<rMax<<endl;}
    double r0   =  iADRSenseWireStore->radiusAtz(successR,itLayer,0.); 
    //    double r0   = layer.radiusAtz(successR, 0.);
    if(!successR) {cout<<"? r0  ="<<r0  <<endl; }
    double r1   =  iADRSenseWireStore->radiusAtz(successR,itLayer,zMax/2.);
    //    double r1   = layer.radiusAtz(successR, zMax/2.);
    if(!successR) {cout<<"? r1  ="<<r1  <<endl; }

//get LAYER'S radMin, radMax  at z=z0, zMin
    double rLAYERMin0 = r0 - layer.radialExtentMinus();
    double rLAYERMax0 = r0 + layer.radialExtentPlus();
    double rLAYERMin1 = rMin - layer.radialExtentMinus();
    double rLAYERMax1 = rMin + layer.radialExtentPlus();
    cout <<endl<<"Layer Access Example:" 
	          <<"====================="
	 <<endl<<"  Layer                     = "<< itLayer
	 <<endl<<"  NumberOfWires             = "<< numberOfWiresLyr
	 <<endl<<"  Radius Layer (z=0)        = "<< radiusLayer     
	 <<endl<<"  Phi0 First   (z=0)        = "<< phi0
	 <<endl<<"  RadiusEast/West local     = "<< radiusLayerEndCap      
	 <<endl<<"  Phi0 First Wire West local= "<< phi0West     
	 <<endl<<"  Phi0 First Wire East local= "<< phi0East
	 <<endl<<"  Phi Cell Width            = "<< phiCellWidth 
	 <<endl<<"  tan(StereoAngle)          = "<< tanStereoAngle
	 <<endl<<"  zmin                      = "<< zMin
	 <<endl<<"  zmax                      = "<< zMax
	 <<endl<<"  r(zmin)                   = "<< rMin
	 <<endl<<"  r(zmax)                   = "<< rMax
	 <<endl<<"  r(z=0)                    = "<< r0
	 <<endl<<"  r(z=zmax/2)               = "<< r1
	 <<endl<<"  Min LAYER radius z=0      = "<< rLAYERMin0 
	 <<endl<<"  Max LAYER radius z=0      = "<< rLAYERMax0 
	 <<endl<<"  Min LAYER radius z=zmin   = "<< rLAYERMin1 
	 <<endl<<"  Max LAYER radius z=zmin   = "<< rLAYERMax1 
	 <<endl;
    outf <<endl<<"Layer Access Example:" 
	 <<"====================="
	 <<endl<<"  Layer                     = "<< itLayer
	 <<endl<<"  NumberOfWires             = "<< numberOfWiresLyr
	 <<endl<<"  Radius Layer (z=0)        = "<< radiusLayer     
	 <<endl<<"  Phi0 First   (z=0)        = "<< phi0
	 <<endl<<"  RadiusEast/West local     = "<< radiusLayerEndCap      
	 <<endl<<"  Phi0 First Wire West local= "<< phi0West     
	 <<endl<<"  Phi0 First Wire East local= "<< phi0East
	 <<endl<<"  Phi Cell Width            = "<< phiCellWidth 
	 <<endl<<"  tan(StereoAngle)          = "<< tanStereoAngle
	 <<endl<<"  zmin                      = "<< zMin
	 <<endl<<"  zmax                      = "<< zMax
	 <<endl<<"  r(zmin)                   = "<< rMin
	 <<endl<<"  r(zmax)                   = "<< rMax
	 <<endl<<"  r(z=0)                    = "<< r0
	 <<endl<<"  r(z=zmax/2)               = "<< r1
	 <<endl<<"  Min LAYER radius z=0      = "<< rLAYERMin0 
	 <<endl<<"  Max LAYER radius z=0      = "<< rLAYERMax0 
	 <<endl<<"  Min LAYER radius z=zmin   = "<< rLAYERMin1 
	 <<endl<<"  Max LAYER radius z=zmin   = "<< rLAYERMax1 
	 <<endl;
   
    if (isAxial) 
      { cout <<endl <<"Layer"<<itLayer<<" is Axial"<<endl;}
    else
      { cout <<"Layer "<<itLayer<<" is Not Axial"<<endl;}
    if (isStereo) 
      { cout <<"Layer "<<itLayer<<" is Stereo"<<endl;
      if (isStereoU) 
	{ cout <<"Layer is U-Stereo"<<endl;}
      else
	{ cout <<"Layer is Not U-Stereo"<<endl;}
      if (isStereoV) 
	{ cout <<"Layer is V-Stereo"<<endl;}
      else
	{ cout <<"Layer is Not V-Stereo"<<endl;}
      }
    else
      { cout <<"Layer"<<itLayer<<" is Not Stereo"<<endl;}       
    
    outf<<endl
	<<"Layer NWires  R [meter]  Phi0         R[meter]  Phi0         PhiCell   tan"<< endl
	<<"              WestEndCap                at z=0               Width       (stereoAngle)"<< endl
	<<"              (local)"         
	<<endl;
    for( itLayer=1;   itLayer <=  iADRSenseWireStore->numberOfLayers(); ++itLayer )      
      {
	const AEWireLayerCylindrical& layer = iADRSenseWireStore->layer(itLayer);
//**    const AEWireLayerRealistic& layer = iADRSenseWireStore->layer(itLayer);
/*         outf<<setw( 5) <<itLayer
	   <<setw( 5) <<layer.numberOfWires()
	   <<setw(12) <<layer.radiusEndCap()
	   <<setw(12) <<layer.phiOfFirstWireInVolume2()
	   <<setw(12) <<layer.radius()
	   <<setw(12) <<layer.phi0()
	   <<setw(12) <<layer.cellPhiWidth()
	   <<setw(12) <<layer.tanStereoAngle()<<endl; */
	  outf<<setw( 5) <<itLayer
	      <<setw( 5) << iADRSenseWireStore->numberOfWiresLyr(itLayer)
	      <<setw(12) <<layer.radiusEndCap()
	      <<setw(12) << layer.phiOfFirstWireInVolume2()
	      <<setw(12) << iADRSenseWireStore->radius(itLayer)
	      <<setw(12) << iADRSenseWireStore->phi0(itLayer)
	      <<setw(12) << iADRSenseWireStore->cellPhiWidth(itLayer)
	      <<setw(12) << iADRSenseWireStore->tanStereoAngle(itLayer)<<endl;
      }
    double xP,yP,zP;
    HepPoint3D position;
    double rP;
    unsigned int irlayer;
    // 
    cout<<"======================="<<endl
	<<"find layer for 3D point"<<endl
	<<"======================="<<endl;
    DABoolean layer_success;
    xP = .1; yP = .15; zP = .13;
    position = HepPoint3D(xP,yP,zP);
    rP   = position.perp();
    irlayer = iADRSenseWireStore->layerNumber(layer_success,position);
    if (layer_success)
      { 
	cout<<endl<<"layerNumber success"<<endl;  
      }
    else	 
      { cout<<endl<<"layerNumber fails"<<endl; 
      }
    cout<<"DRGeom Processor layerNumber(position) ="<<irlayer<<endl;
    cout<<" sense wire layer for point"
	<< setw(10)<<position.x()<<setw(10)<<position.y()
	<<setw(10)<<position.z()<<endl
	<<" radius of point          ="<< setw(10)<< rP<<endl
	<<" layer                    ="<< setw(10)<< irlayer<<endl<<endl;
    
    cout<<endl<<"test point in dead space"<<endl;
    xP = 0.355; yP= 0.; zP =.20;  
    position = HepPoint3D(xP,yP,zP);  rP = position.perp();
    irlayer = iADRSenseWireStore->layerNumber(layer_success,position);
    if (layer_success)
      { 
	cout<<endl<<"layerNumber success"<<endl;  
      }
    else	 
      { cout<<endl<<"layerNumber fails"<<endl; 
      }
    cout<<"DRGeom Processor layerNumber(position) ="<<irlayer<<endl;
    cout<<" sense wire layer for point"
	<< setw(10)<<position.x()
	<<setw(10)<<position.y()
	<<setw(10)<<position.z()<<endl
	<<" radius of point          ="<< setw(10)<< rP<<endl
	<<" layer                    ="<< setw(10)<< irlayer<<endl<<endl;
    
    cout<<endl<<"test point too far in z"<<endl;
    xP = 0.130 ; yP= 0.; zP =.40;  
    position = HepPoint3D(xP,yP,zP);  rP = position.perp();
    irlayer = iADRSenseWireStore->layerNumber(layer_success,position);
    if (layer_success)
      { 
	cout<<"layerNumber success"<<endl;  
      }
    else	 
      { cout<<"layerNumber fails"<<endl; 
      }
    cout<<"DRGeom Processor layerNumber(position) ="<<irlayer<<endl;
    cout<<" sense wire layer for point"
	<< setw(10)<<position.x()
	<< setw(10)<<position.y()
	<< setw(10)<<position.z()<<endl
	<<" radius of point          ="<< setw(10)<< rP<<endl
	<<" layer                    ="<< setw(10)<< irlayer<<endl<<endl;
    
// find wire closest to point HepPoint3D(xP,yP,zP)
// make a point at lyr, wire, zmin/2; 
// recover layer, closest wire
    unsigned int mylayer = 47;
    unsigned int mywire  = 164;
    double mydist = 0.0005;
    double myxdif = mydist/sqrt(2.);
    double myydif = mydist/sqrt(2.);
    const AEWireLayerCylindrical& MyLayer
      = iADRSenseWireStore->layer(mylayer);
    double myZ = MyLayer.zMin()/2.;
    cout<<"=================================="<<endl
	<<"Find closest wire" <<endl
	<<" Input lyr,wire = "<<mylayer<<","<<mywire<<endl
	<<" z              = "<<myZ<<endl
	<<" distance       = "<<mydist<<endl
	<<"=================================="<<endl;
    const AEWireInfo& MyWire = MyLayer.wire(mywire);
    AEWireInfo::WireDirZ MywDZ = MyWire.zPoint(myZ); 
    xP = MywDZ.first.x()+myxdif;
    yP = MywDZ.first.y()+myydif;
    zP = myZ;
    position = HepPoint3D(xP,yP,zP);
    double phiP = atan2(yP,xP);
    //find layer
    irlayer = iADRSenseWireStore->layerNumber(layer_success,position);
    if(!layer_success)
      {
	cout<<endl<<"Failed to find layer, wire for 3D point"<<endl
	    <<setw(10)<<xP<<setw(10)<<yP<<setw(10)<<zP<<endl
	    <<"phiP ="<<phiP<<"="<<180.*phiP/M_PI<<" deg "
	    <<"layer="<<irlayer<<endl;
      }
    else 
      {
	cout<<endl<<"Find layer, wire for 3D point"<<endl
	    <<setw(10)<<xP<<setw(10)<<yP<<setw(10)<<zP<<endl
	    <<"phiP ="<<phiP<<"="<<180.*phiP/M_PI<<" deg "
	    <<"layer="<<irlayer<<endl;
	const AEWireLayerCylindrical& Rlayer
	  = iADRSenseWireStore->layer(irlayer);
	// do search in Processor======================================
	/*
	 cout<<"search in DRGeomProcessor"<<endl;
	 const AEWireInfo& itwire0=Rlayer.wire(0);
	 AEWireInfo::WireDirZ wDZ0 = itwire0.zPoint(position.z());
	 phi0    = atan2(wDZ0.first.y(),wDZ0.first.x());
	 cout<<"phi0 at z"<<phi0<<"="<<180.*phiP/M_PI<<" deg "<<endl;
	 if(phiP < 0) { phiP = phiP + 2.0*M_PI;}    
	 if(phi0 < 0) { phi0 = phi0 + 2.0*M_PI;}
	 cout<<"phi0 at z"<<phi0<<"="<<180.*phiP/M_PI<<" deg "<<endl;
	 double phiDif = phiP-phi0;
	 if (phiDif < 0) { phiDif = phiDif + 2.0*M_PI;}  
	 int iwireM     = phiDif/Rlayer.cellPhiWidth();
	 int iwireMP    = iwireM + 1;
	 int iwireMM    = iwireM - 1;
	 cout<<"cell width="<<Rlayer.cellPhiWidth()
	 <<"="<<Rlayer.cellPhiWidth()*180./M_PI<<" deg "<<endl
	 <<"wire iwireM= "<< iwireM<<endl;
	 if(iwireMP==Rlayer.numberOfWires() ) { iwireMP    = 0;	 }
	 if(iwireMM<0 )  { iwireMM    = Rlayer.numberOfWires()-1; }
	 */       
	      // check distance to wires iwireM, iwireM+1, iwireM-1
       /*
	 const AEWireInfo& itwireM =Rlayer.wire(iwireM);
	 const AEWireInfo& itwireMP=Rlayer.wire(iwireMP);
	 const AEWireInfo& itwireMM=Rlayer.wire(iwireMM);
	 
	 double distanceM  = itwireM.pointDistance(position);
	 double distanceMP = itwireMP.pointDistance(position);
	 double distanceMM = itwireMM.pointDistance(position);
	 cout<<"layer "<< setw(4)<<irlayer<<" distance"<<endl
	 <<"wire  "<< setw(4)<<iwireM <<" "<< distanceM<<endl
	 <<"wire  "<< setw(4)<<iwireMP <<" "<< distanceMP<<endl
	 <<"wire  "<< setw(4)<<iwireMM <<" "<< distanceMM<<endl<<endl;
	 */
// ====================================================================

	cout<<endl<<"search in ADRSenseWireStore::wireObjectClosest"<<endl;
	DABoolean success ;
	const AEWireInfo& wireOS= 
	  iADRSenseWireStore->wireObjectClosest(success,position,irlayer); 
	
	if (success) 
	  { 
	    cout<<"wireObjectClosests gives wire with ends:"<<endl
		<<wireOS.point1()
		<<wireOS.point2()<<endl<<endl;
	  }
	else 
	  {
	    cout<<"NO Success"<<endl;
	  }
	
	cout<<endl<<"search in layer.wireNumberClosest"<<endl;
	unsigned int closestwire = 
	  Rlayer.wireNumberClosest(success,position);
	const AEWireInfo& wireC= Rlayer.wire(closestwire);
	if (success) 
	  {
	    cout<<"wireNumberClosest ="<<closestwire<<" ends:"<<endl
		<<wireC.point1()
		<<wireC.point2()<<endl<<endl;
	  }
	else 
	  { 
	    cout<<"NO Success"<<endl;
	  }
	
	cout<<endl<<"search in layer.wireObjectClosest"<<endl;
	const AEWireInfo& wireOC = Rlayer.wireObjectClosest(success,position);
	if (success) 
	  { 
	    cout<<"wireObjectClosests gives wire with ends:"<<endl
		<<wireOC.point1()
		<<wireOC.point2()<<endl<<endl;
	  }
	else 
	  {
	    cout<<"NO Success"<<endl;
	  }
      }
    
    //----------------------------+
    //Wire Properties Access      |
    //via ADRSenseWireStore       |
    //----------------------------+
    itLayer = 47 ;
    const AEWireLayerCylindrical& ilayer= iADRSenseWireStore->layer(itLayer);
    //**   const AEWireLayerRealistic& ilayer = iADRSenseWireStore->layer(itLayer);
    
    cout <<endl<<"Wire Access Example 1 wireMk/using wire() "<<endl
	 <<"==========================================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    outf <<endl<<"Wire Access Example 1 wireMk/using wire() "<<endl
	 <<"==========================================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    for( iwire=0;   iwire < 4; ++iwire )      
      {
	// use wire information from m_wires[]
	const AEWireInfo& itwire=iADRSenseWireStore->wire(itLayer,iwire);
	//	   const AEWireInfo& itwire=ilayer.wire(iwire);
	//recover wire #
	unsigned int wirenumber=ilayer.wireNumber(itwire);
	cout<<"iwire="<<iwire
	    <<" layer.wireNumber(itwire)="<<wirenumber<<endl;
      }
    
    for( iwire=0;   iwire < 4; ++iwire )      
      {
	// use wire information from m_wires[]
	const AEWireInfo& itwire=iADRSenseWireStore->wire(itLayer,iwire);
	//	  const AEWireInfo& itwire=ilayer.wire(iwire);
	// use wire information from m_wires[], w/copy constr
	//AEWireInfo itwire=ilayer.wire(iwire);
	//  make wire on stack 
	//AEWireInfo itwire = ilayer.wireMk(iwire);
	cout << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  itwire.point1() <<"  "<< itwire.point2() << endl;
	outf << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  itwire.point1() <<"  "<< itwire.point2() << endl;
	
      }

    double z0=0.0;
    double zE;
    double radius0;
    
    AEWireInfo::WireDirZ wDZ;
    cout <<endl<<"Wire Access Example 2 wireMk/using wire() "<<endl
	 <<"=========================================="<<endl
	 <<"Lyr   Wr         (x,y,z,r)              tangent vector   "<<endl;
    
    outf <<endl<<"Wire Access Example 2 wireMk/using wire() "<<endl
	 <<"=========================================="<<endl
	 <<"Lyr   Wr         (x,y,z,r)              tangent vector   "<<endl;
    for( iwire=0;   iwire < 2; ++iwire )      
      {
	// use wire information from m_wires[]
	//	 const AEWireInfo& itwire=ilayer.wire(iwire);
	// use wire information from m_wires[], w/copy constr
	const AEWireInfo& itwire=iADRSenseWireStore->wire(itLayer,iwire);
	 //	 const AEWireInfo itwire=ilayer.wire(iwire);
	 // make wire on stack
	 //AEWireInfo itwire = ilayer.wireMk(iwire);
	zE =  itwire.point1().z();
	wDZ = itwire.zPoint(zE);
	radius0 = wDZ.first.distance();
	phi0    = atan2(wDZ.first.y(),wDZ.first.x());
	cout << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  wDZ.first <<"  "<< wDZ.second  <<  endl
	     << "       z=zE: r(xyz),phi      " 
	     <<setw(15)<<radius0<<"  "<<setw(15)<<phi0<<endl;
	outf << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  wDZ.first <<"  "<< wDZ.second  <<  endl
	     << "       z=zE: r(xyz),phi      "  
	     <<setw(15)<<radius0<<"  "<<setw(15)<<phi0<<endl;

// test finding distance from wire
//=================================
	HepPoint3D position ;
	double xdif = .000100; 
	double ydif = .000100;     
	xP = wDZ.first.x() + xdif;
	yP = wDZ.first.y() + ydif;
	zP = wDZ.first.z();
	position = HepPoint3D(xP,yP,zP);
	double dist  = itwire.pointDistanceQ(position);
	double distL = itwire.pointDistance(position);
	cout<<" distance to sense wire layer for point"<<endl
	    << setw(10)<<position.x()<< setw(10)<<position.y()
	    << setw(10)<<position.z()<<endl
	    <<"expect sqrt(2)*.000100"<<endl
	    <<"linear approximation  "<< distL<<endl;	
	outf<<" distance to sense wire layer for point"<<endl
	    << setw(10)<<position.x()<< setw(10)<<position.y()
	    << setw(10)<<position.z()<<endl
	    <<"expect sqrt(2)*.000100"<<endl
	    <<"linear approximation  "<< distL<<endl;		
	if(dist>=0.)
	  {
	    cout<<"full calculation  =   "<<dist<<endl;
	    outf<<"full calculation  =   "<<dist<<endl;
	  }	
	else
	  {
	    cout<<"full calculation fails "<<dist<<endl;
	    outf<<"full calculation fails "<<dist<<endl;	     
	  }	
	
	wDZ = itwire.zPoint(z0);
	radius0 = wDZ.first.distance();
	phi0    = atan2(wDZ.first.y(),wDZ.first.x());
	cout << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  wDZ.first <<"  "<< wDZ.second  << endl
	     << "       z= 0: r(xyz),phi      " 
	     <<setw(15)<<radius0<<"  "<<setw(15)<<phi0<<endl;
	outf << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  wDZ.first <<"  "<< wDZ.second  << endl
	     << "       z= 0: r(xyz),phi      " 
	     <<setw(15)<<radius0<<"  "<<setw(15)<<phi0<<endl;
	
      }

    //chk ADRSenseWireStore::wireMk(address)======================================
    
    cout <<endl<<"Wire Access via ADRSenseWireStore::wireMk(WireAddress)"<<endl
	 <<"=====================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    
    outf <<endl<<"Wire Access via ADRSenseWireStore::wireMk(WireAddress)"<<endl
	 <<"=====================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    for( iwire=0;   iwire < 4; ++iwire )      
      { 
	const ADRWireAddress& iWireAddress = ADRWireAddress(itLayer,iwire);
	const AEWireInfo iswire =   iADRSenseWireStore->wireMk(iWireAddress);     
	cout << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  iswire.point1() <<"  "<< iswire.point2() << endl;
	outf << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  iswire.point1() <<"  "<< iswire.point2() << endl;
      }	

    outf <<endl<<"Wire Access via ADRSenseWireStore::wire(WireAddress)"<<endl
	 <<"=====================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    cout <<endl<<"Wire Access via ADRSenseWireStore::wire(WireAddress)"<<endl
	 <<"=====================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    
    //chk ADRSenseWireStore::wire=========================================
    for( iwire=0;   iwire < 4; ++iwire )      
      { 
	const ADRWireAddress& iWireAddress = ADRWireAddress(itLayer,iwire);
	const AEWireInfo iswire =   iADRSenseWireStore->wire(iWireAddress);     
	cout << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  iswire.point1() <<"  "<< iswire.point2() << endl;
	outf << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  iswire.point1() <<"  "<< iswire.point2() << endl;
      }
    //chk ADRSenseWireStore::wireMk(cell)=========================================
    
    cout <<endl<<"Wire Access via ADRSenseWireStore::wireMk(cellID)"<<endl
	 <<"=====================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    
    outf <<endl<<"Wire Access via ADRSenseWireStore::wireMk(cellID)"<<endl
	 <<"=====================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    for ( iwire=0;   iwire < 4; ++iwire ) 
      {
	CellID=iADRSenseWireStore->cellID(itLayer,iwire);
	const AEWireInfo iswire =   iADRSenseWireStore->wireMk(CellID);     
	cout << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  iswire.point1() <<"  "<< iswire.point2() << endl;
	outf << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  iswire.point1() <<"  "<< iswire.point2() << endl;
      }	
    //chk ADRSenseWireStore::wire(cell)=========================================
    
    cout <<endl<<"Wire Access via ADRSenseWireStore::wire(cellID)"<<endl
	 <<"=====================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    
    outf <<endl<<"Wire Access via ADRSenseWireStore::wire(cellID)"<<endl
	 <<"=====================" <<endl
	 <<"Lyr   Wr    East (x,y,z-)              West (x,y,z+)    z+(inch)" 
	 << endl ;
    for ( iwire=0;   iwire < 4; ++iwire ) 
      {
	CellID=iADRSenseWireStore->cellID(itLayer,iwire);
	const AEWireInfo iswire =   iADRSenseWireStore->wire(CellID);     
	cout << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  iswire.point1() <<"  "<< iswire.point2() << endl;
	outf << setw(3) << itLayer << setw(4)<< iwire <<" " 
	     <<  iswire.point1() <<"  "<< iswire.point2() << endl;
	if(iwire==3)
	  {	  
	    cout << setw(3) << itLayer << setw(4)<< iwire 
		 << " sagy= "
		 << setw(5)<<iswire.sagy()<<" m"<<endl;
	    outf << setw(3) << itLayer << setw(4)<< iwire 
		 << " sagy= "
		 << setw(5)<<iswire.sagy()<<" m"<<endl;
	  }
      }
    
    }
  else
    {
      report( ERROR, kFacilityString  ) 
	<< "Invalid ADRSenseWireStore pointer!" << endl ;
    }
  return ActionBase::kPassed;
}
//
//
// const member functions
//

//
// static member functions
//