UtilityApps/src/test_cellid_position_converter.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 // @author F.Gaede, DESY
0011 // @date May, 2017
0012 //==========================================================================
0013
0014 // Framework include files
0015 #include "DD4hep/Detector.h"
0016 #include "DD4hep/DDTest.h"
0017
0018 #include "DD4hep/DD4hepUnits.h"
0019 #include "DD4hep/BitFieldCoder.h"
0020 #include "DDRec/CellIDPositionConverter.h"
0021
0022 #include "lcio.h"
0023 #include "IO/LCReader.h"
0024 #include "EVENT/LCEvent.h"
0025 #include "EVENT/LCCollection.h"
0026 #include "EVENT/SimCalorimeterHit.h"
0027
0028 #include <sstream>
0029
0030 using namespace std;
0031 using namespace dd4hep;
0032 using namespace dd4hep::detail;
0033 using namespace dd4hep::rec;
0034 using namespace lcio;
0035
0036
0037 static DDTest test( "cellid_position_converter" ) ;
0038
0039 //=============================================================================
0040
0041 const double epsilon = dd4hep::micrometer ;
0042 const int maxHit = 100 ;
0043
0044
0045 double dist( const Position& p0, const Position& p1 ){
0046   Position p2 = p1 - p0 ;
0047   return p2.r() ;
0048 }
0049
0050
0051 struct TestCounter{
0052   unsigned passed{} ;
0053   unsigned failed{} ;
0054 } ;
0055
0056 struct TestCounters{
0057   TestCounter position{} ;
0058   TestCounter cellid{} ;
0059 };
0060
0061 typedef std::map<std::string, TestCounters > TestMap ;
0062
0063
0064
0065 int main_wrapper(int argc, char** argv ){
0066
0067   if( argc < 3 ) {
0068     std::cout << " usage: test_cellid_position_converter compact.xml lcio_file.slcio" << std::endl ;
0069     exit(1) ;
0070   }
0071
0072   std::string inFile =  argv[1] ;
0073
0074   Detector& description = Detector::getInstance();
0075
0076   description.fromCompact( inFile );
0077
0078   CellIDPositionConverter idposConv( description )  ;
0079
0080
0081   //---------------------------------------------------------------------
0082   //    open lcio file with SimCalorimeterHits
0083   //---------------------------------------------------------------------
0084
0085   std::string lcioFileName = argv[2] ;
0086
0087   LCReader* rdr = LCFactory::getInstance()->createLCReader() ;
0088   rdr->open( lcioFileName ) ;
0089
0090   LCEvent* evt = 0 ;
0091
0092
0093   // use only hits from these collections
0094   std::set< std::string > subset = {} ;
0095   //{"BeamCalCollection" } ; //{"EcalBarrelCollection" } ; //{"HcalBarrelRegCollection"} ;
0096
0097
0098   // ignore all hits from these collections
0099   std::set< std::string > subsetIgnore = {"HCalBarrelRPCHits","HCalECRingRPCHits","HCalEndcapRPCHits" } ;
0100
0101   TestMap tMap ;
0102
0103   while( ( evt = rdr->readNextEvent() ) != 0 ){
0104
0105     const std::vector< std::string >& colNames = *evt->getCollectionNames() ;
0106
0107     for(unsigned icol=0, ncol = colNames.size() ; icol < ncol ; ++icol ){
0108
0109       LCCollection* col =  evt->getCollection( colNames[ icol ] ) ;
0110
0111       std::string typeName = col->getTypeName() ;
0112
0113       if( typeName != lcio::LCIO::SIMCALORIMETERHIT )
0114         continue ;
0115
0116       if( !subset.empty() && subset.find( colNames[icol] ) ==  subset.end()  )
0117     continue ;
0118
0119       if( !subsetIgnore.empty() && subsetIgnore.find( colNames[icol] ) !=  subsetIgnore.end()  )
0120         continue ;
0121
0122       std::cout << "  -- testing collection : " <<  colNames[ icol ] << std::endl ;
0123
0124       std::string cellIDEcoding = col->getParameters().getStringVal("CellIDEncoding") ;
0125
0126       dd4hep::BitFieldCoder idDecoder0( cellIDEcoding ) ;
0127       dd4hep::BitFieldCoder idDecoder1( cellIDEcoding ) ;
0128
0129       int nHit = std::min( col->getNumberOfElements(), maxHit )  ;
0130
0131
0132       for(int i=0 ; i< nHit ; ++i){
0133
0134         SimCalorimeterHit* sHit = (SimCalorimeterHit*) col->getElementAt(i) ;
0135
0136         dd4hep::CellID id0 = sHit->getCellID0() ;
0137         dd4hep::CellID id1 = sHit->getCellID1() ;
0138
0139     dd4hep::CellID id =  idDecoder0.toLong( id0 , id1 ) ;
0140
0141     Position point( sHit->getPosition()[0]* dd4hep::mm , sHit->getPosition()[1]* dd4hep::mm ,  sHit->getPosition()[2]* dd4hep::mm ) ;
0142
0143
0144     // ====== test cellID to position and position to cellID conversion  ================================
0145     DetElement det = idposConv.findDetElement( point ) ;
0146
0147     CellID idFromDecoder = idposConv.cellID( point ) ;
0148
0149     std::stringstream sst ;
0150     sst << " compare ids: " << det.name() << " " <<  idDecoder0.valueString(id) << "  -  " << idDecoder1.valueString(idFromDecoder) ;
0151
0152     test( id, idFromDecoder,  sst.str() ) ;
0153
0154     if( ! strcmp( test.last_test_status() , "PASSED" ) )
0155       tMap[ colNames[icol] ].cellid.passed++ ;
0156     else
0157       tMap[ colNames[icol] ].cellid.failed++ ;
0158
0159     Position pointFromDecoder = idposConv.position( id ) ;
0160
0161     double d = dist(pointFromDecoder, point)  ;
0162     std::stringstream sst1 ;
0163     sst1 << " dist " << d << " ( " <<  point << " ) - ( " << pointFromDecoder << " )  - detElement: "
0164         << det.name() ;
0165
0166     test( d < epsilon , true  , sst1.str()  ) ;
0167
0168     if( ! strcmp( test.last_test_status() , "PASSED" ) )
0169       tMap[ colNames[icol] ].position.passed++ ;
0170     else
0171       tMap[ colNames[icol] ].position.failed++ ;
0172
0173       }
0174     }
0175
0176   }
0177
0178   // print summary
0179
0180   std::cout << "\n ----------------------- summary  ----------------------   " << std::endl ;
0181
0182
0183   for( const auto& res : tMap )  {
0184     const std::string& name = res.first;
0185     unsigned total      = res.second.position.passed+res.second.position.failed ;
0186     unsigned pos_failed = res.second.position.failed ;
0187     unsigned id_failed  = res.second.cellid.failed ;
0188
0189
0190     printf(" %-30s \t  failed position: %5d  failed cellID:  %5d    of total: %5d   \n",
0191            name.c_str(), pos_failed , id_failed, total ) ;
0192
0193   }
0194   std::cout << "\n -------------------------------------------------------- " << std::endl ;
0195
0196
0197   return 0;
0198 }
0199
0200 //=============================================================================
0201 #include "main.h"