UtilityApps/src/test_surfaces.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 //==========================================================================
0011
0012 // Framework include files
0013 #include "DD4hep/Detector.h"
0014
0015 #include "DDRec/Surface.h"
0016 #include "DDRec/DetectorSurfaces.h"
0017 #include "DDRec/SurfaceManager.h"
0018 #include "DDRec/SurfaceHelper.h"
0019 #include "DD4hep/DDTest.h"
0020
0021 #include "DD4hep/DD4hepUnits.h"
0022
0023 #include "lcio.h"
0024 #include "IO/LCReader.h"
0025 #include "EVENT/LCEvent.h"
0026 #include "EVENT/LCCollection.h"
0027 #include "EVENT/SimTrackerHit.h"
0028 #include "UTIL/ILDConf.h"
0029
0030 #include <map>
0031 #include <sstream>
0032
0033 using namespace std ;
0034 using namespace dd4hep ;
0035 using namespace dd4hep::detail;
0036 using namespace dd4hep::rec ;
0037
0038
0039 static DDTest test( "surfaces" ) ;
0040
0041 //=============================================================================
0042
0043 int main_wrapper(int argc, char** argv ){
0044
0045   if( argc < 3 ) {
0046     std::cout << " usage: test_surfaces compact.xml lcio_file.slcio" << std::endl ;
0047     exit(1) ;
0048   }
0049
0050   std::string inFile =  argv[1] ;
0051
0052   Detector& description = Detector::getInstance();
0053
0054   description.fromCompact( inFile );
0055
0056
0057 #if 0
0058   // create a list of all surfaces in the detector:
0059   DetElement world = description.world() ;
0060
0061   SurfaceHelper surfMan(  world ) ;
0062
0063   const SurfaceList& sL = surfMan.surfaceList() ;
0064
0065   // map of surfaces
0066   std::map< dd4hep::long64, Surface* > surfMap ;
0067
0068   for( SurfaceList::const_iterator it = sL.begin() ; it != sL.end() ; ++it ){
0069
0070     Surface* surf =  *it ;
0071
0072     // std::cout << " ------------------------- "
0073     //            << " surface: "  << *surf          << std::endl
0074     //            << " ------------------------- "  << std::endl ;
0075
0076
0077     surfMap[ surf->id() ] = surf ;
0078   }
0079 #else
0080
0081   SurfaceManager& surfMan = *description.extension< SurfaceManager >() ;
0082   const SurfaceMap& surfMap = *surfMan.map( "world" ) ;
0083
0084 #endif
0085
0086   //---------------------------------------------------------------------
0087   //    open lcio file with SimTrackerHits
0088   //---------------------------------------------------------------------
0089
0090   std::string lcioFileName = argv[2] ;
0091
0092   IO::LCReader* rdr = IOIMPL::LCFactory::getInstance()->createLCReader() ;
0093   rdr->open( lcioFileName ) ;
0094
0095   EVENT::LCEvent* evt = 0 ;
0096
0097
0098   while( ( evt = rdr->readNextEvent() ) != 0 ){
0099
0100     const std::vector< std::string >& colNames = *evt->getCollectionNames() ;
0101
0102     for(unsigned icol=0, ncol = colNames.size() ; icol < ncol ; ++icol ){
0103
0104       EVENT::LCCollection* col =  evt->getCollection( colNames[ icol ] ) ;
0105
0106       std::string typeName = col->getTypeName() ;
0107
0108       if( typeName != lcio::LCIO::SIMTRACKERHIT )
0109         continue ;
0110
0111       std::cout << "  -- testing collection : " <<  colNames[ icol ] << std::endl ;
0112
0113       std::string cellIDEcoding = col->getParameters().getStringVal("CellIDEncoding") ;
0114
0115       lcio::BitField64 idDecoder( cellIDEcoding ) ;
0116
0117       int nHit = col->getNumberOfElements() ;
0118
0119       for(int i=0 ; i< nHit ; ++i){
0120
0121         EVENT::SimTrackerHit* sHit = (EVENT::SimTrackerHit*) col->getElementAt(i) ;
0122
0123         dd4hep::FieldID id = sHit->getCellID0() ;
0124
0125         idDecoder.setValue( id ) ;
0126         //      std::cout << " simhit with cellid : " << idDecoder << std::endl ;
0127
0128 #if 0
0129         Surface* surf = surfMap[ id ] ;
0130 #else
0131         SurfaceMap::const_iterator si = surfMap.find( id )  ;
0132     //        Surface* surf = dynamic_cast<Surface*> ( ( si != surfMap.end()  ?  si->second  : 0 )   ) ;
0133         ISurface* surf = ( si != surfMap.end()  ?  si->second  : 0 )  ;
0134 #endif
0135
0136         std::stringstream sst ;
0137         sst << " surface found for id : " << std::hex << id  << std::dec  <<  "  "  << idDecoder.valueString() << std ::endl ;
0138
0139
0140         // ===== test that we have a surface with the correct ID for every hit ======================
0141
0142         test( surf != 0 , true , sst.str() ) ;
0143
0144
0145         if( surf != 0 ){
0146
0147           //  std::cout << " found surface " <<  *surf << std::endl ;
0148
0149           Vector3D point( sHit->getPosition()[0]* dd4hep::mm , sHit->getPosition()[1]* dd4hep::mm ,  sHit->getPosition()[2]* dd4hep::mm ) ;
0150
0151           double dist = surf->distance( point ) ;
0152
0153           bool isInside = surf->insideBounds( point )  ;
0154
0155
0156           sst.str("") ;
0157           sst << " point " << point << " is on surface " ;
0158
0159           // ====== test that hit points are inside their surface ================================
0160
0161           test( isInside , true , sst.str() ) ;
0162
0163           if( ! isInside ) {
0164
0165             std::cout << " found surface " <<  *surf << std::endl
0166                       << " id : " << idDecoder.valueString()
0167                       << " point : " << point
0168                       << " is inside : " <<  isInside
0169                       << " distance from surface : " << dist  << "  (units: cm) " << std::endl
0170                       << std::endl ;
0171           }
0172
0173           // ====== test that slightly moved hit points are inside their surface ================================
0174
0175       // this test does not make too much sense, depending on the position of the surface within the volume
0176
0177           // Vector3D point2 = point + 1e-5 * surf->normal() ;
0178           // sst.str("") ;
0179           // sst << " point2 " << point2 << " is on surface " ;
0180           // isInside = surf->insideBounds( point2 )  ;
0181           // test( isInside , true , sst.str() ) ;
0182
0183           // if( ! isInside ) {
0184
0185           //   std::cout << " found surface " <<  *surf << std::endl
0186           //             << " id : " << idDecoder.valueString()
0187           //             << " point : " << point
0188           //             << " is inside : " <<  isInside
0189           //             << " distance from surface : " << dist/dd4hep::mm << std::endl
0190           //             << std::endl ;
0191
0192           // }
0193
0194           // ====== test that moved hit points are outside their surface ================================
0195
0196           Vector3D point3 = point + 1e-3 * surf->normal() ;
0197           sst.str("") ;
0198           sst << " point3 " << point3 << " is not on surface " ;
0199           isInside = surf->insideBounds( point3)  ;
0200           test( isInside , false , sst.str() ) ;
0201
0202
0203
0204         } else {
0205
0206           std::cout << "ERROR:   no surface found for id: " << idDecoder << std::endl ;
0207         }
0208
0209       }
0210
0211     }
0212
0213
0214   }
0215   return 0;
0216 }
0217
0218 //=============================================================================
0219 #include "main.h"