diff --git a/Reconstruction/SiliconTracking/CMakeLists.txt b/Reconstruction/SiliconTracking/CMakeLists.txt index 27533d0dd47b0353ce1cdc2b085177459dd57f91..025f3449e69af5d0552c1b78a0fc7238f3b8910b 100644 --- a/Reconstruction/SiliconTracking/CMakeLists.txt +++ b/Reconstruction/SiliconTracking/CMakeLists.txt @@ -4,6 +4,7 @@ find_package(GEAR REQUIRED) find_package(GSL REQUIRED ) find_package(LCIO REQUIRED ) find_package(EDM4HEP REQUIRED ) +find_package(DD4hep COMPONENTS DDCore DDRec REQUIRED) gaudi_depends_on_subdirs( Service/GearSvc diff --git a/Service/TrackSystemSvc/CMakeLists.txt b/Service/TrackSystemSvc/CMakeLists.txt index 7e68d7615c49ff58520b09213cab1b83c74b6349..c44d75470fe19657c4ac72eeb622bc8f08401793 100644 --- a/Service/TrackSystemSvc/CMakeLists.txt +++ b/Service/TrackSystemSvc/CMakeLists.txt @@ -6,6 +6,7 @@ find_package(LCIO REQUIRED) find_package(EDM4HEP REQUIRED) #find_package(KalTest REQUIRED) #find_package(KalDet REQUIRED) +find_package(DD4hep COMPONENTS DDCore DDRec REQUIRED) gaudi_depends_on_subdirs(Service/GearSvc Detector/DetInterface Utilities/DataHelper Utilities/KalTest Utilities/KalDet) diff --git a/Service/TrackSystemSvc/src/MarlinKalTest.cc b/Service/TrackSystemSvc/src/MarlinKalTest.cc index 38a4f5fc5c76b1c85827b185ee12b41ba6fe9849..0c0dde8c5a14c4663cce7ae2af288447d375bbae 100644 --- a/Service/TrackSystemSvc/src/MarlinKalTest.cc +++ b/Service/TrackSystemSvc/src/MarlinKalTest.cc @@ -128,7 +128,7 @@ namespace MarlinTrk{ bool SIT_found = false ; try{ - ILDSITKalDetector* sitdet = new ILDSITKalDetector( *_gearMgr ) ; + ILDSITKalDetector* sitdet = new ILDSITKalDetector( *_gearMgr, _geoSvc ) ; // store the measurement layer id's for the active layers this->storeActiveMeasurementModuleIDs(sitdet); _det->Install( *sitdet ) ; @@ -151,7 +151,7 @@ namespace MarlinTrk{ } try{ - ILDSETKalDetector* setdet = new ILDSETKalDetector( *_gearMgr ) ; + ILDSETKalDetector* setdet = new ILDSETKalDetector( *_gearMgr, _geoSvc ) ; // store the measurement layer id's for the active layers this->storeActiveMeasurementModuleIDs(setdet); _det->Install( *setdet ) ; @@ -163,7 +163,7 @@ namespace MarlinTrk{ bool FTD_found = false ; try{ - ILDFTDKalDetector* ftddet = new ILDFTDKalDetector( *_gearMgr ) ; + ILDFTDKalDetector* ftddet = new ILDFTDKalDetector( *_gearMgr, _geoSvc ) ; // store the measurement layer id's for the active layers this->storeActiveMeasurementModuleIDs(ftddet); _det->Install( *ftddet ) ; @@ -186,7 +186,7 @@ namespace MarlinTrk{ } try{ - ILDTPCKalDetector* tpcdet = new ILDTPCKalDetector( *_gearMgr ) ; + ILDTPCKalDetector* tpcdet = new ILDTPCKalDetector( *_gearMgr, _geoSvc ) ; // store the measurement layer id's for the active layers this->storeActiveMeasurementModuleIDs(tpcdet); _det->Install( *tpcdet ) ; diff --git a/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc b/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc index 6a4ec5925048849ee2048de91c2a4cf518ae1805..9f83748cb7a7df9fa1477ca20df7bbb1bc11a7cc 100644 --- a/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc +++ b/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc @@ -122,7 +122,7 @@ namespace MarlinTrk { return this->addHit( trkhit, ml->ConvertLCIOTrkHit(trkhit), ml) ; } else { - std::cout << "MarlinKalTestTrack::addHit: trkhit = " << trkhit.id() << " addr: " << trkhit << " ml = " << ml << std::endl ; + //std::cout << "MarlinKalTestTrack::addHit: trkhit = " << trkhit.id() << " addr: " << trkhit << " ml = " << ml << std::endl ; //streamlog_out( ERROR ) << " MarlinKalTestTrack::addHit - bad inputs " << trkhit << " ml : " << ml << std::endl ; return bad_intputs ; } @@ -609,6 +609,7 @@ namespace MarlinTrk { // get the measurement layer of the current hit const ILDVMeasLayer* ml = dynamic_cast<const ILDVMeasLayer*>( &(kalhit->GetMeasLayer() ) ) ; TVector3 pos = ml->HitToXv(*kalhit); + /* std::cout << "debug: Kaltrack::addAndFit : site discarded! at index : " << ml->GetIndex() << " for type " << ml->GetName() << " chi2increment = " << chi2increment @@ -623,8 +624,8 @@ namespace MarlinTrk { << dynamic_cast<const ILDVMeasLayer*>( &(kalhit->GetMeasLayer() ) )->getCellIDs()[i] << std::endl ; } + */ - #ifdef MARLINTRK_DIAGNOSTICS_ON _ktest->_diagnostics.record_rejected_site(kalhit, temp_site); #endif diff --git a/Service/TrackSystemSvc/src/TrackSystemSvc.cpp b/Service/TrackSystemSvc/src/TrackSystemSvc.cpp index 78bb2779324d89edf505492d67e566667dcf414e..3f1abc0d3489f5a01c883afd613131dcf61f4433 100644 --- a/Service/TrackSystemSvc/src/TrackSystemSvc.cpp +++ b/Service/TrackSystemSvc/src/TrackSystemSvc.cpp @@ -17,18 +17,24 @@ TrackSystemSvc::~TrackSystemSvc(){ MarlinTrk::IMarlinTrkSystem* TrackSystemSvc::getTrackSystem(){ if(!m_trackSystem){ + gear::GearMgr* mgr=0; auto _gear = service<IGearSvc>("GearSvc"); if ( !_gear ) { - error() << "Failed to find GearSvc ..." << endmsg; - return 0; + info() << "Failed to find GearSvc ..." << endmsg; + } + else{ + mgr = _gear->getGearMgr(); } - gear::GearMgr* mgr = _gear->getGearMgr(); auto _geoSvc = service<IGeoSvc>("GeoSvc"); if ( !_geoSvc ) { - error() << "Failed to find GeoSvc ..." << endmsg; + info() << "Failed to find GeoSvc ..." << endmsg; + } + if(mgr==0&&_geoSvc==0){ + fatal() << "Both GearSvc and GeoSvc invalid!" << endmsg; return 0; } + m_trackSystem = new MarlinTrk::MarlinKalTest( *mgr, _geoSvc ) ; } return m_trackSystem; diff --git a/Utilities/KalDet/CMakeLists.txt b/Utilities/KalDet/CMakeLists.txt index d71f37dd3168dc92478f94eb5ac891f1cb418acd..d95a5cc27cffd698cea0cb37100e81e99e086e1b 100644 --- a/Utilities/KalDet/CMakeLists.txt +++ b/Utilities/KalDet/CMakeLists.txt @@ -9,6 +9,7 @@ find_package(LCIO) find_package(GEAR) find_package(ROOT COMPONENTS MathCore) find_package(EDM4HEP REQUIRED) +find_package(DD4hep COMPONENTS DDCore DDRec REQUIRED) gaudi_depends_on_subdirs( Detector/DetInterface @@ -75,5 +76,5 @@ set( KalDetLib_srcs ${LIB_SOURCES} ${COMMON_SOURCES} ) gaudi_add_library(KalDetLib ${KalDetLib_srcs} PUBLIC_HEADERS kaldet - LINK_LIBRARIES GaudiKernel ROOT CLHEP LCIO ${GEAR_LIBRARIES} KalTestLib EDM4HEP::edm4hep EDM4HEP::edm4hepDict + LINK_LIBRARIES GaudiKernel ROOT CLHEP LCIO ${GEAR_LIBRARIES} KalTestLib EDM4HEP::edm4hep EDM4HEP::edm4hepDict ${DD4hep_COMPONENT_LIBRARIES} ) diff --git a/Utilities/KalDet/kaldet/ILDFTDKalDetector.h b/Utilities/KalDet/kaldet/ILDFTDKalDetector.h index 36cf54b8e5a889a09431635740c3e2c3035d6d6a..80d348a8a58a6b916d8e7f901099a3129d5bb1e6 100644 --- a/Utilities/KalDet/kaldet/ILDFTDKalDetector.h +++ b/Utilities/KalDet/kaldet/ILDFTDKalDetector.h @@ -11,6 +11,7 @@ class TNode; class TVector3; +class IGeoSvc; namespace gear{ class GearMgr ; @@ -21,7 +22,7 @@ class ILDFTDKalDetector : public TVKalDetector { public: /** Initialize the FTD from GEAR */ - ILDFTDKalDetector( const gear::GearMgr& gearMgr ); + ILDFTDKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc ); private: diff --git a/Utilities/KalDet/kaldet/ILDSETKalDetector.h b/Utilities/KalDet/kaldet/ILDSETKalDetector.h index 9327aa8a065e8ff6243df53193df544e37389381..6c434cba9d28ee3c89c9d5dbec35dc442eef8789 100644 --- a/Utilities/KalDet/kaldet/ILDSETKalDetector.h +++ b/Utilities/KalDet/kaldet/ILDSETKalDetector.h @@ -16,18 +16,20 @@ namespace gear{ class GearMgr ; } +class IGeoSvc; class ILDSETKalDetector : public TVKalDetector { public: /** Initialize the SET from GEAR */ - ILDSETKalDetector( const gear::GearMgr& gearMgr ); + ILDSETKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc=0 ); private: void setupGearGeom( const gear::GearMgr& gearMgr ) ; + void setupGearGeom( IGeoSvc* geoSvc ); int _nLayers ; double _bZ ; diff --git a/Utilities/KalDet/kaldet/ILDSITKalDetector.h b/Utilities/KalDet/kaldet/ILDSITKalDetector.h index 639e1b27bb76512ad05c8570814b218a9480cb5b..d9eb3059b680eb16823cef7de3b315659c056272 100644 --- a/Utilities/KalDet/kaldet/ILDSITKalDetector.h +++ b/Utilities/KalDet/kaldet/ILDSITKalDetector.h @@ -16,18 +16,20 @@ namespace gear{ class GearMgr ; } +class IGeoSvc; class ILDSITKalDetector : public TVKalDetector { public: /** Initialize the SIT from GEAR */ - ILDSITKalDetector( const gear::GearMgr& gearMgr ); + ILDSITKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc ); private: void setupGearGeom( const gear::GearMgr& gearMgr ) ; + void setupGearGeom( IGeoSvc* geoSvc ); int _nLayers ; double _bZ ; diff --git a/Utilities/KalDet/kaldet/ILDTPCKalDetector.h b/Utilities/KalDet/kaldet/ILDTPCKalDetector.h index e5ab8db97c723357f2fbd8e964ab319ec7102ca6..a249f22ce1a7133dc88135ec233c785b54d777b7 100644 --- a/Utilities/KalDet/kaldet/ILDTPCKalDetector.h +++ b/Utilities/KalDet/kaldet/ILDTPCKalDetector.h @@ -14,12 +14,13 @@ namespace gear{ class GearMgr ; } +class IGeoSvc; class ILDTPCKalDetector : public TVKalDetector { public: /** Initialize the TPC from GEAR */ - ILDTPCKalDetector( const gear::GearMgr& gearMgr ); + ILDTPCKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc=0 ); private: diff --git a/Utilities/KalDet/src/ild/ftd/ILDFTDKalDetector.cc b/Utilities/KalDet/src/ild/ftd/ILDFTDKalDetector.cc index 3ac020812669ae175d17c9f9df6c983ce506790e..799bddf94354086f1054107522e9df3e82722f21 100644 --- a/Utilities/KalDet/src/ild/ftd/ILDFTDKalDetector.cc +++ b/Utilities/KalDet/src/ild/ftd/ILDFTDKalDetector.cc @@ -5,6 +5,10 @@ #include <sstream> +#include "DetInterface/IGeoSvc.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" + #include "gear/GEAR.h" #include "gear/BField.h" #include "gearimpl/Util.h" @@ -19,18 +23,23 @@ #include <UTIL/ILDConf.h> // #include "streamlog/streamlog.h" - +#include "CLHEP/Units/SystemOfUnits.h" #include "TVector3.h" -ILDFTDKalDetector::ILDFTDKalDetector( const gear::GearMgr& gearMgr ) : +ILDFTDKalDetector::ILDFTDKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc ) : TVKalDetector(300), _nDisks(0) // SJA:FIXME initial size, 300 looks reasonable for ILD, though this would be better stored as a const somewhere { // streamlog_out(DEBUG1) << "ILDFTDKalDetector building FTD detector using GEAR " << std::endl ; - MaterialDataBase::Instance().registerForService(gearMgr); - setupGearGeom( gearMgr ) ; - + MaterialDataBase::Instance().registerForService(gearMgr, geoSvc); + if(geoSvc){ + setupGearGeom( geoSvc ) ; + } + else{ + setupGearGeom( gearMgr ); + } + this->build_staggered_design(); @@ -380,7 +389,7 @@ void ILDFTDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ } - + //std::cout << "=============FTD strip angle: " << strip_angle_deg << "==============" << std::endl; _nDisks = ftdlayers.getNLayers() ; // just do the first disk for now _FTDgeo.resize(_nDisks); @@ -419,6 +428,11 @@ void ILDFTDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ _FTDgeo[disk].stripAngle = 0.0 ; } + + //std::cout << _FTDgeo[disk].nPetals << " " << _FTDgeo[disk].dphi << " " << _FTDgeo[disk].phi0 << " " << _FTDgeo[disk].alpha << " " + // << _FTDgeo[disk].rInner << " " << _FTDgeo[disk].height << " " << _FTDgeo[disk].innerBaseLength << " " << _FTDgeo[disk].outerBaseLength << " " + // << _FTDgeo[disk].senThickness << " " << _FTDgeo[disk].supThickness << " " << _FTDgeo[disk].senZPos_even_front << " " << _FTDgeo[disk].senZPos_odd_front << " " + // << _FTDgeo[disk].isDoubleSided << " " << _FTDgeo[disk].isStripReadout << " " << _FTDgeo[disk].nSensors << " " << _FTDgeo[disk].stripAngle << std::endl; //////////////////////////////////////////////////////////////////////////////////////////////////////// // Assertions //////////////////////////////////////////////////////////////////////////////////// @@ -499,3 +513,116 @@ void ILDFTDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ } + +void ILDFTDKalDetector::setupGearGeom( IGeoSvc* geoSvc ){ + dd4hep::DetElement world = geoSvc->getDD4HepGeo(); + dd4hep::DetElement ftd; + const std::map<std::string, dd4hep::DetElement>& subs = world.children(); + for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){ + if(it->first!="FTD") continue; + ftd = it->second; + } + dd4hep::rec::ZDiskPetalsData* ftdData = nullptr; + try{ + ftdData = ftd.extension<dd4hep::rec::ZDiskPetalsData>(); + } + catch(std::runtime_error& e){ + std::cout << e.what() << " " << ftdData << std::endl; + throw GaudiException(e.what(), "FATAL", StatusCode::FAILURE); + } + + const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + _bZ = field.z(); + + double strip_angle_deg = ftdData->angleStrip/CLHEP::degree; + bool strip_angle_present = true; + + std::vector<dd4hep::rec::ZDiskPetalsData::LayerLayout>& ftdlayers = ftdData->layers; + _nDisks = ftdlayers.size() ; + + _FTDgeo.resize(_nDisks); + + double eps = 1.0e-08; + //std::cout << "=============FTD strip angle: " << strip_angle_deg << "==============" << std::endl; + for(int disk=0; disk< _nDisks; ++disk){ + dd4hep::rec::ZDiskPetalsData::LayerLayout& ftdlayer = ftdlayers[disk]; + _FTDgeo[disk].nPetals = ftdlayer.petalNumber ; + _FTDgeo[disk].dphi = 2*M_PI / _FTDgeo[disk].nPetals ; + _FTDgeo[disk].phi0 = ftdlayer.phi0 ; + _FTDgeo[disk].alpha = ftdlayer.alphaPetal ; + _FTDgeo[disk].rInner = ftdlayer.distanceSensitive*CLHEP::cm ; + _FTDgeo[disk].height = ftdlayer.lengthSensitive*CLHEP::cm ; + _FTDgeo[disk].innerBaseLength = ftdlayer.widthInnerSensitive*CLHEP::cm ; + _FTDgeo[disk].outerBaseLength = ftdlayer.widthOuterSensitive*CLHEP::cm ; + _FTDgeo[disk].senThickness = ftdlayer.thicknessSensitive*CLHEP::cm ; + _FTDgeo[disk].supThickness = ftdlayer.thicknessSupport*CLHEP::cm ; + + _FTDgeo[disk].senZPos_even_front = ftdlayer.zPosition*CLHEP::cm - ftdlayer.zOffsetSensitive*CLHEP::cm;//getSensitiveZposition(disk, 0, 1) ; + _FTDgeo[disk].senZPos_odd_front = ftdlayer.zPosition*CLHEP::cm - ftdlayer.zOffsetSensitive*CLHEP::cm - 2*ftdlayer.zOffsetSupport*CLHEP::cm;//getSensitiveZposition(disk, 1, 1) ; + + _FTDgeo[disk].isDoubleSided = ftdlayer.typeFlags[dd4hep::rec::ZDiskPetalsData::SensorType::DoubleSided]; + _FTDgeo[disk].isStripReadout = !((bool)ftdlayer.typeFlags[dd4hep::rec::ZDiskPetalsData::SensorType::Pixel]); + _FTDgeo[disk].nSensors = ftdlayer.sensorsPerPetal; + + + if (strip_angle_present) { + _FTDgeo[disk].stripAngle = strip_angle_deg * M_PI/180 ; + } else { + _FTDgeo[disk].stripAngle = 0.0 ; + } + //std::cout << _FTDgeo[disk].nPetals << " " << _FTDgeo[disk].dphi << " " << _FTDgeo[disk].phi0 << " " << _FTDgeo[disk].alpha << " " + // << _FTDgeo[disk].rInner << " " << _FTDgeo[disk].height << " " << _FTDgeo[disk].innerBaseLength << " " << _FTDgeo[disk].outerBaseLength << " " + // << _FTDgeo[disk].senThickness << " " << _FTDgeo[disk].supThickness << " " << _FTDgeo[disk].senZPos_even_front << " " << _FTDgeo[disk].senZPos_odd_front << " " + // << _FTDgeo[disk].isDoubleSided << " " << _FTDgeo[disk].isStripReadout << " " << _FTDgeo[disk].nSensors << " " << _FTDgeo[disk].stripAngle << std::endl; + + assert( _FTDgeo[disk].nPetals%2 == 0 ); + assert( fabs( ftdlayer.widthInnerSupport - ftdlayer.widthInnerSensitive ) < eps ); + assert( fabs( ftdlayer.widthOuterSupport - ftdlayer.widthOuterSensitive ) < eps ); + assert( fabs( ftdlayer.lengthSupport - ftdlayer.lengthSensitive ) < eps ); + assert( fabs( ftdlayer.distanceSupport - ftdlayer.distanceSensitive ) < eps ); + if( _FTDgeo[disk].isDoubleSided ) assert( _FTDgeo[disk].nSensors%2 == 0 ); + assert( fabs( ftdlayer.alphaPetal ) < eps ); + /* + if( _FTDgeo[disk].isDoubleSided ){ + + for( int iPetal=0; iPetal< _FTDgeo[disk].nPetals; iPetal++){ + + int sensors1Side = _FTDgeo[disk].nSensors/2; + for( int iSensor=2; iSensor <= sensors1Side; iSensor++ ){ + + assert( fabs( ftdlayers.getSensitiveZposition( disk, iPetal, iSensor) - ftdlayers.getSensitiveZposition( disk, iPetal, iSensor-1) ) < eps ); + + } + + for( int iSensor=sensors1Side + 2; iSensor <= _FTDgeo[disk].nSensors; iSensor++ ){ + + assert( fabs( ftdlayers.getSensitiveZposition( disk, iPetal, iSensor) - ftdlayers.getSensitiveZposition( disk, iPetal, iSensor-1) ) < eps ); + assert( fabs( ftdlayers.getSensitiveZposition( disk, iPetal, iSensor) ) - fabs( ftdlayers.getSensitiveZposition( disk, iPetal, iSensor-sensors1Side ) ) >= 0. ); + } + } + } + else{ + + for( int iPetal=0; iPetal< _FTDgeo[disk].nPetals; iPetal++){ + + for(int iSensor=2; iSensor <= _FTDgeo[disk].nSensors; iSensor++ ){ + + assert( fabs( ftdlayers.getSensitiveZposition( disk, iPetal, iSensor) - ftdlayers.getSensitiveZposition( disk, iPetal, iSensor-1) ) < eps ); + + } + } + } + + for( int iPetal=0; iPetal< _FTDgeo[disk].nPetals; iPetal++){ + + int side = ftdlayers.getSupportZposition( disk, iPetal ) > 0? 1 : -1; + + double endSensitive = ftdlayers.getSensitiveZposition( disk, iPetal, 1) + side * _FTDgeo[disk].senThickness/2.; + double endSupport = ftdlayers.getSupportZposition( disk, iPetal ) -side * _FTDgeo[disk].supThickness/2.; + + assert( fabs( endSensitive- endSupport ) < eps ); + + } + */ + } +} diff --git a/Utilities/KalDet/src/ild/ftd/ILDFTDKalDetector.h b/Utilities/KalDet/src/ild/ftd/ILDFTDKalDetector.h index 36cf54b8e5a889a09431635740c3e2c3035d6d6a..deee9061226776bfd3986fb2cc6c21b65ab469e9 100644 --- a/Utilities/KalDet/src/ild/ftd/ILDFTDKalDetector.h +++ b/Utilities/KalDet/src/ild/ftd/ILDFTDKalDetector.h @@ -11,17 +11,17 @@ class TNode; class TVector3; +class IGeoSvc; namespace gear{ class GearMgr ; } - class ILDFTDKalDetector : public TVKalDetector { public: /** Initialize the FTD from GEAR */ - ILDFTDKalDetector( const gear::GearMgr& gearMgr ); + ILDFTDKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc ); private: @@ -80,6 +80,7 @@ private: void setupGearGeom( const gear::GearMgr& gearMgr ) ; + void setupGearGeom( IGeoSvc* geoSvc ); int _nDisks ; double _bZ ; diff --git a/Utilities/KalDet/src/ild/set/ILDSETKalDetector.cc b/Utilities/KalDet/src/ild/set/ILDSETKalDetector.cc index a8ae3eda05e97f1c3b659f9213e505a67f1db102..f1f7c6751f3a6ff1f76719298c9d5f34a05cd717 100644 --- a/Utilities/KalDet/src/ild/set/ILDSETKalDetector.cc +++ b/Utilities/KalDet/src/ild/set/ILDSETKalDetector.cc @@ -9,12 +9,17 @@ #include <UTIL/BitField64.h> #include <UTIL/ILDConf.h> +#include "DetInterface/IGeoSvc.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" + #include <gear/GEAR.h> #include "gear/BField.h" #include <gear/ZPlanarParameters.h> #include <gear/ZPlanarLayerLayout.h> #include "gearimpl/Util.h" +#include "CLHEP/Units/SystemOfUnits.h" #include "TMath.h" #include "math.h" @@ -22,7 +27,7 @@ // #include "streamlog/streamlog.h" -ILDSETKalDetector::ILDSETKalDetector( const gear::GearMgr& gearMgr ) +ILDSETKalDetector::ILDSETKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc ) : TVKalDetector(300) // SJA:FIXME initial size, 300 looks reasonable for ILD, though this would be better stored as a const somewhere { @@ -35,7 +40,12 @@ ILDSETKalDetector::ILDSETKalDetector( const gear::GearMgr& gearMgr ) TMaterial & silicon = *MaterialDataBase::Instance().getMaterial("silicon"); TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("carbon"); - this->setupGearGeom(gearMgr) ; + if(geoSvc){ + setupGearGeom( geoSvc ) ; + } + else{ + setupGearGeom( gearMgr ); + } if (_isStripDetector) { // streamlog_out(DEBUG4) << "\t\t building SET detector as STRIP Detector." << std::endl ; @@ -241,7 +251,7 @@ void ILDSETKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ // if this is not done then the exposed areas of the support would leave a carbon - air boundary, // which if traversed in the reverse direction to the next boundary then the track would be propagated through carbon // for a significant distance - + //std::cout << "=============SET strip angle: " << strip_angle_deg << "==============" << std::endl; for( int layer=0; layer < _nLayers; ++layer){ @@ -271,7 +281,9 @@ void ILDSETKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ } else { _SETgeo[layer].stripAngle = 0.0 ; } - + //std::cout << _SETgeo[layer].nLadders << " " << _SETgeo[layer].phi0 << " "<< _SETgeo[layer].dphi << " " << _SETgeo[layer].senRMin << " " << _SETgeo[layer].supRMin << " " + // << _SETgeo[layer].length << " " << _SETgeo[layer].width << " " << _SETgeo[layer].offset << " " << _SETgeo[layer].senThickness << " " << _SETgeo[layer].supThickness << " " + // << _SETgeo[layer].nSensorsPerLadder << " " << _SETgeo[layer].sensorLength << std::endl; // streamlog_out(DEBUG0) << " layer = " << layer << std::endl; // streamlog_out(DEBUG0) << " nSensorsPerLadder = " << _SETgeo[layer].nSensorsPerLadder << std::endl; // streamlog_out(DEBUG0) << " sensorLength = " << _SETgeo[layer].sensorLength << std::endl; @@ -279,8 +291,57 @@ void ILDSETKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ // streamlog_out(DEBUG0) << " _isStripDetector = " << _isStripDetector << std::endl; } - - - - +} + +void ILDSETKalDetector::setupGearGeom( IGeoSvc* geoSvc ){ + dd4hep::DetElement world = geoSvc->getDD4HepGeo(); + dd4hep::DetElement set; + const std::map<std::string, dd4hep::DetElement>& subs = world.children(); + for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){ + if(it->first!="SET") continue; + set = it->second; + } + dd4hep::rec::ZPlanarData* setData = nullptr; + try{ + setData = set.extension<dd4hep::rec::ZPlanarData>(); + } + catch(std::runtime_error& e){ + std::cout << e.what() << " " << setData << std::endl; + throw GaudiException(e.what(), "FATAL", StatusCode::FAILURE); + } + + const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + _bZ = field.z(); + + std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& setlayers = setData->layers; + _nLayers = setlayers.size(); + _SETgeo.resize(_nLayers); + + double strip_angle_deg = setData->angleStrip/CLHEP::degree; + //std::cout << "=============SET strip angle: " << strip_angle_deg << "==============" << std::endl; + for( int layer=0; layer < _nLayers; ++layer){ + dd4hep::rec::ZPlanarData::LayerLayout& pSETLayerLayout = setlayers[layer]; + + _SETgeo[layer].nLadders = pSETLayerLayout.ladderNumber; + _SETgeo[layer].phi0 = pSETLayerLayout.phi0; + _SETgeo[layer].dphi = 2*M_PI / _SETgeo[layer].nLadders; + _SETgeo[layer].senRMin = pSETLayerLayout.distanceSensitive*CLHEP::cm; + _SETgeo[layer].supRMin = pSETLayerLayout.distanceSupport*CLHEP::cm; + _SETgeo[layer].length = pSETLayerLayout.zHalfSensitive*2.0*CLHEP::cm; // note: gear for historical reasons uses the halflength + _SETgeo[layer].width = pSETLayerLayout.widthSensitive*CLHEP::cm; + _SETgeo[layer].offset = pSETLayerLayout.offsetSensitive*CLHEP::cm; + _SETgeo[layer].senThickness = pSETLayerLayout.thicknessSensitive*CLHEP::cm; + _SETgeo[layer].supThickness = pSETLayerLayout.thicknessSupport*CLHEP::cm; + _SETgeo[layer].nSensorsPerLadder = pSETLayerLayout.sensorsPerLadder; + _SETgeo[layer].sensorLength = _SETgeo[layer].length / _SETgeo[layer].nSensorsPerLadder; + + if (_isStripDetector) { + _SETgeo[layer].stripAngle = strip_angle_deg * M_PI/180 ; + } else { + _SETgeo[layer].stripAngle = 0.0 ; + } + //std::cout << _SETgeo[layer].nLadders << " " << _SETgeo[layer].phi0 << " "<< _SETgeo[layer].dphi << " " << _SETgeo[layer].senRMin << " " << _SETgeo[layer].supRMin << " " + // << _SETgeo[layer].length << " " << _SETgeo[layer].width << " " << _SETgeo[layer].offset << " " << _SETgeo[layer].senThickness << " " << _SETgeo[layer].supThickness << " " + // << _SETgeo[layer].nSensorsPerLadder << " " << _SETgeo[layer].sensorLength << " " << pSETLayerLayout.lengthSensor << std::endl; + } } diff --git a/Utilities/KalDet/src/ild/sit/ILDSITKalDetector.cc b/Utilities/KalDet/src/ild/sit/ILDSITKalDetector.cc index e58b9d6e773e9afabad21e3179be18e2891f4f69..1c409d64a20b13bca6a638dc003c256d5977b3b2 100644 --- a/Utilities/KalDet/src/ild/sit/ILDSITKalDetector.cc +++ b/Utilities/KalDet/src/ild/sit/ILDSITKalDetector.cc @@ -9,6 +9,11 @@ #include <UTIL/BitField64.h> #include <UTIL/ILDConf.h> +#include "DetInterface/IGeoSvc.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" +#include "CLHEP/Units/SystemOfUnits.h" + #include <gear/GEAR.h> #include "gear/BField.h" #include <gear/ZPlanarParameters.h> @@ -22,21 +27,25 @@ // #include "streamlog/streamlog.h" -ILDSITKalDetector::ILDSITKalDetector( const gear::GearMgr& gearMgr ) +ILDSITKalDetector::ILDSITKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc ) : TVKalDetector(300) // SJA:FIXME initial size, 300 looks reasonable for ILD, though this would be better stored as a const somewhere { // streamlog_out(DEBUG4) << "ILDSITKalDetector building SIT detector using GEAR " << std::endl ; - MaterialDataBase::Instance().registerForService(gearMgr); - + MaterialDataBase::Instance().registerForService(gearMgr, geoSvc); + TMaterial & air = *MaterialDataBase::Instance().getMaterial("air"); TMaterial & silicon = *MaterialDataBase::Instance().getMaterial("silicon"); TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("carbon"); - - this->setupGearGeom(gearMgr) ; - + if(geoSvc){ + this->setupGearGeom(geoSvc); + } + else{ + this->setupGearGeom(gearMgr) ; + } + if (_isStripDetector) { // streamlog_out(DEBUG4) << "\t\t building SIT detector as STRIP Detector." << std::endl ; } else { @@ -241,7 +250,7 @@ void ILDSITKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ // if this is not done then the exposed areas of the support would leave a carbon - air boundary, // which if traversed in the reverse direction to the next boundary then the track would be propagated through carbon // for a significant distance - + //std::cout << "=============SIT strip angle: " << strip_angle_deg << "==============" << std::endl; for( int layer=0; layer < _nLayers; ++layer){ @@ -271,7 +280,9 @@ void ILDSITKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ } else { _SITgeo[layer].stripAngle = 0.0 ; } - + //std::cout << _SITgeo[layer].nLadders << " " << _SITgeo[layer].phi0 << " "<< _SITgeo[layer].dphi << " " << _SITgeo[layer].senRMin << " " << _SITgeo[layer].supRMin << " " + // << _SITgeo[layer].length << " " << _SITgeo[layer].width << " " << _SITgeo[layer].offset << " " << _SITgeo[layer].senThickness << " " << _SITgeo[layer].supThickness << " " + // << _SITgeo[layer].nSensorsPerLadder << " " << _SITgeo[layer].sensorLength << std::endl; // streamlog_out(DEBUG0) << " layer = " << layer << std::endl; // streamlog_out(DEBUG0) << " nSensorsPerLadder = " << _SITgeo[layer].nSensorsPerLadder << std::endl; // streamlog_out(DEBUG0) << " sensorLength = " << _SITgeo[layer].sensorLength << std::endl; @@ -280,7 +291,61 @@ void ILDSITKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ } - - - +} + +void ILDSITKalDetector::setupGearGeom( IGeoSvc* geoSvc ){ + + dd4hep::DetElement world = geoSvc->getDD4HepGeo(); + dd4hep::DetElement sit; + const std::map<std::string, dd4hep::DetElement>& subs = world.children(); + for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){ + if(it->first!="SIT") continue; + sit = it->second; + } + dd4hep::rec::ZPlanarData* sitData = nullptr; + try{ + sitData = sit.extension<dd4hep::rec::ZPlanarData>(); + } + catch(std::runtime_error& e){ + std::cout << e.what() << " " << sitData << std::endl; + throw GaudiException(e.what(), "FATAL", StatusCode::FAILURE); + } + + const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + _bZ = field.z(); + + std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& sitlayers = sitData->layers; + _nLayers = sitlayers.size(); + _SITgeo.resize(_nLayers); + + double strip_angle_deg = sitData->angleStrip/CLHEP::degree; + if(strip_angle_deg!=0){ + _isStripDetector = true; + } + //std::cout << "=============SIT strip angle: " << strip_angle_deg << "==============" << std::endl; + for( int layer=0; layer < _nLayers; ++layer){ + dd4hep::rec::ZPlanarData::LayerLayout& pSITLayerLayout = sitlayers[layer]; + + _SITgeo[layer].nLadders = pSITLayerLayout.ladderNumber; + _SITgeo[layer].phi0 = pSITLayerLayout.phi0; + _SITgeo[layer].dphi = 2*M_PI / _SITgeo[layer].nLadders; + _SITgeo[layer].senRMin = pSITLayerLayout.distanceSensitive*CLHEP::cm; + _SITgeo[layer].supRMin = pSITLayerLayout.distanceSupport*CLHEP::cm; + _SITgeo[layer].length = pSITLayerLayout.zHalfSensitive*2.0*CLHEP::cm; // note: gear for historical reasons uses the halflength + _SITgeo[layer].width = pSITLayerLayout.widthSensitive*CLHEP::cm; + _SITgeo[layer].offset = pSITLayerLayout.offsetSensitive*CLHEP::cm; + _SITgeo[layer].senThickness = pSITLayerLayout.thicknessSensitive*CLHEP::cm; + _SITgeo[layer].supThickness = pSITLayerLayout.thicknessSupport*CLHEP::cm; + _SITgeo[layer].nSensorsPerLadder = pSITLayerLayout.sensorsPerLadder; + _SITgeo[layer].sensorLength = _SITgeo[layer].length / _SITgeo[layer].nSensorsPerLadder; + + if (_isStripDetector) { + _SITgeo[layer].stripAngle = strip_angle_deg * M_PI/180 ; + } else { + _SITgeo[layer].stripAngle = 0.0 ; + } + //std::cout << _SITgeo[layer].nLadders << " " << _SITgeo[layer].phi0 << " "<< _SITgeo[layer].dphi << " " << _SITgeo[layer].senRMin << " " << _SITgeo[layer].supRMin << " " + // << _SITgeo[layer].length << " " << _SITgeo[layer].width << " " << _SITgeo[layer].offset << " " << _SITgeo[layer].senThickness << " " << _SITgeo[layer].supThickness << " " + // << _SITgeo[layer].nSensorsPerLadder << " " << _SITgeo[layer].sensorLength << " " << pSITLayerLayout.lengthSensor << std::endl; + } } diff --git a/Utilities/KalDet/src/ild/sit/ILDSITKalDetector.h b/Utilities/KalDet/src/ild/sit/ILDSITKalDetector.h index 639e1b27bb76512ad05c8570814b218a9480cb5b..d9eb3059b680eb16823cef7de3b315659c056272 100644 --- a/Utilities/KalDet/src/ild/sit/ILDSITKalDetector.h +++ b/Utilities/KalDet/src/ild/sit/ILDSITKalDetector.h @@ -16,18 +16,20 @@ namespace gear{ class GearMgr ; } +class IGeoSvc; class ILDSITKalDetector : public TVKalDetector { public: /** Initialize the SIT from GEAR */ - ILDSITKalDetector( const gear::GearMgr& gearMgr ); + ILDSITKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc ); private: void setupGearGeom( const gear::GearMgr& gearMgr ) ; + void setupGearGeom( IGeoSvc* geoSvc ); int _nLayers ; double _bZ ; diff --git a/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc b/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc index 3b59a6af44c926aaa2d88a9b8217b8e98bdec813..98a6f247b86b2f4cee58936846535205c0f5f96c 100644 --- a/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc +++ b/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc @@ -16,6 +16,10 @@ #include <sstream> #include <cmath> +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" +#include "CLHEP/Units/SystemOfUnits.h" + #include "gear/GEAR.h" #include "gear/BField.h" #include "gearimpl/Util.h" @@ -26,28 +30,33 @@ ILDSupportKalDetector::ILDSupportKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc ) : TVKalDetector(10) { - + Double_t bz; + std::vector<double> z, rInner, rOuter; // streamlog_out(DEBUG1) << "ILDSupportKalDetector building beampipe using GEAR " << std::endl ; - - const gear::GearParameters& pBeamPipe = gearMgr.getGearParameters("BeamPipe"); - const Double_t bz = gearMgr.getBField().at( gear::Vector3D( 0.,0.,0.) ).z() ; - - // switch gear to dd4hep::rec, once validated, gear will been removed. - const dd4hep::rec::ConicalSupportData* pBeamPipeData = geoSvc->getBeamPipeData(); - const std::vector<dd4hep::rec::ConicalSupportData::Section>& sections = pBeamPipeData->sections; - std::cout << "======================BeamPipe===================" << std::endl; - for(int i=0;i<sections.size();i++){ - std::cout << sections[i].zPos << " " << sections[i].rInner << " " << sections[i].rOuter << std::endl; + if(geoSvc){ + const dd4hep::rec::ConicalSupportData* pBeamPipeData = geoSvc->getBeamPipeData(); + const std::vector<dd4hep::rec::ConicalSupportData::Section>& sections = pBeamPipeData->sections; + const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + bz = field.z(); + for(int i=0;i<sections.size();i++){ + z.push_back(sections[i].zPos); + rInner.push_back(sections[i].rInner); + rOuter.push_back(sections[i].rOuter); + } } - std::cout << "vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv" << std::endl; - // the beampipe is supposed to be a chain of cut cones (cut, that means the spike is cut off, also called a cone frustum). - // as they are connected, RStart[i] == REnd[i-1]. With this all we need are the z values and the radii at the place. - const std::vector<double> z = pBeamPipe.getDoubleVals("Z"); - const std::vector<double> rInner = pBeamPipe.getDoubleVals("RInner"); //inner radius of the cone - const std::vector<double> rOuter = pBeamPipe.getDoubleVals("ROuter"); //outer radius of the cone - for(int i=0;i<sections.size();i++){ - std::cout << z[i] << " " << rInner[i] << " " << rOuter[i] << std::endl; + else{ + const gear::GearParameters& pBeamPipe = gearMgr.getGearParameters("BeamPipe"); + bz = gearMgr.getBField().at( gear::Vector3D( 0.,0.,0.) ).z() ; + + // the beampipe is supposed to be a chain of cut cones (cut, that means the spike is cut off, also called a cone frustum). + // as they are connected, RStart[i] == REnd[i-1]. With this all we need are the z values and the radii at the place. + z = pBeamPipe.getDoubleVals("Z"); + rInner = pBeamPipe.getDoubleVals("RInner"); //inner radius of the cone + rOuter = pBeamPipe.getDoubleVals("ROuter"); //outer radius of the cone } + //for(int i=0;i<sections.size();i++){ + // std::cout << z[i] << " " << rInner[i] << " " << rOuter[i] << std::endl; + //} MaterialDataBase::Instance().registerForService(gearMgr, geoSvc); TMaterial & beam = *MaterialDataBase::Instance().getMaterial("beam"); TMaterial & air = *MaterialDataBase::Instance().getMaterial("air"); diff --git a/Utilities/KalDet/src/ild/tpc/ILDTPCKalDetector.cc b/Utilities/KalDet/src/ild/tpc/ILDTPCKalDetector.cc index cb41bf467d2a208fbc8b6ff0587f08bf5aca2ba7..fbb5f380b9cf0b9d703e0bb334d835dee1321134 100644 --- a/Utilities/KalDet/src/ild/tpc/ILDTPCKalDetector.cc +++ b/Utilities/KalDet/src/ild/tpc/ILDTPCKalDetector.cc @@ -10,6 +10,11 @@ #include <sstream> +#include "DetInterface/IGeoSvc.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" +#include "CLHEP/Units/SystemOfUnits.h" + #include "gear/GEAR.h" #include "gear/BField.h" #include "gear/TPCParameters.h" @@ -22,37 +27,69 @@ // #include "streamlog/streamlog.h" -ILDTPCKalDetector::ILDTPCKalDetector( const gear::GearMgr& gearMgr ) : +ILDTPCKalDetector::ILDTPCKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc ) : TVKalDetector(250) // SJA:FIXME initial size, 250 looks reasonable for ILD, though this would be better stored as a const somewhere { - + Double_t bz; + Int_t nlayers; + Double_t lhalf, rstep, rmin, rtub, outerr, inthick, outthick; // streamlog_out(DEBUG1) << "ILDTPCKalDetector building TPC detector using GEAR " << std::endl ; - - const gear::TPCParameters& tpcParams = gearMgr.getTPCParameters(); - - const gear::PadRowLayout2D& pL = tpcParams.getPadLayout() ; - - // streamlog_out(DEBUG1) << "ILDTPCKalDetector - got padlayout with nLayers = " << pL.getNRows() << std::endl ; + if(geoSvc){ + dd4hep::DetElement world = geoSvc->getDD4HepGeo(); + dd4hep::DetElement tpc; + const std::map<std::string, dd4hep::DetElement>& subs = world.children(); + for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){ + if(it->first!="TPC") continue; + tpc = it->second; + } + dd4hep::rec::FixedPadSizeTPCData* tpcData = nullptr; + try{ + tpcData = tpc.extension<dd4hep::rec::FixedPadSizeTPCData>(); + } + catch(std::runtime_error& e){ + std::cout << e.what() << " " << tpcData << std::endl; + throw GaudiException(e.what(), "FATAL", StatusCode::FAILURE); + } - const Double_t bz = gearMgr.getBField().at( gear::Vector3D( 0.,0.,0.) ).z() ; - - static const Int_t nlayers = pL.getNRows() ; // n rows - static const Double_t lhalf = tpcParams.getMaxDriftLength() ; // half length - - static const Double_t rstep = pL.getRowHeight(0) ; // step length of radius - - // assuming that this is the radius of the first measurment layer .... - static const Double_t rmin = tpcParams.getPlaneExtent()[0] + rstep/2. ; // minimum radius - - // streamlog_out( DEBUG0 ) << tpcParams << std::endl ; + const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + bz = field.z(); + nlayers = tpcData->maxRow; + lhalf = tpcData->driftLength*CLHEP::cm; + rstep = tpcData->padHeight*CLHEP::cm; + rmin = tpcData->rMinReadout*CLHEP::cm + 0.5*rstep; + rtub = tpcData->rMin*CLHEP::cm; + outerr = tpcData->rMax*CLHEP::cm; + inthick = tpcData->innerWallThickness*CLHEP::cm; + outthick = tpcData->outerWallThickness*CLHEP::cm; + //std::cout << "TPC: " << nlayers << " " << lhalf << " " << rstep << " " << rmin << " " << rtub << " " << outerr << " " << inthick << " " << outthick << std::endl; + } + else{ + const gear::TPCParameters& tpcParams = gearMgr.getTPCParameters(); - static const Double_t rtub = tpcParams.getDoubleVal("tpcInnerRadius") ; // inner r of support tube - static const Double_t outerr = tpcParams.getDoubleVal("tpcOuterRadius") ; // outer radius of TPC + const gear::PadRowLayout2D& pL = tpcParams.getPadLayout() ; + + // std::cout << "ILDTPCKalDetector - got padlayout with nLayers = " << pL.getNRows() << std::endl ; + + bz = gearMgr.getBField().at( gear::Vector3D( 0.,0.,0.) ).z() ; - static const Double_t inthick = tpcParams.getDoubleVal("tpcInnerWallThickness") ; // thickness of inner shell - static const Double_t outthick = tpcParams.getDoubleVal("tpcOuterWallThickness") ; // thickness of outer shell + nlayers = pL.getNRows() ; // n rows + lhalf = tpcParams.getMaxDriftLength() ; // half length + rstep = pL.getRowHeight(0) ; // step length of radius + + // assuming that this is the radius of the first measurment layer .... + rmin = tpcParams.getPlaneExtent()[0] + rstep/2. ; // minimum radius + + // std::cout << tpcParams << std::endl ; + + rtub = tpcParams.getDoubleVal("tpcInnerRadius") ; // inner r of support tube + outerr = tpcParams.getDoubleVal("tpcOuterRadius") ; // outer radius of TPC + + inthick = tpcParams.getDoubleVal("tpcInnerWallThickness") ; // thickness of inner shell + outthick = tpcParams.getDoubleVal("tpcOuterWallThickness") ; // thickness of outer shell + //std::cout << "TPC: " << nlayers << " " << lhalf << " " << rstep << " " << rmin << " " << rtub << " " << outerr << " " << inthick << " " << outthick << std::endl; + } - MaterialDataBase::Instance().registerForService(gearMgr); + MaterialDataBase::Instance().registerForService(gearMgr, geoSvc); TMaterial & air = *MaterialDataBase::Instance().getMaterial("air"); TMaterial & tpcgas = *MaterialDataBase::Instance().getMaterial("tpcgas"); // TMaterial & aluminium = *MaterialDataBase::Instance().getMaterial("aluminium"); diff --git a/Utilities/KalDet/src/ild/tpc/ILDTPCKalDetector.h b/Utilities/KalDet/src/ild/tpc/ILDTPCKalDetector.h index e5ab8db97c723357f2fbd8e964ab319ec7102ca6..a249f22ce1a7133dc88135ec233c785b54d777b7 100644 --- a/Utilities/KalDet/src/ild/tpc/ILDTPCKalDetector.h +++ b/Utilities/KalDet/src/ild/tpc/ILDTPCKalDetector.h @@ -14,12 +14,13 @@ namespace gear{ class GearMgr ; } +class IGeoSvc; class ILDTPCKalDetector : public TVKalDetector { public: /** Initialize the TPC from GEAR */ - ILDTPCKalDetector( const gear::GearMgr& gearMgr ); + ILDTPCKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc=0 ); private: diff --git a/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc b/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc index a820db54985ae8936acecc432c08577a3f9324de..edf608f8d762c5763290671cf44a573a9c288e07 100644 --- a/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc +++ b/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc @@ -10,6 +10,11 @@ #include <UTIL/BitField64.h> #include <UTIL/ILDConf.h> +#include "DetInterface/IGeoSvc.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" +#include "CLHEP/Units/SystemOfUnits.h" + #include <gear/GEAR.h> #include "gear/BField.h" #include <gearimpl/ZPlanarParametersImpl.h> @@ -44,8 +49,12 @@ ILDVXDKalDetector::ILDVXDKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geo _vxd_Cryostat.exists = false; - this->setupGearGeom(gearMgr, geoSvc) ; - + if(geoSvc){ + this->setupGearGeom(geoSvc) ; + } + else{ + this->setupGearGeom(gearMgr) ; + } //--The Ladder structure (realistic ladder)-- int nLadders; @@ -315,18 +324,84 @@ ILDVXDKalDetector::ILDVXDKalDetector( const gear::GearMgr& gearMgr, IGeoSvc* geo SetOwner(); } +void ILDVXDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ + const gear::VXDParameters& pVXDDetMain = gearMgr.getVXDParameters(); + const gear::VXDLayerLayout& pVXDLayerLayout = pVXDDetMain.getVXDLayerLayout(); + _bZ = gearMgr.getBField().at( gear::Vector3D( 0.,0.,0.) ).z() ; + _nLayers = pVXDLayerLayout.getNLayers(); + _VXDgeo.resize(_nLayers); -void ILDVXDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc){ + for( int layer=0; layer < _nLayers; ++layer){ + _VXDgeo[layer].nLadders = pVXDLayerLayout.getNLadders(layer); + _VXDgeo[layer].phi0 = pVXDLayerLayout.getPhi0(layer); + _VXDgeo[layer].dphi = 2*M_PI / _VXDgeo[layer].nLadders; + _VXDgeo[layer].senRMin = pVXDLayerLayout.getSensitiveDistance(layer); + _VXDgeo[layer].supRMin = pVXDLayerLayout.getLadderDistance(layer); + _VXDgeo[layer].length = pVXDLayerLayout.getSensitiveLength(layer) * 2.0 ; // note: gear for historical reasons uses the halflength + _VXDgeo[layer].width = pVXDLayerLayout.getSensitiveWidth(layer); + _VXDgeo[layer].offset = pVXDLayerLayout.getSensitiveOffset(layer); + _VXDgeo[layer].senThickness = pVXDLayerLayout.getSensitiveThickness(layer); + _VXDgeo[layer].supThickness = pVXDLayerLayout.getLadderThickness(layer); + //std::cout << layer << ": " << _VXDgeo[layer].nLadders << " " << _VXDgeo[layer].phi0 << " " << _VXDgeo[layer].dphi << " " << _VXDgeo[layer].senRMin + // << " " << _VXDgeo[layer].supRMin << " " << _VXDgeo[layer].length << " " << _VXDgeo[layer].width << " " << _VXDgeo[layer].offset + // << " " << _VXDgeo[layer].senThickness << " " << _VXDgeo[layer].supThickness << std::endl; + } + + _relative_position_of_measurement_surface = 0.5 ; + + try { + _relative_position_of_measurement_surface = pVXDDetMain.getDoubleVal( "relative_position_of_measurement_surface" ); + } + catch (gear::UnknownParameterException& e) {} + try { + const gear::GearParameters& pVXDInfra = gearMgr.getGearParameters("VXDInfra"); + _vxd_Cryostat.alRadius = pVXDInfra.getDoubleVal( "CryostatAlRadius" ); + _vxd_Cryostat.alThickness = pVXDInfra.getDoubleVal( "CryostatAlThickness" ); + _vxd_Cryostat.alInnerR = pVXDInfra.getDoubleVal( "CryostatAlInnerR" ); + _vxd_Cryostat.alZEndCap = pVXDInfra.getDoubleVal( "CryostatAlZEndCap" ); + _vxd_Cryostat.alHalfZ = pVXDInfra.getDoubleVal( "CryostatAlHalfZ" ); + + _vxd_Cryostat.shellInnerR = pVXDDetMain.getShellInnerRadius(); + _vxd_Cryostat.shellThickness = pVXDDetMain.getShellOuterRadius() - _vxd_Cryostat.shellInnerR; + _vxd_Cryostat.shelllHalfZ = pVXDDetMain.getShellHalfLength(); + + _vxd_Cryostat.exists = true; + //std::cout << "VXDInfra: " << _vxd_Cryostat.alRadius << " " << _vxd_Cryostat.alThickness << " " << _vxd_Cryostat.alInnerR << " " << _vxd_Cryostat.alZEndCap << " " + // << _vxd_Cryostat.alHalfZ << " " << _vxd_Cryostat.shellInnerR << " " << _vxd_Cryostat.shellThickness << " " << _vxd_Cryostat.shelllHalfZ << std::endl; + } + catch (gear::UnknownParameterException& e) { + std::cout << e.what() << std::endl ; + _vxd_Cryostat.exists = false; + + } +} + + +void ILDVXDKalDetector::setupGearGeom( IGeoSvc* geoSvc){ + /* + dd4hep::DetElement world = geoSvc->getDD4HepGeo(); + dd4hep::DetElement vxd; + const std::map<std::string, dd4hep::DetElement>& subs = world.children(); + for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){ + if(it->first!="VXD") continue; + vxd = it->second; + } + dd4hep::rec::ZPlanarData* vxdData = nullptr; + try{ + vxdData = vxd.extension<dd4hep::rec::ZPlanarData>(); + } + catch(std::runtime_error& e){ + std::cout << e.what() << " " << vxdData << std::endl; + throw GaudiException(e.what(), "FATAL", StatusCode::FAILURE); + } + */ + const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + _bZ = field.z(); - const gear::VXDParameters& pVXDDetMainObject = gearMgr.getVXDParameters(); - const gear::VXDParameters* pVXDDetMain2 = &pVXDDetMainObject; const gear::ZPlanarParametersImpl* pVXDDetMain = geoSvc->getVXDParameters(); const gear::VXDLayerLayout& pVXDLayerLayout = pVXDDetMain->getVXDLayerLayout(); - const gear::VXDLayerLayout& pVXDLayerLayout2 = pVXDDetMain2->getVXDLayerLayout(); - - _bZ = gearMgr.getBField().at( gear::Vector3D( 0.,0.,0.) ).z() ; _nLayers = pVXDLayerLayout.getNLayers(); _VXDgeo.resize(_nLayers); @@ -347,42 +422,16 @@ void ILDVXDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr, IGeoSvc* ge _VXDgeo[layer].offset = pVXDLayerLayout.getSensitiveOffset(layer); _VXDgeo[layer].senThickness = pVXDLayerLayout.getSensitiveThickness(layer); _VXDgeo[layer].supThickness = pVXDLayerLayout.getLadderThickness(layer); - std::cout << layer << ": " << _VXDgeo[layer].nLadders << " " << _VXDgeo[layer].phi0 << " " << _VXDgeo[layer].dphi << " " << _VXDgeo[layer].senRMin - << " " << _VXDgeo[layer].supRMin << " " << _VXDgeo[layer].length << " " << _VXDgeo[layer].width << " " << _VXDgeo[layer].offset - << " " << _VXDgeo[layer].senThickness << " " << _VXDgeo[layer].supThickness << std::endl; - std::cout << layer << ": " << pVXDLayerLayout2.getNLadders(layer) << " " << pVXDLayerLayout2.getPhi0(layer) << " " << 2*M_PI/pVXDLayerLayout2.getNLadders(layer) - << " " << pVXDLayerLayout2.getSensitiveDistance(layer) << " " << pVXDLayerLayout2.getLadderDistance(layer) - << " " << pVXDLayerLayout2.getSensitiveLength(layer)*2 << " " << pVXDLayerLayout2.getSensitiveWidth(layer) - << " " << pVXDLayerLayout2.getSensitiveOffset(layer) << " " << pVXDLayerLayout2.getSensitiveThickness(layer) - << " " << pVXDLayerLayout2.getLadderThickness(layer) << std::endl; + //std::cout << layer << ": " << _VXDgeo[layer].nLadders << " " << _VXDgeo[layer].phi0 << " " << _VXDgeo[layer].dphi << " " << _VXDgeo[layer].senRMin + // << " " << _VXDgeo[layer].supRMin << " " << _VXDgeo[layer].length << " " << _VXDgeo[layer].width << " " << _VXDgeo[layer].offset + // << " " << _VXDgeo[layer].senThickness << " " << _VXDgeo[layer].supThickness << std::endl; } // by default, we put the measurement surface in the middle of the sensitive // layer, this can optionally be changed, e.g. in the case of the FPCCD where the // epitaxial layer is 15 mu thick (in a 50 mu wafer) _relative_position_of_measurement_surface = 0.5 ; - - try { - - _relative_position_of_measurement_surface = pVXDDetMain->getDoubleVal( "relative_position_of_measurement_surface" ); - - // streamlog_out(DEBUG) << " ILDVXDKalDetector::setupGearGeom: relative_position_of_measurement_surface parameter is provided : " << _relative_position_of_measurement_surface << std::endl ; - - } catch (gear::UnknownParameterException& e) {} - // Cryostat try { - - //const gear::GearParameters& pVXDInfraObject = gearMgr.getGearParameters("VXDInfra"); - //const gear::GearParameters* pVXDInfra = &pVXDInfraObject; - //change Gear to GeoSvc - //const gear::GearParametersImpl* pVXDInfra = geoSvc->getDetParameters("VXDInfra"); - //_vxd_Cryostat.alRadius = pVXDInfra->getDoubleVal( "CryostatAlRadius" ); - //_vxd_Cryostat.alThickness = pVXDInfra->getDoubleVal( "CryostatAlThickness" ); - //_vxd_Cryostat.alInnerR = pVXDInfra->getDoubleVal( "CryostatAlInnerR" ); - //_vxd_Cryostat.alZEndCap = pVXDInfra->getDoubleVal( "CryostatAlZEndCap" ); - //_vxd_Cryostat.alHalfZ = pVXDInfra->getDoubleVal( "CryostatAlHalfZ" ); - //change GearParametersImpl to get directly - //const std::map<std::string,double>& vxdInfra = geoSvc->getDetParameters("VXDInfra"); _vxd_Cryostat.alRadius = geoSvc->getDetParameter("VXDInfra","CryostatAlRadius"); _vxd_Cryostat.alThickness = geoSvc->getDetParameter("VXDInfra","CryostatAlThickness"); _vxd_Cryostat.alInnerR = geoSvc->getDetParameter("VXDInfra","CryostatAlInnerR"); @@ -394,10 +443,10 @@ void ILDVXDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr, IGeoSvc* ge _vxd_Cryostat.shelllHalfZ = pVXDDetMain->getShellHalfLength(); _vxd_Cryostat.exists = true; + //std::cout << "VXDInfra: " << _vxd_Cryostat.alRadius << " " << _vxd_Cryostat.alThickness << " " << _vxd_Cryostat.alInnerR << " " << _vxd_Cryostat.alZEndCap << " " + // << _vxd_Cryostat.alHalfZ << " " << _vxd_Cryostat.shellInnerR << " " << _vxd_Cryostat.shellThickness << " " << _vxd_Cryostat.shelllHalfZ << std::endl; } - //catch (gear::UnknownParameterException& e) { catch (std::runtime_error& e) { - //streamlog_out( WARNING ) << "ILDVXDKalDetector Cryostat values not found in GEAR file " << std::endl ; std::cout << e.what() << std::endl ; _vxd_Cryostat.exists = false; diff --git a/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.h b/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.h index c0205e149f55b4ca2340cdfc7549b17fd020cdbc..05840c8c85477a4083404a22dfcadb652d7a0e60 100644 --- a/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.h +++ b/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.h @@ -27,7 +27,8 @@ public: private: - void setupGearGeom( const gear::GearMgr& gearMgr, IGeoSvc* geoSvc) ; + void setupGearGeom( const gear::GearMgr& gearMgr ); + void setupGearGeom( IGeoSvc* geoSvc) ; int _nLayers ; double _bZ ; diff --git a/Utilities/KiTrack/CMakeLists.txt b/Utilities/KiTrack/CMakeLists.txt index 6ebfefe261de9528ab433639a37252e6dc1b95ed..22d2f25da413adfae3d0013bd3d5136046f79ce5 100644 --- a/Utilities/KiTrack/CMakeLists.txt +++ b/Utilities/KiTrack/CMakeLists.txt @@ -5,6 +5,7 @@ find_package(ROOT REQUIRED) find_package(GSL REQUIRED) find_package(EDM4HEP REQUIRED) find_package(LCIO REQUIRED) +find_package(DD4hep COMPONENTS DDCore DDRec REQUIRED) gaudi_depends_on_subdirs(Service/TrackSystemSvc Utilities/DataHelper)