diff --git a/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.cpp b/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.cpp
index 773f1998a17af850fd5f06d42b1ef3d5d0d5a79a..b2972273323733713eae3e3eecd54f6ea72521dd 100644
--- a/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.cpp
+++ b/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.cpp
@@ -17,13 +17,13 @@ DECLARE_COMPONENT(KalTestTool)
StatusCode KalTestTool::initialize() {
StatusCode sc;
always() << m_fitterName << endmsg;
- if (m_fitterName=="KalTest") {
+ if (m_fitterName=="KalTest" || m_fitterName=="DDKalTest") {
auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc");
if (!_trackSystemSvc) {
error() << "Failed to find TrackSystemSvc ..." << endmsg;
return StatusCode::FAILURE;
}
- m_factoryMarlinTrk = _trackSystemSvc->getTrackSystem(this);
+ m_factoryMarlinTrk = _trackSystemSvc->getTrackSystem(this, m_fitterName.value());
m_factoryMarlinTrk->setOption(MarlinTrk::IMarlinTrkSystem::CFG::useQMS, m_useQMS);
m_factoryMarlinTrk->setOption(MarlinTrk::IMarlinTrkSystem::CFG::usedEdx, m_usedEdx);
m_factoryMarlinTrk->setOption(MarlinTrk::IMarlinTrkSystem::CFG::useSmoothing, m_useSmoothing);
@@ -65,7 +65,7 @@ int KalTestTool::Fit(edm4hep::MutableTrack track, std::vector<edm4hep::TrackerHi
return 0;
}
- if (m_fitterName=="KalTest") {
+ if (m_fitterName=="KalTest" || m_fitterName=="DDKalTest") {
debug() << "start..." << endmsg;
std::shared_ptr<MarlinTrk::IMarlinTrack> marlinTrack(m_factoryMarlinTrk->createTrack());
debug() << "created MarlinKalTestTrack" << endmsg;
@@ -87,7 +87,7 @@ int KalTestTool::Fit(edm4hep::MutableTrack track, std::vector<edm4hep::TrackerHi
return 0;
}
- if (m_fitterName=="KalTest") {
+ if (m_fitterName=="KalTest" || m_fitterName=="DDKalTest") {
debug() << "start..." << endmsg;
std::shared_ptr<MarlinTrk::IMarlinTrack> marlinTrack(m_factoryMarlinTrk->createTrack());
debug() << "created MarlinKalTestTrack" << endmsg;
diff --git a/Service/TrackSystemSvc/CMakeLists.txt b/Service/TrackSystemSvc/CMakeLists.txt
index 67710976b6990942ab95d781ccf8fbf9ad1ca90a..2baf051a903ff9b62d89d8ce446079bcce946b5f 100644
--- a/Service/TrackSystemSvc/CMakeLists.txt
+++ b/Service/TrackSystemSvc/CMakeLists.txt
@@ -8,11 +8,15 @@ gaudi_add_library(TrackSystemSvcLib
src/LCIOTrackPropagators.cc
src/MarlinKalTest.cc
src/MarlinKalTestTrack.cc
+ src/MarlinDDKalTest.cc
+ src/MarlinDDKalTestTrack.cc
src/MarlinTrkUtils.cc
- LINK DataHelperLib
- KalTestLib
- KalDetLib
+ LINK DataHelperLib
+ ${KalTest_LIBRARIES}
+ KalDetLib
+ ${DDKalTest_LIBRARIES}
+ ${streamlog_LIBRARIES}
GearSvc
Gaudi::GaudiKernel
${ROOT_LIBRARIES}
@@ -22,6 +26,14 @@ gaudi_add_library(TrackSystemSvcLib
EDM4HEP::edm4hep EDM4HEP::edm4hepDict
)
+target_include_directories(TrackSystemSvcLib
+ PUBLIC ${KalTest_INCLUDE_DIRS}
+ ${DDKalTest_INCLUDE_DIRS}
+ ${streamlog_INCLUDE_DIRS}
+)
+
+add_dependencies(TrackSystemSvcLib DDKalTest)
+
gaudi_add_module(TrackSystemSvc
SOURCES src/TrackSystemSvc.cpp
LINK TrackSystemSvcLib
diff --git a/Service/TrackSystemSvc/include/TrackSystemSvc/IMarlinTrack.h b/Service/TrackSystemSvc/include/TrackSystemSvc/IMarlinTrack.h
index 99121fbfab145048b5b1b0fddba2bb19e287b67f..193ef29debcff912961194337a35ad5b71c4fedd 100644
--- a/Service/TrackSystemSvc/include/TrackSystemSvc/IMarlinTrack.h
+++ b/Service/TrackSystemSvc/include/TrackSystemSvc/IMarlinTrack.h
@@ -218,8 +218,11 @@ namespace MarlinTrk{
* returning intersection point in global coordinates via reference
*/
virtual int intersectionWithDetElement( int detEementID, edm4hep::TrackerHit& hit, edm4hep::Vector3d& point, int mode=modeClosest ) = 0 ;
-
-
+
+ /** Dump this track to a string for debugging - implementation dependant.
+ */
+ virtual std::string toString() ;
+
protected:
private:
diff --git a/Service/TrackSystemSvc/include/TrackSystemSvc/ITrackSystemSvc.h b/Service/TrackSystemSvc/include/TrackSystemSvc/ITrackSystemSvc.h
index 88a2d31dd7a85403efdfd81c790e7f93bc876b27..87f135a1479758a80c1c7c5b6c3c2a8d110a3671 100644
--- a/Service/TrackSystemSvc/include/TrackSystemSvc/ITrackSystemSvc.h
+++ b/Service/TrackSystemSvc/include/TrackSystemSvc/ITrackSystemSvc.h
@@ -13,7 +13,7 @@ class ITrackSystemSvc: virtual public IService {
virtual ~ITrackSystemSvc() = default;
//Get the track manager
- virtual MarlinTrk::IMarlinTrkSystem* getTrackSystem(void* address=0) = 0;
+ virtual MarlinTrk::IMarlinTrkSystem* getTrackSystem(void* address=0, std::string type="KalTest") = 0;
virtual void removeTrackSystem(void* address=0) = 0;
};
diff --git a/Service/TrackSystemSvc/src/HelixTrack.cc b/Service/TrackSystemSvc/src/HelixTrack.cc
index 8901f5034ad92ba556baca32b90490cf8eac81a6..63aa06e4b18b458eb1d47aa5d3cbedb010a5a94d 100644
--- a/Service/TrackSystemSvc/src/HelixTrack.cc
+++ b/Service/TrackSystemSvc/src/HelixTrack.cc
@@ -3,7 +3,7 @@
#include "TrackSystemSvc/HelixTrack.h"
#include <cmath>
#include <TVector3.h>
-#include <kaltest/THelicalTrack.h>
+#include "kaltest/THelicalTrack.h"
#include <edm4hep/Vector3d.h> //plcio/DoubleThree.h>
//#include "streamlog/streamlog.h"
diff --git a/Service/TrackSystemSvc/src/IMarlinTrack.cc b/Service/TrackSystemSvc/src/IMarlinTrack.cc
index 494c171698b459352d22b5609393bffae3558ebd..d9c1ef6ab1afb9202e7842bb4da50369e6d562f0 100644
--- a/Service/TrackSystemSvc/src/IMarlinTrack.cc
+++ b/Service/TrackSystemSvc/src/IMarlinTrack.cc
@@ -34,5 +34,41 @@ namespace MarlinTrk{
default: return "UNKNOWN" ;
}
}
+
+ std::string IMarlinTrack::toString() {
+
+ std::stringstream str ;
+
+ int ndf ;
+ double chi2 ;
+ std::vector<std::pair<edm4hep::TrackerHit, double> > hits ;
+ edm4hep::TrackState ts ;
+ edm4hep::TrackerHit lastHit ;
+
+ getTrackerHitAtPositiveNDF( lastHit ) ;
+
+ getHitsInFit( hits ) ;
+
+ getTrackState( ts, chi2, ndf ) ;
+
+ str << " ------------------- MarlinDDKalTestTrack: ------------------------ " << std::endl ;
+
+ str << " ndf: " << ndf << std::endl
+ << " chi2: " << chi2 << std::endl
+ << " number of hits in fit : " << hits.size() << std::endl
+ << " last constraned hit id : " << ( lastHit.isAvailable() ? lastHit.id().collectionID*10000000 + lastHit.id().index : -9999 ) << std::endl ;
+
+ for( unsigned i=0,n= hits.size() ; i<n ; ++i ){
+
+ str << " hit at index: " << i << " " << hits[i].first << std::endl
+ << " track state : " << ts << std::endl ;
+ }
+
+ str << " --------------------- " << std::endl ;
+ str << " current track state :" << &ts << std::endl ;
+ str << " --------------------- " << std::endl ;
+
+ return str.str() ;
+ }
}
diff --git a/Service/TrackSystemSvc/src/MarlinDDKalTest.cc b/Service/TrackSystemSvc/src/MarlinDDKalTest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8772a384f7f3e0e8149315ea439f3b0e448191bd
--- /dev/null
+++ b/Service/TrackSystemSvc/src/MarlinDDKalTest.cc
@@ -0,0 +1,517 @@
+#include "MarlinDDKalTest.h"
+#include "MarlinDDKalTestTrack.h"
+
+#include "kaltest/TKalDetCradle.h"
+#include "kaltest/TVKalDetector.h"
+#include "kaltest/THelicalTrack.h"
+
+#include "DDKalTest/DDVMeasLayer.h"
+#include "DDKalTest/DDKalDetector.h"
+#include "DDKalTest/DDCylinderMeasLayer.h"
+
+// #include "DDKalTest/DDSupportKalDetector.h"
+// #include "DDKalTest/DDVXDKalDetector.h"
+
+//SJA:FIXME: only needed for storing the modules in the layers map
+#include <UTIL/BitField64.h>
+#include "UTIL/LCTrackerConf.h"
+
+#include "DDRec/SurfaceManager.h"
+
+#include <algorithm>
+#include <string>
+#include <math.h>
+#include <cmath>
+#include <fstream>
+
+#include <utility>
+
+#include "streamlog/streamlog.h"
+
+//#include "DDKalTest/DDMeasurementSurfaceStoreFiller.h"
+
+namespace MarlinTrk{
+
+
+ MarlinDDKalTest::MarlinDDKalTest() :
+ _ipLayer(NULL) {
+
+ streamlog_out( DEBUG4 ) << " MarlinDDKalTest - initializing the detector ..." << std::endl ;
+
+ _det = new TKalDetCradle ; // from kaltest. TKalDetCradle inherits from TObjArray ...
+ _det->SetOwner( true ) ; // takes care of deleting subdetector in the end ...
+
+ is_initialised = false;
+
+ this->registerOptions() ;
+
+ streamlog_out( DEBUG4 ) << " MarlinDDKalTest - established " << std::endl ;
+ }
+
+ MarlinDDKalTest::~MarlinDDKalTest(){
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+ _diagnostics.end();
+#endif
+ for( auto* ddKalDet : _detectors ) { delete ddKalDet; }
+ delete _det ;
+ }
+
+
+ void MarlinDDKalTest::setOption(unsigned CFGOption, bool val) {
+
+ IMarlinTrkSystem::setOption( CFGOption, val ) ;
+
+ switch( CFGOption ) {
+
+ case IMarlinTrkSystem::CFG::useQMS :
+ this->includeMultipleScattering( val ) ;
+ break ;
+ case IMarlinTrkSystem::CFG::usedEdx :
+ this->includeEnergyLoss( val ) ;
+ break ;
+ // IMarlinTrkSystem::CFG::useSmoothing handled directly in MarlinDDKalTestTrack
+ }
+
+ }
+
+
+ void MarlinDDKalTest::init() {
+
+
+ this->includeMultipleScattering( getOption(IMarlinTrkSystem::CFG::useQMS) ) ;
+
+ this->includeEnergyLoss( getOption(IMarlinTrkSystem::CFG::usedEdx) ) ;
+
+ streamlog_out( DEBUG5 ) << " -------------------------------------------------------------------------------- " << std::endl ;
+ streamlog_out( DEBUG5 ) << " MarlinDDKalTest::init() called with the following options : " << std::endl ;
+ streamlog_out( DEBUG5 ) << this->getOptions() ;
+ streamlog_out( DEBUG5 ) << " -------------------------------------------------------------------------------- " << std::endl ;
+
+ if( is_initialised ) {
+
+ streamlog_out( DEBUG5 ) << " MarlinDDKalTest::init() - already initialized - only options are set .. " << std::endl ;
+
+ return ;
+ }
+
+
+ streamlog_out( DEBUG5 ) << " ##################### MarlinDDKalTest::init() - initializing " << std::endl ;
+
+ dd4hep::Detector& lcdd = dd4hep::Detector::getInstance();
+
+ double minS = 1.e99;
+ DDCylinderMeasLayer* ipLayer = 0 ;
+
+
+ // for the tracking we get all tracking detectors and all passive detectors (beam pipe,...)
+
+ std::vector< dd4hep::DetElement> detectors = lcdd.detectors( "tracker" ) ;
+ const std::vector< dd4hep::DetElement>& passiveDets = lcdd.detectors( "passive" ) ;
+ const std::vector< dd4hep::DetElement>& calos = lcdd.detectors( "calorimeter" ) ;
+
+ detectors.reserve( detectors.size() + passiveDets.size() + calos.size() ) ;
+
+ std::copy( passiveDets.begin() , passiveDets.end() , std::back_inserter( detectors ) ) ;
+
+ for ( std::vector< dd4hep::DetElement>::const_iterator it=calos.begin() ; it != calos.end() ; ++it ){
+
+ // for ( std::vector< dd4hep::DetElement>::const_iterator it=passiveDets.begin() ; it != passiveDets.end() ; ++it ){
+
+ std::string name = it->name() ;
+ std::transform( name.begin() , name.end() , name.begin() , ::tolower ) ;
+ if( name.find( "ecal" ) != std::string::npos ){
+
+ detectors.push_back( *it ) ;
+ }
+ }
+
+
+
+ for ( std::vector< dd4hep::DetElement>::iterator it=detectors.begin() ; it != detectors.end() ; ++it ){
+
+ dd4hep::DetElement det = *it ;
+
+ streamlog_out( DEBUG5 ) << " MarlinDDKalTest::init() - creating DDKalDetector for : " << det.name() << std::endl ;
+
+ _detectors.push_back( new DDKalDetector( det ) );
+ DDKalDetector* kalDet = _detectors.back();
+
+ this->storeActiveMeasurementModuleIDs( kalDet ) ;
+
+ _det->Install( *kalDet ) ;
+
+
+ Int_t nLayers = kalDet->GetEntriesFast() ;
+
+
+ // --- keep the cylinder meas layer with smallest sorting policy (radius) as ipLayer
+ // fixme: this should be implemented in a more explicit way ...
+ for( int i=0; i < nLayers; ++i ) {
+ const TVSurface* tvs = static_cast<const TVSurface*>( kalDet->At( i ) );
+
+ double s = tvs->GetSortingPolicy() ;
+ if( s < minS && dynamic_cast< DDCylinderMeasLayer* > ( kalDet->At( i) ) ) {
+ minS = s ;
+ ipLayer = dynamic_cast< DDCylinderMeasLayer* > ( kalDet->At( i) ) ;
+ }
+ }
+
+ if( streamlog_level( DEBUG5 ) ) { // dump surfaces to text file
+
+ std::map< double, DDVMeasLayer*> smap ;
+ std::ofstream file ;
+ std::stringstream s ; s << "DDKalTest_" << det.name() << "_surfaces.txt" ;
+ file.open( s.str().c_str() , std::ofstream::out ) ;
+ lcio::BitField64 bf( UTIL::LCTrackerCellID::encoding_string() ) ;
+
+ for( unsigned i=0,N=kalDet->GetEntriesFast() ; i<N ;++i){
+ DDVMeasLayer* ml = dynamic_cast<DDVMeasLayer*> ( kalDet->At( i ) ) ;
+ TVSurface* surf = dynamic_cast<TVSurface*> ( kalDet->At( i ) ) ;
+ smap[ surf->GetSortingPolicy() ] = ml ;
+ }
+ for( std::map<double,DDVMeasLayer*>::iterator itm=smap.begin() ; itm!=smap.end() ; ++itm){
+ bf.setValue( itm->second->getCellIDs()[0] ) ;
+ file << " " << std::scientific << std::setw(10) << itm->first << "\t" << bf.valueString() << *itm->second->surface() << "\n" ;
+ }
+ file.close() ;
+ }
+
+ }
+
+
+
+ if( ipLayer) {
+
+ _ipLayer = ipLayer ;
+
+ streamlog_out( MESSAGE ) << " MarlinDDKalTest: install IP layer at radius : " << minS << std::endl ;
+ }
+ //-------------------------------------------------------------------------------
+
+
+ _det->Close() ; // close the cradle
+ //done in Close() _det->Sort() ; // sort meas. layers from inside to outside
+
+ streamlog_out( DEBUG4 ) << " MarlinDDKalTest - number of layers = " << _det->GetEntriesFast() << std::endl ;
+
+
+ if( streamlog_level( DEBUG ) ) {
+
+ lcio::BitField64 bf( UTIL::LCTrackerCellID::encoding_string() ) ;
+
+ for( unsigned i=0,N=_det->GetEntriesFast() ; i<N ;++i){
+
+ DDVMeasLayer* ml = dynamic_cast<DDVMeasLayer*> ( _det->At( i ) ) ;
+
+ bf.setValue( ml->getLayerID() ) ;
+
+ TVSurface* s = dynamic_cast<TVSurface*> ( _det->At( i ) ) ;
+
+ streamlog_out( DEBUG ) << " *** meas. layer : " << bf.valueString() << " sorting: " << s->GetSortingPolicy() << std::endl ;
+ }
+
+ }
+
+ is_initialised = true;
+
+ }
+
+ MarlinTrk::IMarlinTrack* MarlinDDKalTest::createTrack() {
+
+ if ( ! is_initialised ) {
+
+ std::stringstream errorMsg;
+ errorMsg << "MarlinDDKalTest::createTrack: Fitter not initialised. MarlinDDKalTest::init() must be called before MarlinDDKalTest::createTrack()" << std::endl ;
+ throw MarlinTrk::Exception(errorMsg.str());
+
+ }
+
+ return new MarlinDDKalTestTrack(this) ;
+
+ }
+
+ void MarlinDDKalTest::includeMultipleScattering( bool msOn ) {
+
+ streamlog_out( DEBUG2 ) << " **** MarlinDDKalTest::includeMultipleScattering( " << msOn << " ) called " << std::endl ;
+
+ if( msOn == true ) {
+ _det->SwitchOnMS();
+ }
+ else{
+ _det->SwitchOffMS();
+ }
+
+ }
+
+ void MarlinDDKalTest::includeEnergyLoss( bool energyLossOn ) {
+
+ streamlog_out( DEBUG2 ) << " **** MarlinDDKalTest::includeEnergyLoss( " << energyLossOn << " ) called " << std::endl ;
+
+ if( energyLossOn == true ) {
+ _det->SwitchOnDEDX();
+ }
+ else{
+ _det->SwitchOffDEDX();
+ }
+
+ }
+
+ void MarlinDDKalTest::getSensitiveMeasurementModulesForLayer( int layerID, std::vector< const DDVMeasLayer *>& measmodules) const {
+
+ if( ! measmodules.empty() ) {
+
+ std::stringstream errorMsg;
+ errorMsg << "MarlinDDKalTest::getSensitiveMeasurementModulesForLayer vector passed as second argument is not empty " << std::endl ;
+ throw MarlinTrk::Exception(errorMsg.str());
+
+ }
+
+ streamlog_out( DEBUG0 ) << "MarlinDDKalTest::getSensitiveMeasurementModulesForLayer: layerID = " << layerID << std::endl;
+
+ std::multimap<Int_t, const DDVMeasLayer *>::const_iterator it; //Iterator to be used along with ii
+
+
+
+ // for(it = _active_measurement_modules_by_layer.begin(); it != _active_measurement_modules_by_layer.end(); ++it) {
+ // streamlog_out( DEBUG0 ) << "Key = "<< ttdecodeILD(it->first) <<" Value = "<<it->second << std::endl ;
+ // }
+
+
+ std::pair<std::multimap<Int_t, const DDVMeasLayer *>::const_iterator, std::multimap<Int_t, const DDVMeasLayer *>::const_iterator> ii;
+
+ // set the module and sensor bit ranges to zero as these are not used in the map
+ lcio::BitField64 bf( UTIL::LCTrackerCellID::encoding_string() ) ;
+ bf.setValue( layerID ) ;
+ bf[lcio::LCTrackerCellID::module()] = 0 ;
+ bf[lcio::LCTrackerCellID::sensor()] = 0 ;
+ layerID = bf.lowWord();
+
+ ii = this->_active_measurement_modules_by_layer.equal_range(layerID); // set the first and last entry in ii;
+
+ for(it = ii.first; it != ii.second; ++it) {
+ // streamlog_out( DEBUG0 ) <<"Key = "<< it->first <<" Value = "<<it->second << std::endl ;
+ measmodules.push_back( it->second ) ;
+ }
+
+ }
+
+ void MarlinDDKalTest::getSensitiveMeasurementModules( int moduleID , std::vector< const DDVMeasLayer *>& measmodules ) const {
+
+ if( ! measmodules.empty() ) {
+
+ std::stringstream errorMsg;
+ errorMsg << "MarlinDDKalTest::getSensitiveMeasurementLayer vector passed as second argument is not empty " << std::endl ;
+ throw MarlinTrk::Exception(errorMsg.str());
+
+ }
+
+ std::pair<std::multimap<int, const DDVMeasLayer *>::const_iterator, std::multimap<Int_t, const DDVMeasLayer *>::const_iterator> ii;
+ ii = this->_active_measurement_modules.equal_range(moduleID); // set the first and last entry in ii;
+
+ std::multimap<int,const DDVMeasLayer *>::const_iterator it; //Iterator to be used along with ii
+
+
+ for(it = ii.first; it != ii.second; ++it) {
+ // std::cout<<"Key = "<<it->first<<" Value = "<<it->second << std::endl ;
+ measmodules.push_back( it->second ) ;
+ }
+ }
+
+
+ void MarlinDDKalTest::storeActiveMeasurementModuleIDs(TVKalDetector* detector) {
+
+ Int_t nLayers = detector->GetEntriesFast() ;
+
+ for( int i=0; i < nLayers; ++i ) {
+
+ const DDVMeasLayer* ml = dynamic_cast<const DDVMeasLayer*>( detector->At( i ) );
+
+ if( ! ml ) {
+ std::stringstream errorMsg;
+ errorMsg << "MarlinDDKalTest::storeActiveMeasurementLayerIDs dynamic_cast to DDVMeasLayer* failed " << std::endl ;
+ throw MarlinTrk::Exception(errorMsg.str());
+ }
+
+ if( ml->IsActive() ) {
+
+ // then get all the sensitive element id's assosiated with this DDVMeasLayer and store them in the map
+ std::vector<int>::const_iterator it = ml->getCellIDs().begin();
+
+ while ( it!=ml->getCellIDs().end() ) {
+
+ int sensitive_element_id = *it;
+ this->_active_measurement_modules.insert(std::pair<int,const DDVMeasLayer*>( sensitive_element_id, ml ));
+ ++it;
+
+ }
+
+ int subdet_layer_id = ml->getLayerID() ;
+
+ this->_active_measurement_modules_by_layer.insert(std::pair<int ,const DDVMeasLayer*>(subdet_layer_id,ml));
+
+ streamlog_out(DEBUG0) << "MarlinDDKalTest::storeActiveMeasurementLayerIDs added active layer with "
+ << " LayerID = " << subdet_layer_id << " and DetElementIDs " ;
+
+ for (it = ml->getCellIDs().begin(); it!=ml->getCellIDs().end(); ++it) {
+
+ streamlog_out(DEBUG0) << " : " << *it ;
+
+ }
+
+ streamlog_out(DEBUG0) << std::endl;
+
+
+
+
+ }
+
+ }
+
+ }
+
+ const DDVMeasLayer* MarlinDDKalTest::getLastMeasLayer(THelicalTrack const& hel, TVector3 const& point) const {
+
+ THelicalTrack helix = hel;
+
+ double deflection_to_point = 0 ;
+ helix.MoveTo( point, deflection_to_point , 0 , 0) ;
+
+ bool isfwd = ((helix.GetKappa() > 0 && deflection_to_point < 0) || (helix.GetKappa() <= 0 && deflection_to_point > 0)) ? true : false;
+
+ int mode = isfwd ? -1 : +1 ;
+
+ // streamlog_out( DEBUG4 ) << " MarlinDDKalTest - getLastMeasLayer deflection to point = " << deflection_to_point << " kappa = " << helix.GetKappa() << " mode = " << mode << std::endl ;
+ // streamlog_out( DEBUG4 ) << " Point to move to:" << std::endl;
+ // point.Print();
+
+ int nsufaces = _det->GetEntriesFast();
+
+ const TVSurface* ml_retval = nullptr;
+ double min_deflection = DBL_MAX;
+
+ for(int i=0; i<nsufaces; ++i) {
+
+ double defection_angle = 0 ;
+ TVector3 crossing_point ;
+
+ const TVSurface *sfp = static_cast<const TVSurface *>(_det->At(i)); // surface at destination
+
+ int does_cross = sfp->CalcXingPointWith(helix, crossing_point, defection_angle, mode) ;
+
+ if( does_cross ) {
+
+ const double deflection = fabs( deflection_to_point - defection_angle ) ;
+
+ if( deflection < min_deflection ) {
+
+ // streamlog_out( DEBUG4 ) << " MarlinDDKalTest - crossing found for suface = " << ml.GetMLName()
+ // << std::endl
+ // << " min_deflection = " << min_deflection
+ // << " deflection = " << deflection
+ // << " deflection angle = " << defection_angle
+ // << std::endl
+ // << " x = " << crossing_point.X()
+ // << " y = " << crossing_point.Y()
+ // << " z = " << crossing_point.Z()
+ // << " r = " << crossing_point.Perp()
+ // << std::endl ;
+
+ min_deflection = deflection ;
+ ml_retval = sfp;
+ }
+
+ }
+
+ }
+
+ return dynamic_cast<const DDVMeasLayer *>(ml_retval);
+ }
+
+ const DDVMeasLayer* MarlinDDKalTest::findMeasLayer( edm4hep::TrackerHit trkhit) const {
+
+ const TVector3 hit_pos( trkhit.getPosition()[0], trkhit.getPosition()[1], trkhit.getPosition()[2]) ;
+
+ return this->findMeasLayer( trkhit.getCellID(), hit_pos ) ;
+
+ }
+
+ const DDVMeasLayer* MarlinDDKalTest::findMeasLayer( int detElementID, const TVector3& point) const {
+
+ const DDVMeasLayer* ml = 0; // return value
+
+ std::vector<const DDVMeasLayer*> meas_modules ;
+
+ // search for the list of measurement layers associated with this CellID
+ this->getSensitiveMeasurementModules( detElementID, meas_modules ) ;
+
+ if( meas_modules.size() == 0 ) { // no measurement layers found
+
+ UTIL::BitField64 encoder( UTIL::LCTrackerCellID::encoding_string() ) ;
+ encoder.setValue(detElementID) ;
+
+ std::stringstream errorMsg;
+ errorMsg << "MarlinDDKalTest::findMeasLayer module id unkown: moduleID = " << detElementID
+ << " [" << encoder.valueString() << "]" << std::endl ;
+ throw MarlinTrk::Exception(errorMsg.str());
+
+ }
+ else if (meas_modules.size() == 1) { // one to one mapping
+
+ ml = meas_modules[0] ;
+
+ }
+ else { // layer has been split
+
+ bool surf_found(false);
+ //std::cout << "MarlinDDKalTest::findMeasLayer module size " << meas_modules.size() << std::endl;
+ // loop over the measurement layers associated with this CellID and find the correct one using the position of the hit
+ for( unsigned int i=0; i < meas_modules.size(); ++i) {
+
+ const TVSurface* surf = 0;
+ //std::cout << meas_modules[i]->surface()->id() << " " << meas_modules[i]->surface()->type() << std::endl;
+ if( ! (surf = dynamic_cast<const TVSurface*> ( meas_modules[i] )) ) {
+ std::stringstream errorMsg;
+ errorMsg << "MarlinDDKalTest::findMeasLayer dynamic_cast failed for surface type: moduleID = " << detElementID << std::endl ;
+ throw MarlinTrk::Exception(errorMsg.str());
+ }
+
+ bool hit_on_surface = surf->IsOnSurface(point);
+
+ if( (!surf_found) && hit_on_surface ){
+
+ ml = meas_modules[i] ;
+ surf_found = true ;
+
+ }
+ else if( surf_found && hit_on_surface ) { // only one surface should be found, if not throw
+
+ std::stringstream errorMsg;
+ errorMsg << "MarlinDDKalTest::findMeasLayer point found to be on two surfaces: moduleID = " << detElementID << std::endl ;
+ throw MarlinTrk::Exception(errorMsg.str());
+ }
+
+ }
+ if( ! surf_found ){ // print out debug info
+ streamlog_out(DEBUG1) << "MarlinDDKalTest::findMeasLayer point not found to be on any surface matching moduleID = "
+ << detElementID
+ << ": x = " << point.x()
+ << " y = " << point.y()
+ << " z = " << point.z()
+ << std::endl ;
+ }
+ else{
+ streamlog_out(DEBUG1) << "MarlinDDKalTest::findMeasLayer point found to be on surface matching moduleID = "
+ << detElementID
+ << ": x = " << point.x()
+ << " y = " << point.y()
+ << " z = " << point.z()
+ << std::endl ;
+ }
+ }
+
+ return ml ;
+
+ }
+
+} // end of namespace MarlinTrk
diff --git a/Service/TrackSystemSvc/src/MarlinDDKalTest.h b/Service/TrackSystemSvc/src/MarlinDDKalTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d064966a2000a0ca7d897bb0ac889614f7f2a9b
--- /dev/null
+++ b/Service/TrackSystemSvc/src/MarlinDDKalTest.h
@@ -0,0 +1,139 @@
+#ifndef MarlinDDKalTest_h
+#define MarlinDDKalTest_h
+
+#include "TrackSystemSvc/IMarlinTrkSystem.h"
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+#include "DiagnosticsController.h"
+#endif
+
+#include "DDKalTest/DDKalDetector.h"
+
+//LCIO:
+#include "lcio.h"
+#include "UTIL/BitField64.h"
+//#include "UTIL/LCTOOLS.h"
+//#include <LCRTRelations.h>
+
+#include "streamlog/streamlog.h"
+
+#include "TObjArray.h"
+#include "TVector3.h"
+
+#include <cmath>
+#include <vector>
+
+
+class TKalDetCradle ;
+class TVKalDetector ;
+class DDVMeasLayer ;
+class THelicalTrack ;
+
+class DDCylinderMeasLayer;
+
+namespace edm4hep{
+ class TrackerHit ;
+}
+
+namespace MarlinTrk{
+
+ /** Interface to KaltTest Kalman fitter - instantiates and holds the detector geometry.
+ */
+ class MarlinDDKalTest : public MarlinTrk::IMarlinTrkSystem {
+
+ public:
+
+ friend class MarlinDDKalTestTrack;
+
+ // define some configuration constants
+ static const bool FitBackward = kIterBackward ;
+ static const bool FitForward = kIterForward ;
+ static const bool OrderOutgoing = true ;
+ static const bool OrderIncoming = false ;
+
+ /// Default c'tor
+ MarlinDDKalTest() ;
+ MarlinDDKalTest(const MarlinDDKalTest&) = delete;
+ MarlinDDKalTest const& operator=(const MarlinDDKalTest&) = delete;
+
+ /** d'tor */
+ ~MarlinDDKalTest() ;
+
+ /** Sets the specified option ( one of the constants defined in IMarlinTrkSystem::CFG )
+ * to the given value. Override here to re-configure E-loss and QMS
+ * after the initialization.
+ */
+ virtual void setOption(unsigned CFGOption, bool val) ;
+
+ /** initialise track fitter system */
+ void init() ;
+
+ /// the name of the implementation
+ virtual std::string name() { return "DDKalTest" ; }
+
+ /** instantiate its implementation of the IMarlinTrack */
+ IMarlinTrack* createTrack();
+
+
+ protected:
+
+ /** take multiple scattering into account during the fit */
+ void includeMultipleScattering( bool on ) ;
+
+ /** take energy loss into account during the fit */
+ void includeEnergyLoss( bool on ) ;
+
+ /** Store active measurement module IDs for a given TVKalDetector needed for navigation */
+ void storeActiveMeasurementModuleIDs(TVKalDetector* detector);
+
+ /** Store active measurement module IDs needed for navigation */
+ void getSensitiveMeasurementModules( int detElementID, std::vector< const DDVMeasLayer *>& measmodules) const;
+
+ /** Store active measurement module IDs needed for navigation */
+ void getSensitiveMeasurementModulesForLayer( int layerID, std::vector<const DDVMeasLayer *>& measmodules) const;
+
+ // void init(bool MSOn, bool EnergyLossOn) ;
+ bool is_initialised=false;
+
+ //** find the measurment layer for a given hit
+ const DDVMeasLayer* findMeasLayer( edm4hep::TrackerHit trkhit) const ;
+ //** find the measurment layer for a given det element ID and point in space
+ const DDVMeasLayer* findMeasLayer( int detElementID, const TVector3& point) const ;
+
+ // get the last layer crossed by the helix when extrapolating from the present position to the pca to point
+ const DDVMeasLayer* getLastMeasLayer(THelicalTrack const& helix, TVector3 const& point) const ;
+
+ const DDCylinderMeasLayer* getIPLayer() const { return _ipLayer; }
+
+
+ // members:
+
+ const DDCylinderMeasLayer* _ipLayer=nullptr;
+
+ TKalDetCradle* _det=nullptr; // the detector cradle
+
+ std::multimap< int,const DDVMeasLayer *> _active_measurement_modules{};
+
+ std::multimap< int,const DDVMeasLayer *> _active_measurement_modules_by_layer{};
+
+ std::vector< DDKalDetector* > _detectors{};
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+
+ private:
+ MarlinTrk::DiagnosticsController _diagnostics;
+
+ public:
+
+ /** Return the pointer to the Diagnositics Object. Forseen for internal diagnostics, only available when complied with MARLINTRK_DIAGNOSTICS_ON defined.
+ */
+ virtual void * getDiagnositicsPointer() { return &_diagnostics ; }
+
+#endif
+
+
+ } ;
+
+} // end of namespace MarlinTrk
+
+#endif
diff --git a/Service/TrackSystemSvc/src/MarlinDDKalTestTrack.cc b/Service/TrackSystemSvc/src/MarlinDDKalTestTrack.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b98a9eb498a00fa90a66d0143f8d248eed0a58e4
--- /dev/null
+++ b/Service/TrackSystemSvc/src/MarlinDDKalTestTrack.cc
@@ -0,0 +1,1724 @@
+#include "MarlinDDKalTestTrack.h"
+
+#include "MarlinDDKalTest.h"
+#include "TrackSystemSvc/IMarlinTrkSystem.h"
+
+#include <kaltest/TKalDetCradle.h>
+#include <kaltest/TKalTrack.h>
+#include <kaltest/TKalTrackState.h>
+#include "kaltest/TKalTrackSite.h"
+#include "TKalFilterCond.h"
+
+//#include <lcio.h>
+#include "edm4hep/TrackerHit.h"
+//#include <EVENT/TrackerHitPlane.h>
+
+#include <UTIL/BitField64.h>
+#include <UTIL/Operators.h>
+#include "UTIL/LCTrackerConf.h"
+
+#include "DDKalTest/DDCylinderMeasLayer.h" // needed for dedicated IP Layer
+#include "DDKalTest/DDCylinderHit.h"
+#include "DDKalTest/DDPlanarHit.h"
+//#include "DDKalTest/DDPlanarStripHit.h"
+
+#include <sstream>
+
+#include "streamlog/streamlog.h"
+
+#include "podio/podioVersion.h"
+
+/** Helper class for defining a filter condition based on the delta chi2 in the AddAndFilter step.
+ */
+class KalTrackFilter : public TKalFilterCond{
+
+public:
+
+ /** C'tor - takes as optional argument the maximum allowed delta chi2 for adding the hit (in IsAccepted() )
+ */
+ KalTrackFilter(double maxDeltaChi2 = DBL_MAX) : _maxDeltaChi2( maxDeltaChi2 ), _passed_last_filter_step(true) {
+ }
+ virtual ~KalTrackFilter() {}
+
+ virtual Bool_t IsAccepted(const TKalTrackSite &site) {
+
+ double deltaChi2 = site.GetDeltaChi2();
+
+ streamlog_out( DEBUG1 ) << " KalTrackFilter::IsAccepted called ! deltaChi2 = " << std::scientific << deltaChi2 << " _maxDeltaChi2 = " << _maxDeltaChi2 << std::endl;
+
+ _passed_last_filter_step = deltaChi2 < _maxDeltaChi2;
+
+ return ( _passed_last_filter_step ) ;
+ }
+
+ void resetFilterStatus() { _passed_last_filter_step = true; }
+ bool passedLastFilterStep() const { return _passed_last_filter_step; }
+
+protected:
+
+ double _maxDeltaChi2 ;
+ bool _passed_last_filter_step;
+
+} ;
+
+namespace MarlinTrk {
+
+ //---------------------------------------------------------------------------------------------------------------
+
+ std::string cellIDString( int detElementID ) {
+ lcio::BitField64 bf( UTIL::LCTrackerCellID::encoding_string() ) ;
+ bf.setValue( detElementID ) ;
+ return bf.valueString() ;
+ }
+
+ //---------------------------------------------------------------------------------------------------------------
+
+
+ MarlinDDKalTestTrack::MarlinDDKalTestTrack( MarlinDDKalTest* ktest)
+ : _ktest(ktest),
+#if PODIO_BUILD_VERSION < PODIO_VERSION(0, 17, 4)
+ _trackHitAtPositiveNDF(edm4hep::TrackerHit(0)) {
+#else
+ _trackHitAtPositiveNDF(edm4hep::TrackerHit::makeEmpty()) {
+#endif
+
+ _kaltrack = new TKalTrack() ;
+ _kaltrack->SetOwner() ;
+
+ _kalhits = new TObjArray() ;
+ _kalhits->SetOwner() ;
+
+ _initialised = false ;
+ _fitDirection = false ;
+ _smoothed = false ;
+
+ _hitIndexAtPositiveNDF = 0;
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+ _ktest->_diagnostics.new_track(this) ;
+#endif
+
+ //streamlog::out.init(std::cout, "MarlinDDKalTestTrack");
+ //streamlog::out.addLevelName<streamlog::DEBUG1>();
+ }
+
+
+ MarlinDDKalTestTrack::~MarlinDDKalTestTrack(){
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+ _ktest->_diagnostics.end_track() ;
+#endif
+
+ delete _kaltrack ;
+ delete _kalhits ;
+ }
+
+ void MarlinDDKalTestTrack::setMass(double mass) { _kaltrack->SetMass( mass ) ; }
+
+ double MarlinDDKalTestTrack::getMass() { return _kaltrack->GetMass() ; }
+
+
+ int MarlinDDKalTestTrack::addHit(edm4hep::TrackerHit& trkhit) {
+ streamlog_out(DEBUG1) << "MarlinDDKalTestTrack::addHit: trkhit = " << trkhit.id() << std::endl;
+ return this->addHit( trkhit, _ktest->findMeasLayer( trkhit )) ;
+
+ }
+
+ int MarlinDDKalTestTrack::addHit(edm4hep::TrackerHit& trkhit, const DDVMeasLayer* ml) {
+
+ streamlog_out(DEBUG1) << "MarlinDDKalTestTrack::addHit: trkhit = " << trkhit.id() << " pos: " << trkhit.getPosition() << " ml = " << ml << " " << ml->GetName() << std::endl ;
+
+ if( trkhit.isAvailable() && ml ) {
+ return this->addHit( trkhit, ml->ConvertLCIOTrkHit(trkhit), ml) ;
+ }
+ else {
+ streamlog_out( ERROR ) << " MarlinDDKalTestTrack::addHit - bad inputs " << trkhit << " ml : " << ml << std::endl ;
+ return bad_intputs ;
+ }
+
+ }
+
+ int MarlinDDKalTestTrack::addHit( edm4hep::TrackerHit& trkhit, DDVTrackHit* kalhit, const DDVMeasLayer* ml) {
+
+ if( kalhit && ml ) {
+ _kalhits->Add(kalhit ) ; // Add hit and set surface found
+ _lcio_hits_to_kaltest_hits[trkhit] = kalhit ; // add hit to map relating lcio and kaltest hits
+ // _kaltest_hits_to_lcio_hits[kalhit] = trkhit ; // add hit to map relating kaltest and lcio hits
+ }
+ else {
+ delete kalhit;
+ return bad_intputs ;
+ }
+
+ streamlog_out(DEBUG1) << "MarlinDDKalTestTrack::addHit: hit added "
+ << "number of hits for track = " << _kalhits->GetEntries()
+ << std::endl ;
+
+ return success ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::initialise( bool fitDirection ) {;
+
+
+ //SJA:FIXME: check here if the track is already initialised, and for now don't allow it to be re-initialised
+ // if the track is going to be re-initialised then we would need to do it directly on the first site
+ if ( _initialised ) {
+
+ throw MarlinTrk::Exception("Track fit already initialised");
+
+ }
+
+ if (_kalhits->GetEntries() < 3) {
+
+ streamlog_out( ERROR) << "<<<<<< MarlinDDKalTestTrack::initialise: Shortage of Hits! nhits = "
+ << _kalhits->GetEntries() << " >>>>>>>" << std::endl;
+ return error ;
+
+ }
+
+ _fitDirection = fitDirection ;
+
+ // establish the hit order
+ Int_t i1, i2, i3; // (i1,i2,i3) = (1st,mid,last) hit to filter
+ if (_fitDirection == kIterBackward) {
+ i3 = 0 ; // fg: first index is 0 and not 1
+ i1 = _kalhits->GetEntries() - 1;
+ i2 = i1 / 2;
+ } else {
+ i1 = 0 ;
+ i3 = _kalhits->GetEntries() - 1;
+ i2 = i3 / 2;
+ }
+
+
+
+ TVTrackHit *startingHit = dynamic_cast<TVTrackHit *>(_kalhits->At(i1));
+
+ // ---------------------------
+ // Create an initial start site for the track using the first hit
+ // ---------------------------
+ // set up a dummy hit needed to create initial site
+
+ TVTrackHit* pDummyHit = 0;
+
+ if ( (pDummyHit = dynamic_cast<DDCylinderHit *>( startingHit )) ) {
+ pDummyHit = (new DDCylinderHit(*static_cast<DDCylinderHit*>( startingHit )));
+ }
+ else if ( (pDummyHit = dynamic_cast<DDPlanarHit *>( startingHit )) ) {
+ pDummyHit = (new DDPlanarHit(*static_cast<DDPlanarHit*>( startingHit )));
+ }
+ // else if ( DDPlanarStripHit_DIM == 2 && (pDummyHit = dynamic_cast<DDPlanarStripHit *>( startingHit )) ) {
+ // pDummyHit = (new DDPlanarStripHit(*static_cast<DDPlanarStripHit*>( startingHit )));
+ // }
+ else {
+ streamlog_out( ERROR) << "<<<<<<<<< MarlinDDKalTestTrack::initialise: dynamic_cast failed for hit type >>>>>>>" << std::endl;
+ return error ;
+ }
+
+ TVTrackHit& dummyHit = *pDummyHit;
+
+ //SJA:FIXME: this constants should go in a header file
+ // give the dummy hit huge errors so that it does not contribute to the fit
+ dummyHit(0,1) = 1.e16; // give a huge error to d
+ dummyHit(1,1) = 1.e16; // give a huge error to z
+
+ // use dummy hit to create initial site
+ TKalTrackSite& initialSite = *new TKalTrackSite(dummyHit);
+
+ initialSite.SetHitOwner();// site owns hit
+ initialSite.SetOwner(); // site owns states
+
+ // ---------------------------
+ // Create initial helix
+ // ---------------------------
+
+ TVTrackHit &h1 = *dynamic_cast<TVTrackHit *>(_kalhits->At(i1)); // first hit
+ TVTrackHit &h2 = *dynamic_cast<TVTrackHit *>(_kalhits->At(i2)); // middle hit
+ TVTrackHit &h3 = *dynamic_cast<TVTrackHit *>(_kalhits->At(i3)); // last hit
+ TVector3 x1 = h1.GetMeasLayer().HitToXv(h1);
+ TVector3 x2 = h2.GetMeasLayer().HitToXv(h2);
+ TVector3 x3 = h3.GetMeasLayer().HitToXv(h3);
+
+ if ( h1.GetDimension() == 1 || h2.GetDimension() == 1 || h3.GetDimension() == 1 ) {
+
+ throw MarlinTrk::Exception("Track fit cannot be initialised from 1 Dimensional hits. Use method MarlinDDKalTestTrack::initialise( bool fitDirection )");
+
+ }
+
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::initialise Create initial helix from hits: \n "
+ << "P1 x = " << x1.x() << " y = " << x1.y() << " z = " << x1.z() << " r = " << x1.Perp() << "\n "
+ << "P2 x = " << x2.x() << " y = " << x2.y() << " z = " << x2.z() << " r = " << x2.Perp() << "\n "
+ << "P3 x = " << x3.x() << " y = " << x3.y() << " z = " << x3.z() << " r = " << x3.Perp() << "\n "
+ << "Bz = " << h1.GetBfield() << " direction = " << _fitDirection
+ << std::endl;
+
+ // create helix using 3 global space points
+ THelicalTrack helstart(x1, x2, x3, h1.GetBfield(), _fitDirection); // initial helix
+
+ // ---------------------------
+ // Set up initial track state ... could try to use lcio track parameters ...
+ // ---------------------------
+
+ static TKalMatrix initialState(kSdim,1) ;
+ initialState(0,0) = 0.0 ; // dr
+ initialState(1,0) = helstart.GetPhi0() ; // phi0
+ initialState(2,0) = helstart.GetKappa() ; // kappa
+ initialState(3,0) = 0.0 ; // dz
+ initialState(4,0) = helstart.GetTanLambda() ; // tan(lambda)
+ if (kSdim == 6) initialState(5,0) = 0.; // t0
+
+
+ // ---------------------------
+ // Set up initial Covariance Matrix with very large errors
+ // ---------------------------
+
+ TKalMatrix Cov(kSdim,kSdim);
+
+ // make sure everything is initialised to zero
+ for (int i=0; i<kSdim*kSdim; ++i) {
+ Cov.GetMatrixArray()[i] = 0.0;
+ }
+
+// for (Int_t i=0; i<kSdim; i++) {
+// // fg: if the error is too large the initial helix parameters might be changed extremely by the first three (or so) hits,
+// // such that the fit will not work because the helix curls away and does not hit the next layer !!!
+// Cov(i,i) = 1.e2 ; // initialise diagonal elements of dummy error matrix
+// }
+
+ // prefer translation over rotation of the trackstate early in the fit
+
+ Cov(0,0) = 1.e6 ; // d0
+ Cov(1,1) = 1.e2 ; // dphi0
+ Cov(2,2) = 1.e1 ; // dkappa
+ Cov(3,3) = 1.e6 ; // dz
+ Cov(4,4) = 1.e1 ; // dtanL
+ if (kSdim == 6) Cov(5,5) = 1.e2; // t0
+
+ // Add initial states to the site
+ initialSite.Add(new TKalTrackState(initialState,Cov,initialSite,TVKalSite::kPredicted));
+ initialSite.Add(new TKalTrackState(initialState,Cov,initialSite,TVKalSite::kFiltered));
+
+ // add the initial site to the track: that is, give the track initial parameters and covariance
+ // matrix at the starting measurement layer
+ _kaltrack->Add(&initialSite);
+
+ _initialised = true ;
+
+ streamlog_out( DEBUG2 ) << " track parameters used for init : " << std::scientific << std::setprecision(6)
+ << "\t D0 " << 0.0
+ << "\t Phi :" << toBaseRange( helstart.GetPhi0() + M_PI/2. )
+ << "\t Omega " << 1. /helstart.GetRho()
+ << "\t Z0 " << 0.0
+ << "\t tan(Lambda) " << helstart.GetTanLambda()
+
+ << "\t pivot : [" << helstart.GetPivot().X() << ", " << helstart.GetPivot().Y() << ", " << helstart.GetPivot().Z()
+ << " - r: " << std::sqrt( helstart.GetPivot().X()*helstart.GetPivot().X()+helstart.GetPivot().Y()*helstart.GetPivot().Y() ) << "]"
+ << std::endl ;
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+
+ // convert to LICO parameters first
+
+ double d0 = 0.0 ;
+ double phi = toBaseRange( helstart.GetPhi0() + M_PI/2. );
+ double omega = 1. /helstart.GetRho() ;
+ double z0 = 0.0 ;
+ double tanLambda = helstart.GetTanLambda() ;
+
+// Cov.Print();
+
+ _ktest->_diagnostics.set_intial_track_parameters(d0,
+ phi,
+ omega,
+ z0,
+ tanLambda,
+ helstart.GetPivot().X(),
+ helstart.GetPivot().Y(),
+ helstart.GetPivot().Z(),
+ Cov);
+
+
+
+#endif
+
+
+ return success ;
+
+ }
+
+ int MarlinDDKalTestTrack::initialise(const edm4hep::TrackState& ts, double /*bfield_z*/, bool fitDirection ) {
+
+ // the bfield_z is not taken from the argument but from the first hit
+ // should consider changing the interface ...
+
+ if (_kalhits->GetEntries() == 0) {
+
+ streamlog_out( ERROR) << "<<<<<< MarlinDDKalTestTrack::Initialise: Number of Hits is Zero. Cannot Initialise >>>>>>>" << std::endl;
+ return error ;
+
+ }
+
+ //SJA:FIXME: check here if the track is already initialised, and for now don't allow it to be re-initialised
+ // if the track is going to be re-initialised then we would need to do it directly on the first site
+ if ( _initialised ) {
+
+ throw MarlinTrk::Exception("Track fit already initialised");
+
+ }
+
+ streamlog_out( DEBUG2 ) << "MarlinDDKalTestTrack::initialise using TrackState: track parameters used for init : " << ts << std::endl;
+
+
+ _fitDirection = fitDirection ;
+
+ // get Bz from first hit
+ TVTrackHit &h1 = *dynamic_cast<TVTrackHit *>(_kalhits->At(0));
+ double Bz = h1.GetBfield() ;
+
+ // for GeV, Tesla, R in mm
+ double alpha = Bz * 2.99792458E-4 ;
+
+ double kappa = ( Bz == 0.0 ? DBL_MAX : ts.omega / alpha ) ;
+
+ THelicalTrack helix( -ts.D0,
+ toBaseRange( ts.phi - M_PI/2. ),
+ kappa,
+ ts.Z0,
+ ts.tanLambda,
+ ts.referencePoint[0],
+ ts.referencePoint[1],
+ ts.referencePoint[2],
+ Bz );
+
+ TMatrixD cov(5,5) ;
+
+ auto covLCIO = ts.covMatrix;
+
+ cov( 0 , 0 ) = covLCIO[ 0] ; // d0, d0
+ cov( 0 , 1 ) = - covLCIO[ 1] ; // d0, phi
+ cov( 0 , 2 ) = - covLCIO[ 3] / alpha ; // d0, kappa
+ cov( 0 , 3 ) = - covLCIO[ 6] ; // d0, z0
+ cov( 0 , 4 ) = - covLCIO[10] ; // d0, tanl
+
+ cov( 1 , 0 ) = - covLCIO[ 1] ; // phi, d0
+ cov( 1 , 1 ) = covLCIO[ 2] ; // phi, phi
+ cov( 1 , 2 ) = covLCIO[ 4] / alpha ; // phi, kappa
+ cov( 1 , 3 ) = covLCIO[ 7] ; // phi, z0
+ cov( 1 , 4 ) = covLCIO[11] ; // tanl, phi
+
+ cov( 2 , 0 ) = - covLCIO[ 3] / alpha ; // kappa, d0
+ cov( 2 , 1 ) = covLCIO[ 4] / alpha ; // kappa, phi
+ cov( 2 , 2 ) = covLCIO[ 5] / (alpha * alpha) ; // kappa, kappa
+ cov( 2 , 3 ) = covLCIO[ 8] / alpha ; // kappa, z0
+ cov( 2 , 4 ) = covLCIO[12] / alpha ; // kappa, tanl
+
+ cov( 3 , 0 ) = - covLCIO[ 6] ; // z0, d0
+ cov( 3 , 1 ) = covLCIO[ 7] ; // z0, phi
+ cov( 3 , 2 ) = covLCIO[ 8] / alpha ; // z0, kappa
+ cov( 3 , 3 ) = covLCIO[ 9] ; // z0, z0
+ cov( 3 , 4 ) = covLCIO[13] ; // z0, tanl
+
+ cov( 4 , 0 ) = - covLCIO[10] ; // tanl, d0
+ cov( 4 , 1 ) = covLCIO[11] ; // tanl, phi
+ cov( 4 , 2 ) = covLCIO[12] / alpha ; // tanl, kappa
+ cov( 4 , 3 ) = covLCIO[13] ; // tanl, z0
+ cov( 4 , 4 ) = covLCIO[14] ; // tanl, tanl
+
+// cov.Print();
+
+ // move the helix to either the position of the last hit or the first depending on initalise_at_end
+
+ // default case initalise_at_end
+ int index = _kalhits->GetEntries() - 1 ;
+ // or initialise at start
+ if( _fitDirection == IMarlinTrack::forward ){
+ index = 0 ;
+ }
+
+ TVTrackHit* kalhit = dynamic_cast<TVTrackHit *>(_kalhits->At(index));
+
+ double dphi;
+
+ TVector3 initial_pivot ;
+
+ // Leave the pivot at the origin for a 1-dim hit
+ if (kalhit->GetDimension() > 1) {
+
+ initial_pivot = kalhit->GetMeasLayer().HitToXv(*kalhit) ;
+ }
+ else{
+ initial_pivot = TVector3(0.0,0.0,0.0);
+ }
+
+
+ // ---------------------------
+ // Create an initial start site for the track using the hit
+ // ---------------------------
+ // set up a dummy hit needed to create initial site
+
+ TVTrackHit* pDummyHit = 0;
+
+ if ( (pDummyHit = dynamic_cast<DDCylinderHit *>( kalhit )) ) {
+ pDummyHit = (new DDCylinderHit(*static_cast<DDCylinderHit*>( kalhit )));
+
+ }
+ else if ( (pDummyHit = dynamic_cast<DDPlanarHit *>( kalhit )) ) {
+ pDummyHit = (new DDPlanarHit(*static_cast<DDPlanarHit*>( kalhit )));
+ //}
+ //else if ( (pDummyHit = dynamic_cast<DDPlanarStripHit *>( kalhit )) ) {
+ //pDummyHit = (new DDPlanarStripHit(*static_cast<DDPlanarStripHit*>( kalhit )));
+ if( pDummyHit->GetDimension() == 1 ) {
+
+ const TVMeasLayer *ml = &pDummyHit->GetMeasLayer();
+
+ const TVSurface* surf = dynamic_cast<const TVSurface*>(ml);
+
+ if (surf) {
+ double phi;
+
+ surf->CalcXingPointWith(helix, initial_pivot, phi);
+ streamlog_out( DEBUG ) << " MarlinDDKalTestTrack::initialise - CalcXingPointWith called for 1d hit ... " << std::endl ;
+
+ } else {
+ streamlog_out( ERROR) << "<<<<<<<<< MarlinDDKalTestTrack::initialise: dynamic_cast failed for TVSurface >>>>>>>" << std::endl;
+ return error ;
+ }
+
+ }
+ }
+ else {
+ streamlog_out( ERROR) << "<<<<<<<<< MarlinDDKalTestTrack::initialise: dynamic_cast failed for hit type >>>>>>>" << std::endl;
+ return error ;
+ }
+
+ TVTrackHit& dummyHit = *pDummyHit;
+
+ //SJA:FIXME: this constants should go in a header file
+ // give the dummy hit huge errors so that it does not contribute to the fit
+ dummyHit(0,1) = 1.e16; // give a huge error to d
+
+ if(dummyHit.GetDimension()>1) dummyHit(1,1) = 1.e16; // give a huge error to z
+
+ // use dummy hit to create initial site
+ TKalTrackSite& initialSite = *new TKalTrackSite(dummyHit);
+
+ initialSite.SetHitOwner();// site owns hit
+ initialSite.SetOwner(); // site owns states
+
+ // ---------------------------
+ // Set up initial track state
+ // ---------------------------
+
+ helix.MoveTo( initial_pivot, dphi, 0, &cov );
+
+ static TKalMatrix initialState(kSdim,1) ;
+ initialState(0,0) = helix.GetDrho() ; // d0
+ initialState(1,0) = helix.GetPhi0() ; // phi0
+ initialState(2,0) = helix.GetKappa() ; // kappa
+ initialState(3,0) = helix.GetDz(); // dz
+ initialState(4,0) = helix.GetTanLambda() ; // tan(lambda)
+ if (kSdim == 6) initialState(5,0) = 0.; // t0
+
+ // make sure that the pivot is in the right place
+ initialSite.SetPivot(initial_pivot);
+
+ // ---------------------------
+ // Set up initial Covariance Matrix
+ // ---------------------------
+
+ TKalMatrix covK(kSdim,kSdim) ;
+ for(int i=0;i<5;++i) {
+ for(int j=0;j<5;++j) {
+ covK[i][j] = cov[i][j] ;
+ }
+ }
+ if (kSdim == 6) covK(5,5) = 1.e6; // t0
+
+// covK.Print();
+
+ // Add initial states to the site
+ initialSite.Add(new TKalTrackState(initialState,covK,initialSite,TVKalSite::kPredicted));
+ initialSite.Add(new TKalTrackState(initialState,covK,initialSite,TVKalSite::kFiltered));
+
+ // add the initial site to the track: that is, give the track initial parameters and covariance
+ // matrix at the starting measurement layer
+ _kaltrack->Add(&initialSite);
+
+ _initialised = true ;
+
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+
+ // convert to LICO parameters first
+
+ double d0 = - helix.GetDrho() ;
+ double phi = toBaseRange( helix.GetPhi0() + M_PI/2. );
+ double omega = 1. /helix.GetRho() ;
+ double z0 = helix.GetDz() ;
+ double tanLambda = helix.GetTanLambda() ;
+
+ _ktest->_diagnostics.set_intial_track_parameters(d0,
+ phi,
+ omega,
+ z0,
+ tanLambda,
+ helix.GetPivot().X(),
+ helix.GetPivot().Y(),
+ helix.GetPivot().Z(),
+ covK);
+
+
+
+#endif
+
+ return success ;
+
+ }
+
+ int MarlinDDKalTestTrack::addAndFit( DDVTrackHit* kalhit, double& chi2increment, TKalTrackSite*& site, double maxChi2Increment) {
+
+ streamlog_out(DEBUG1) << "MarlinDDKalTestTrack::addAndFit called : maxChi2Increment = " << std::scientific << maxChi2Increment << std::endl ;
+
+ if ( ! _initialised ) {
+
+ throw MarlinTrk::Exception("Track fit not initialised");
+
+ }
+
+ // here do dynamic cast repeatedly in DEBUG statement as this will be stripped out any way for production code
+ // otherwise we have to do the cast outside of the DEBUG statement and it won't be stripped out
+ streamlog_out( DEBUG1 ) << "Kaltrack::addAndFit : add site to track at index : "
+ << (dynamic_cast<const DDVMeasLayer*>( &(kalhit->GetMeasLayer() ) ))->GetIndex()
+ << " for type "
+ << dynamic_cast<const DDVMeasLayer*>( &(kalhit->GetMeasLayer() ) )->GetName() ;
+ streamlog_out( DEBUG0 ) << " with CellIDs:";
+
+ for (unsigned int i = 0; i < (dynamic_cast<const DDVMeasLayer*>( &(kalhit->GetMeasLayer() ) )->getNCellIDs());++i) {
+ streamlog_out( DEBUG0 ) << " : "
+ << dynamic_cast<const DDVMeasLayer*>( &(kalhit->GetMeasLayer() ) )->getCellIDs()[i] ;
+
+ }
+
+ streamlog_out( DEBUG1 ) << std::endl ;
+
+ TKalTrackSite* temp_site = new TKalTrackSite(*kalhit); // create new site for this hit
+
+ KalTrackFilter filter( maxChi2Increment );
+ filter.resetFilterStatus();
+
+ temp_site->SetFilterCond( &filter ) ;
+
+
+ // this is the only point at which a hit is actually filtered
+ // and it is here that we can get the GetDeltaChi2 vs the maxChi2Increment
+ // it will always be possible to get the delta chi2 so long as we have a link to the sites ...
+ // although calling smooth will natrually update delta chi2.
+
+
+ if (!_kaltrack->AddAndFilter(*temp_site)) {
+
+ chi2increment = temp_site->GetDeltaChi2() ;
+ // get the measurement layer of the current hit
+ const DDVMeasLayer* ml = dynamic_cast<const DDVMeasLayer*>( &(kalhit->GetMeasLayer() ) ) ;
+ TVector3 pos = ml->HitToXv(*kalhit);
+ streamlog_out( DEBUG2 ) << "Kaltrack::addAndFit : site discarded! at index : " << ml->GetIndex()
+ << " for type " << ml->GetName()
+ << " chi2increment = " << chi2increment
+ << " maxChi2Increment = " << maxChi2Increment
+ << " x = " << pos.x()
+ << " y = " << pos.y()
+ << " z = " << pos.z()
+ << " with CellIDs: " << std::endl;
+ std::cout << "Kaltrack::addAndFit : site discarded! at index : " << ml->GetIndex()
+ << " for type " << ml->GetName()
+ << " chi2increment = " << chi2increment
+ << " maxChi2Increment = " << maxChi2Increment
+ << " x = " << pos.x()
+ << " y = " << pos.y()
+ << " z = " << pos.z()
+ << " with CellIDs: " << std::endl;
+
+ for (unsigned int i = 0; i < (dynamic_cast<const DDVMeasLayer*>( &(kalhit->GetMeasLayer() ) )->getNCellIDs());++i) {
+ streamlog_out( DEBUG1 ) << " CellID = "
+ << dynamic_cast<const DDVMeasLayer*>( &(kalhit->GetMeasLayer() ) )->getCellIDs()[i]
+ << std::endl ;
+ std::cout << " CellID = "
+ << dynamic_cast<const DDVMeasLayer*>( &(kalhit->GetMeasLayer() ) )->getCellIDs()[i]
+ << std::endl ;
+ }
+
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+ _ktest->_diagnostics.record_rejected_site(kalhit, temp_site);
+#endif
+
+ delete temp_site; // delete site if filter step failed
+
+
+ // compiling the code below with the cmake option CMAKE_BUILD_TYPE=Debug
+ // and with LDFLAGS=-Wl,--no-undefined
+ // causes an undefined reference error
+ // the problem gets fixed using the if/else statement below
+
+ // this fails
+ //return filter.usedForLastFilterStep() ? site_fails_chi2_cut : site_discarded ;
+
+ // this also fails..
+ //bool rc = filter.usedForLastFilterStep() ;
+ //return (rc ? site_fails_chi2_cut : site_discarded);
+
+ // but this works ?!!
+ //return ( true ? site_fails_chi2_cut : site_discarded);
+
+ // and this also works..
+ streamlog_out(DEBUG2) << " addAndFit : Site passed Chi2 filter step ? " << filter.passedLastFilterStep() << std::endl;
+
+ if( filter.passedLastFilterStep() == false ) {
+ return site_fails_chi2_cut ;
+ } else {
+ return site_discarded ;
+ }
+
+ }
+
+ site = temp_site;
+ chi2increment = site->GetDeltaChi2() ;
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+ _ktest->_diagnostics.record_site(kalhit, site);
+#endif
+
+ return success ;
+
+ }
+
+ int MarlinDDKalTestTrack::addAndFit( edm4hep::TrackerHit& trkhit, double& chi2increment, double maxChi2Increment) {
+
+ if (! trkhit.isAvailable() ) {
+ streamlog_out( ERROR) << "MarlinDDKalTestTrack::addAndFit( edm4hep::TrackerHit& trkhit, double& chi2increment, double maxChi2Increment): trkhit == 0" << std::endl;
+ return bad_intputs ;
+ }
+
+ const DDVMeasLayer* ml = _ktest->findMeasLayer( trkhit ) ;
+
+ if( ml == 0 ){
+ // fg: not sure if ml should ever be 0 - but it seems to happen,
+ // if point is not on surface and more than one surface exists ...
+
+ streamlog_out( ERROR ) << ">>>>>>>>>>> no measurment layer found for trkhit cellid0 : "
+ << cellIDString( trkhit.getCellID() ) << " at "
+ << edm4hep::Vector3d( trkhit.getPosition() ) << std::endl ;
+
+ return IMarlinTrack::bad_intputs ;
+ }
+
+ DDVTrackHit* kalhit = ml->ConvertLCIOTrkHit(trkhit) ;
+
+ if( kalhit == 0 ){ //fg: ml->ConvertLCIOTrkHit returns 0 if hit not on surface !!!
+ return IMarlinTrack::bad_intputs ;
+ }
+
+ TKalTrackSite* site = 0 ;
+ int error_code = this->addAndFit( kalhit, chi2increment, site, maxChi2Increment);
+
+ if( error_code != success ){
+
+ delete kalhit;
+
+ // if the hit fails for any reason other than the Chi2 cut record the Chi2 contibution as DBL_MAX
+ if( error_code != site_fails_chi2_cut ) {
+ chi2increment = DBL_MAX;
+ }
+
+ _outlier_chi2_values.push_back(std::make_pair(trkhit, chi2increment));
+
+ streamlog_out( DEBUG2 ) << ">>>>>>>>>>> addAndFit Number of Outliers : "
+ << _outlier_chi2_values.size() << std::endl;
+
+ return error_code ;
+ }
+ else {
+ this->addHit( trkhit, kalhit, ml ) ;
+ _hit_used_for_sites[trkhit] = site ;
+ _hit_chi2_values.push_back(std::make_pair(trkhit, chi2increment));
+ }
+
+ // set the values for the point at which the fit becomes constained
+ if( !_trackHitAtPositiveNDF.isAvailable() && _kaltrack->GetNDF() >= 0){
+
+ _trackHitAtPositiveNDF = trkhit;
+ _hitIndexAtPositiveNDF = _kaltrack->IndexOf( site );
+
+ streamlog_out( DEBUG2 ) << ">>>>>>>>>>> Fit is now constrained at : "
+ << cellIDString( trkhit.getCellID() )
+ << " pos " << edm4hep::Vector3d( trkhit.getPosition() )
+ << " trkhit = " << _trackHitAtPositiveNDF
+ << " index of kalhit = " << _hitIndexAtPositiveNDF
+ << " NDF = " << _kaltrack->GetNDF()
+ << std::endl;
+
+ }
+
+ return success ;
+
+ }
+
+
+
+ int MarlinDDKalTestTrack::testChi2Increment( edm4hep::TrackerHit& trkhit, double& chi2increment ) {
+
+ if (! trkhit.isAvailable() ) {
+ streamlog_out( ERROR) << "MarlinDDKalTestTrack::addAndFit( edm4hep::TrackerHit& trkhit, double& chi2increment, double maxChi2Increment): trkhit == 0" << std::endl;
+ return IMarlinTrack::bad_intputs ;
+ }
+
+ const DDVMeasLayer* ml = _ktest->findMeasLayer( trkhit ) ;
+
+ if( ml == 0 ){
+ // fg: not sure if ml should ever be 0 - but it seems to happen,
+ // if point is not on surface and more than one surface exists ...
+
+ streamlog_out( ERROR ) << ">>>>>>>>>>> no measurment layer found for trkhit cellid0 : "
+ << cellIDString( trkhit.getCellID() ) << " at "
+ << edm4hep::Vector3d( trkhit.getPosition() ) << std::endl ;
+
+ return IMarlinTrack::bad_intputs ;
+
+ }
+
+ DDVTrackHit* kalhit = ml->ConvertLCIOTrkHit(trkhit) ;
+
+ if( kalhit == 0 ){ //fg: ml->ConvertLCIOTrkHit returns 0 if hit not on surface !!!
+ return IMarlinTrack::bad_intputs ;
+ }
+
+
+ TKalTrackSite* site = 0 ;
+ int error_code = this->addAndFit( kalhit, chi2increment, site, -DBL_MAX); // using -DBL_MAX here ensures the hit will never be added to the fit
+
+ delete kalhit;
+
+ return error_code;
+
+ }
+
+
+
+ int MarlinDDKalTestTrack::fit( double maxChi2Increment ) {
+
+ // SJA:FIXME: what do we do about calling fit after we have already added hits and filtered
+ // I guess this would created new sites when addAndFit is called
+ // one option would be to remove the sites
+ // need to check where the sites are stored ... probably in the KalTrackSystem
+ //
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::fit() called " << std::endl ;
+
+ if ( ! _initialised ) {
+
+ throw MarlinTrk::Exception("Track fit not initialised");
+
+ }
+
+ // ---------------------------
+ // Prepare hit iterrator for adding hits to kaltrack
+ // ---------------------------
+
+ TIter next(_kalhits, _fitDirection);
+
+ // ---------------------------
+ // Start Kalman Filter
+ // ---------------------------
+
+ DDVTrackHit *kalhit = 0;
+
+ while ( (kalhit = dynamic_cast<DDVTrackHit *>( next() ) ) ) {
+
+ double chi2increment;
+ TKalTrackSite* site=nullptr;
+ int error_code = this->addAndFit( kalhit, chi2increment, site, maxChi2Increment );
+
+
+ edm4hep::TrackerHit trkhit = kalhit->getLCIOTrackerHit();
+
+ if( error_code == 0 ){ // add trkhit to map associating trkhits and sites
+ _hit_used_for_sites[trkhit] = site;
+ _hit_chi2_values.push_back(std::make_pair(trkhit, chi2increment));
+
+ // set the values for the point at which the fit becomes constained
+ if( !_trackHitAtPositiveNDF.isAvailable() && _kaltrack->GetNDF() >= 0){
+
+ _trackHitAtPositiveNDF = trkhit;
+ _hitIndexAtPositiveNDF = _kaltrack->IndexOf( site );
+
+ streamlog_out( DEBUG2 ) << ">>>>>>>>>>> Fit is now constrained at : "
+ << cellIDString( trkhit.getCellID() )
+ << " pos " << edm4hep::Vector3d( trkhit.getPosition() )
+ << " trkhit = " << _trackHitAtPositiveNDF
+ << " index of kalhit = " << _hitIndexAtPositiveNDF
+ << " NDF = " << _kaltrack->GetNDF()
+ << std::endl;
+ }
+
+ }
+ else { // hit rejected by the filter, so store in the list of rejected hits
+
+ // if the hit fails for any reason other than the Chi2 cut record the Chi2 contibution as DBL_MAX
+ if( error_code != site_fails_chi2_cut ) {
+ chi2increment = DBL_MAX;
+ }
+
+ _outlier_chi2_values.push_back(std::make_pair(trkhit, chi2increment));
+ streamlog_out( DEBUG2 ) << ">>>>>>>>>>> fit(): Number of Outliers : " << _outlier_chi2_values.size() << std::endl;
+ std::cout << ">>>>>>>>>>> fit(): Number of Outliers : " << _outlier_chi2_values.size() << std::endl;
+ _hit_not_used_for_sites.push_back(trkhit) ;
+
+ }
+
+ } // end of Kalman filter
+
+ if( _ktest->getOption( MarlinTrk::IMarlinTrkSystem::CFG::useSmoothing ) ){
+ streamlog_out( DEBUG2 ) << "Perform Smoothing for All Previous Measurement Sites " << std::endl ;
+ int error = this->smooth() ;
+
+ if( error != success ) return error ;
+
+ }
+
+ //return _hit_used_for_sites.empty() == false ? success : all_sites_fail_fit ;
+ if( _hit_used_for_sites.empty() == false )
+ {
+ return success ;
+ }
+ else{
+ return all_sites_fail_fit ;
+ }
+
+ }
+
+
+ /** smooth all track states
+ */
+ int MarlinDDKalTestTrack::smooth(){
+
+ streamlog_out( DEBUG2 ) << "MarlinDDKalTestTrack::smooth() " << std::endl ;
+
+ //fg: we should actually smooth all sites - it is then up to the user which smoothed tracks state to take
+ // for any furthter extrapolation/propagation ...
+
+ if( !_smoothed )
+ _kaltrack->SmoothAll() ;
+
+ //SJA:FIXME: in the current implementation it is only possible to smooth back to the 4th site.
+ // This is due to the fact that the covariance matrix is not well defined at the first 3 measurement sites filtered.
+
+ // _kaltrack->SmoothBackTo( _hitIndexAtPositiveNDF + 1 ) ;
+
+ _smoothed = true ;
+
+ return success ;
+
+ }
+
+
+ /** smooth track states from the last filtered hit back to the measurement site associated with the given hit
+ */
+ int MarlinDDKalTestTrack::smooth( edm4hep::TrackerHit& trkhit ) {
+
+ streamlog_out( DEBUG2 ) << "MarlinDDKalTestTrack::smooth( edm4hep::TrackerHit& " << trkhit << " ) " << std::endl ;
+
+
+ if (! trkhit.isAvailable() ) {
+ return bad_intputs ;
+ }
+
+
+ std::map<edm4hep::TrackerHit,TKalTrackSite*>::const_iterator it;
+
+
+ TKalTrackSite* site = 0 ;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code ;
+
+ int index = _kaltrack->IndexOf( site );
+
+ _kaltrack->SmoothBackTo( index ) ;
+
+ _smoothed = true ;
+
+ return success ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::getTrackState( edm4hep::TrackState& ts, double& chi2, int& ndf ) {
+
+ streamlog_out( DEBUG2 ) << "MarlinDDKalTestTrack::getTrackState( edm4hep::TrackState& ts ) " << std::endl ;
+
+ // use the last filtered track state
+ const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
+
+ this->ToLCIOTrackState( site, ts, chi2, ndf );
+
+ return success ;
+
+
+ }
+
+
+ int MarlinDDKalTestTrack::getTrackState( edm4hep::TrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf ) {
+
+ streamlog_out( DEBUG2 ) << "MarlinDDKalTestTrack::getTrackState( edm4hep::TrackerHit& trkhit, edm4hep::TrackState& ts ) using hit: " << trkhit << " with cellID0 = " << trkhit.getCellID() << std::endl ;
+
+ TKalTrackSite* site = 0 ;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code ;
+
+ streamlog_out( DEBUG1 ) << "MarlinDDKalTestTrack::getTrackState: site " << site << std::endl;
+
+ this->ToLCIOTrackState( *site, ts, chi2, ndf );
+
+ return success ;
+ }
+
+
+ int MarlinDDKalTestTrack::getHitsInFit( std::vector<std::pair<edm4hep::TrackerHit, double> >& hits ) {
+
+ std::copy( _hit_chi2_values.begin() , _hit_chi2_values.end() , std::back_inserter( hits ) ) ;
+
+ // this needs more thought. What about when the hits are added using addAndFit?
+
+ // need to check the order so that we can return the list ordered in time
+ // as they will be added to _hit_chi2_values in the order of fitting
+ // not in the order of time
+//
+// if( _fitDirection == IMarlinTrack::backward ){
+// std::reverse_copy( _hit_chi2_values.begin() , _hit_chi2_values.end() , std::back_inserter( hits ) ) ;
+// } else {
+// std::copy( _hit_chi2_values.begin() , _hit_chi2_values.end() , std::back_inserter( hits ) ) ;
+// }
+
+ return success ;
+
+ }
+
+ int MarlinDDKalTestTrack::getOutliers( std::vector<std::pair<edm4hep::TrackerHit, double> >& hits ) {
+
+ std::copy( _outlier_chi2_values.begin() , _outlier_chi2_values.end() , std::back_inserter( hits ) ) ;
+
+
+ // this needs more thought. What about when the hits are added using addAndFit?
+// // need to check the order so that we can return the list ordered in time
+// // as they will be added to _hit_chi2_values in the order of fitting
+// // not in the order of time
+//
+// if( _fitDirection == IMarlinTrack::backward ){
+// std::reverse_copy( _outlier_chi2_values.begin() , _outlier_chi2_values.end() , std::back_inserter( hits ) ) ;
+// } else {
+// std::copy( _outlier_chi2_values.begin() , _outlier_chi2_values.end() , std::back_inserter( hits ) ) ;
+// }
+
+
+ return success ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::getNDF( int& ndf ){
+
+ if( _initialised == false ) {
+ return error;
+ } else {
+
+ ndf = _kaltrack->GetNDF();
+ return success;
+
+ }
+
+ }
+
+ int MarlinDDKalTestTrack::getTrackerHitAtPositiveNDF( edm4hep::TrackerHit& trkhit ) {
+
+ trkhit = _trackHitAtPositiveNDF;
+ return success;
+
+ }
+
+
+ int MarlinDDKalTestTrack::extrapolate( const edm4hep::Vector3d& point, edm4hep::TrackState& ts, double& chi2, int& ndf ){
+
+ const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
+
+ return this->extrapolate( point, site, ts, chi2, ndf ) ;
+
+ }
+
+ int MarlinDDKalTestTrack::extrapolate( const edm4hep::Vector3d& point, edm4hep::TrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf ) {
+
+ TKalTrackSite* site = 0 ;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code;
+
+ return this->extrapolate( point, *site, ts, chi2, ndf ) ;
+
+ }
+
+ int MarlinDDKalTestTrack::extrapolate( const edm4hep::Vector3d& point, const TKalTrackSite& site ,edm4hep::TrackState& ts, double& chi2, int& ndf ){
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::extrapolate( const edm4hep::Vector3d& point, edm4hep::TrackState& ts, double& chi2, int& ndf ) called " << std::endl ;
+
+ TKalTrackState& trkState = (TKalTrackState&) site.GetCurState(); // this segfaults if no hits are present
+
+ THelicalTrack helix = trkState.GetHelix() ;
+ double dPhi ;
+
+ const TVector3 tpoint( point.x, point.y, point.z ) ;
+
+ Int_t sdim = trkState.GetDimension(); // dimensions of the track state, it will be 5 or 6
+ TKalMatrix sv(sdim,1);
+
+ // now move to the point
+ TKalMatrix DF(sdim,sdim);
+ DF.UnitMatrix();
+ helix.MoveTo( tpoint , dPhi , &DF , 0) ; // move helix to desired point, and get propagator matrix
+
+ TMatrixD c0(trkState.GetCovMat());
+
+ TKalMatrix DFt = TKalMatrix(TMatrixD::kTransposed, DF);
+ c0 = DF * c0 * DFt ; // update the covariance matrix
+
+ this->ToLCIOTrackState( helix, c0, ts, chi2, ndf );
+
+ return success;
+
+ }
+
+
+ int MarlinDDKalTestTrack::extrapolateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+
+ const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
+
+ return this->extrapolateToLayer( layerID, site, ts, chi2, ndf, detElementID, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::extrapolateToLayer( int layerID, edm4hep::TrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+
+ TKalTrackSite* site = 0;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code ;
+
+ return this->extrapolateToLayer( layerID, *site, ts, chi2, ndf, detElementID, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::extrapolateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::extrapolateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID ) called " << std::endl ;
+
+ edm4hep::Vector3d crossing_point ;
+ const DDVMeasLayer* ml = 0;
+
+ int error_code = this->intersectionWithLayer( layerID, site, crossing_point, detElementID, ml, mode ) ;
+
+ if( error_code != 0 ) return error_code ;
+
+ return this->extrapolate( crossing_point, site, ts, chi2, ndf ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::extrapolateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+
+ const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
+
+ return this->extrapolateToDetElement( detElementID, site, ts, chi2, ndf, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::extrapolateToDetElement( int detElementID, edm4hep::TrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+
+ TKalTrackSite* site = 0;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code ;
+
+ return this->extrapolateToDetElement( detElementID, *site, ts, chi2, ndf, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::extrapolateToDetElement( int detElementID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::extrapolateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) called " << std::endl ;
+
+ edm4hep::Vector3d crossing_point ;
+
+ const DDVMeasLayer* ml = 0;
+ int error_code = this->intersectionWithDetElement( detElementID, site, crossing_point, ml, mode ) ;
+
+ if( error_code != 0 ) return error_code ;
+
+ return this->extrapolate( crossing_point, site, ts, chi2, ndf ) ;
+
+ }
+
+
+
+ int MarlinDDKalTestTrack::propagate( const edm4hep::Vector3d& point, edm4hep::TrackState& ts, double& chi2, int& ndf ){
+
+ const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
+
+ // check if the point is inside the beampipe
+ // SJA:FIXME: really we should also check if the PCA to the point is also less than R
+ const DDVMeasLayer* ml = ( (_ktest->getIPLayer() ) && sqrt(point.x*point.x+point.y*point.y+point.z*point.z) < _ktest->getIPLayer()->GetR()) ? _ktest->getIPLayer() : 0;
+
+ return this->propagate( point, site, ts, chi2, ndf, ml ) ;
+
+ }
+
+ int MarlinDDKalTestTrack::propagate( const edm4hep::Vector3d& point, edm4hep::TrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf ){
+
+ TKalTrackSite* site = 0;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code ;
+
+ // check if the point is inside the beampipe
+ // SJA:FIXME: really we should also check if the PCA to the point is also less than R
+
+ const DDVMeasLayer* ml = _ktest->getIPLayer();
+
+ if ( ml )
+ if (sqrt(point.x*point.x+point.y*point.y+point.z*point.z) > _ktest->getIPLayer()->GetR()) ml = NULL;
+
+// const DDVMeasLayer* ml = (point.r() < _ktest->getIPLayer()->GetR()) ? _ktest->getIPLayer() : 0;
+
+ return this->propagate( point, *site, ts, chi2, ndf, ml ) ;
+
+ }
+
+ int MarlinDDKalTestTrack::propagate( const edm4hep::Vector3d& point, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, const DDVMeasLayer* ml ){
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::propagate( const edm4hep::Vector3d& point, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf ) called " << std::endl ;
+
+ const TVector3 tpoint( point.x, point.y, point.z ) ;
+
+ TKalTrackState& trkState = (TKalTrackState&) site.GetCurState(); // this segfaults if no hits are present
+
+ THelicalTrack helix = trkState.GetHelix() ;
+ double dPhi = 0.0;
+
+
+ Int_t sdim = trkState.GetDimension(); // dimensions of the track state, it will be 5 or 6
+ TKalMatrix sv(sdim,1);
+
+ TKalMatrix F(sdim,sdim); // propagator matrix to be returned by transport function
+ F.UnitMatrix(); // set the propagator matrix to the unit matrix
+
+ TKalMatrix Q(sdim,sdim); // noise matrix to be returned by transport function
+ Q.Zero();
+ TVector3 x0; // intersection point to be returned by transport
+
+ TMatrixD c0(trkState.GetCovMat());
+
+ // the last layer crossed by the track before point
+ if( ! ml ){
+ ml = _ktest->getLastMeasLayer(helix, tpoint);
+ }
+
+ if ( ml ) {
+
+ _ktest->_det->Transport(site, *ml, x0, sv, F, Q ) ; // transport to last layer cross before point
+
+ // given that we are sure to have intersected the layer ml as this was provided via getLastMeasLayer, x0 will lie on the layer
+ // this could be checked with the method isOnSurface
+ // so F will be the propagation matrix from the current location to the last surface and Q will be the noise matrix up to this point
+
+
+ TKalMatrix Ft = TKalMatrix(TMatrixD::kTransposed, F);
+ c0 = F * c0 * Ft + Q; // update covaraince matrix and add the MS assosiated with moving to tvml
+
+ helix.MoveTo( x0 , dPhi , 0 , 0 ) ; // move the helix to tvml
+
+
+ }
+ else { // the current site is at the last surface before the point to propagate to
+
+ ml = dynamic_cast<const DDVMeasLayer*>(&(site.GetHit().GetMeasLayer())) ;
+
+ }
+
+ // get whether the track is incomming or outgoing at the last surface
+ const TVSurface *sfp = dynamic_cast<const TVSurface *>(ml); // last surface
+
+ TMatrixD dxdphi = helix.CalcDxDphi(0); // tangent vector at last surface
+ TVector3 dxdphiv(dxdphi(0,0),dxdphi(1,0),dxdphi(2,0)); // convert matirix diagonal to vector
+// Double_t cpa = helix.GetKappa(); // get pt
+
+ Bool_t isout = -dPhi*dxdphiv.Dot(sfp->GetOutwardNormal(x0)) < 0 ? kTRUE : kFALSE; // out-going or in-coming at the destination surface
+
+ // now move to the point
+ TKalMatrix DF(sdim,sdim);
+ DF.UnitMatrix();
+ helix.MoveTo( tpoint , dPhi , &DF , 0) ; // move helix to desired point, and get propagator matrix
+
+ TKalMatrix Qms(sdim, sdim);
+ ml->CalcQms(isout, helix, dPhi, Qms); // calculate MS for the final step through the present material
+
+ TKalMatrix DFt = TKalMatrix(TMatrixD::kTransposed, DF);
+ c0 = DF * c0 * DFt + Qms ; // update the covariance matrix
+
+
+ this->ToLCIOTrackState( helix, c0, ts, chi2, ndf );
+
+ return success;
+
+ }
+
+
+ int MarlinDDKalTestTrack::propagateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+
+ const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
+
+ return this->propagateToLayer( layerID, site, ts, chi2, ndf, detElementID, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::propagateToLayer( int layerID, edm4hep::TrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+
+ TKalTrackSite* site = 0;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code ;
+
+ return this->propagateToLayer( layerID, *site, ts, chi2, ndf, detElementID, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::propagateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::propagateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID ) called " << std::endl;
+
+ edm4hep::Vector3d crossing_point ;
+
+ const DDVMeasLayer* ml = 0;
+
+ int error_code = this->intersectionWithLayer( layerID, site, crossing_point, detElementID, ml, mode) ;
+
+ if( error_code != success ) return error_code ;
+
+ return this->propagate( crossing_point, site, ts, chi2, ndf , ml) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::propagateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+
+ const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
+
+ return this->propagateToDetElement( detElementID, site, ts, chi2, ndf, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::propagateToDetElement( int detElementID, edm4hep::TrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+
+ TKalTrackSite* site = 0;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code ;
+
+ return this->propagateToDetElement( detElementID, *site, ts, chi2, ndf, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::propagateToDetElement( int detElementID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::propagateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) called " << std::endl ;
+
+ edm4hep::Vector3d crossing_point ;
+
+ const DDVMeasLayer* ml = 0;
+ int error_code = this->intersectionWithDetElement( detElementID, site, crossing_point, ml, mode ) ;
+
+ if( error_code != 0 ) return error_code ;
+
+ return this->propagate( crossing_point, site, ts, chi2, ndf, ml ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::intersectionWithDetElement( int detElementID, edm4hep::Vector3d& point, int mode ) {
+
+ const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
+
+ const DDVMeasLayer* ml = 0;
+ return this->intersectionWithDetElement( detElementID, site, point, ml, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::intersectionWithDetElement( int detElementID, edm4hep::TrackerHit& trkhit, edm4hep::Vector3d& point, int mode ) {
+
+ TKalTrackSite* site = 0;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code ;
+
+ const DDVMeasLayer* ml = 0;
+ return this->intersectionWithDetElement( detElementID, *site, point, ml, mode ) ;
+
+ }
+
+ int MarlinDDKalTestTrack::intersectionWithDetElement( int detElementID, const TKalTrackSite& site, edm4hep::Vector3d& point, const DDVMeasLayer*& ml, int mode ) {
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::intersectionWithDetElement( int detElementID, const TKalTrackSite& site, edm4hep::Vector3d& point, const DDVMeasLayer*& ml, int mode) called " << std::endl;
+
+ std::vector<const DDVMeasLayer*> meas_modules ;
+ _ktest->getSensitiveMeasurementModules( detElementID, meas_modules ) ;
+
+ if( meas_modules.size() == 0 ) {
+
+ std::stringstream errorMsg;
+ errorMsg << "MarlinDDKalTestTrack::intersectionWithDetElement detector element id unkown: detElementID = "
+ << cellIDString( detElementID ) << std::endl ;
+
+ throw MarlinTrk::Exception(errorMsg.str());
+
+ }
+
+ int dummy_detElementID; // not returned here as this is a single element as far as the outside world is concerned. Could check they are equal if we wanted ...
+ int error_code = this->findIntersection( meas_modules, site, point, dummy_detElementID, ml, mode ) ;
+
+ if ( error_code == success ) {
+ streamlog_out(DEBUG1) << "MarlinDDKalTestTrack::intersectionWithDetElement intersection with detElementID = "
+ << cellIDString( detElementID )
+ << ": at " << point
+ << std::endl ;
+ }
+
+ else if ( error_code == no_intersection ) {
+ ml = 0;
+
+ streamlog_out(DEBUG1) << "MarlinDDKalTestTrack::intersectionWithDetElement No intersection with detElementID = "
+ << cellIDString( detElementID )
+ << std::endl ;
+ }
+
+ return error_code ;
+
+ }
+
+ int MarlinDDKalTestTrack::intersectionWithLayer( int layerID, edm4hep::Vector3d& point, int& detElementID, int mode ) {
+
+ const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
+ const DDVMeasLayer* ml = 0;
+ return this->intersectionWithLayer( layerID, site, point, detElementID, ml, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::intersectionWithLayer( int layerID, edm4hep::TrackerHit& trkhit, edm4hep::Vector3d& point, int& detElementID, int mode ) {
+
+ TKalTrackSite* site = 0;
+ int error_code = getSiteFromLCIOHit(trkhit, site);
+
+ if( error_code != success ) return error_code ;
+
+ const DDVMeasLayer* ml = 0;
+ return this->intersectionWithLayer( layerID, *site, point, detElementID, ml, mode ) ;
+
+ }
+
+
+ int MarlinDDKalTestTrack::intersectionWithLayer( int layerID, const TKalTrackSite& site, edm4hep::Vector3d& point, int& detElementID, const DDVMeasLayer*& ml, int mode ) {
+
+ streamlog_out(DEBUG2) << "MarlinDDKalTestTrack::intersectionWithLayer( int layerID, const TKalTrackSite& site, edm4hep::Vector3d& point, int& detElementID, int mode) called layerID = " << layerID << std::endl;
+
+ std::vector<DDVMeasLayer const*> meas_modules ;
+ _ktest->getSensitiveMeasurementModulesForLayer( layerID, meas_modules ) ;
+
+ if( meas_modules.size() == 0 ) {
+ std::cout << "MarlinDDKalTestTrack::intersectionWithLayer layer id unknown: layerID = " << layerID << std::endl;
+ streamlog_out(DEBUG5)<< "MarlinDDKalTestTrack::intersectionWithLayer layer id unknown: layerID = " << cellIDString( layerID ) << std::endl ;
+ return no_intersection;
+
+ }
+
+ // int index_of_intersected;
+ int error_code = this->findIntersection( meas_modules, site, point, detElementID, ml, mode ) ;
+
+ if ( error_code == success ) {
+ streamlog_out(DEBUG1) << "MarlinDDKalTestTrack::intersectionWithLayer intersection with layerID = "
+ << layerID
+ << ": at " << point
+ << " detElementID = " << detElementID
+ << " " << cellIDString( detElementID )
+ << std::endl;
+ }
+ else if ( error_code == no_intersection ) {
+
+ ml = 0;
+ streamlog_out(DEBUG1) << "MarlinDDKalTestTrack::intersectionWithLayer No intersection with layerID = "
+ << layerID
+ << " " << cellIDString( layerID )
+ << std::endl;
+ }
+
+ return error_code;
+ }
+
+
+ int MarlinDDKalTestTrack::findIntersection( const DDVMeasLayer& meas_module, const TKalTrackSite& site, edm4hep::Vector3d& point, double& dphi, int& detElementID, int mode ) {
+
+
+ TKalTrackState& trkState = (TKalTrackState&) site.GetCurState();
+
+
+ //--------- DEBUG --------------
+ // TKalTrackState* tsSmoothed = ( &((TVKalSite&)site).GetState(TVKalSite::kSmoothed) != 0 ?
+ // &(TKalTrackState&) ((TVKalSite&)site).GetState( TVKalSite::kSmoothed ) : 0 ) ;
+ // if( tsSmoothed == &trkState )
+ // streamlog_out(DEBUG2) << "************ MarlinDDKalTestTrack::intersectionWithLayer : using smoothed TrackState !!!!! " << std::endl ;
+
+ // TKalTrackState* tsFiltered = ( &((TVKalSite&)site).GetState(TVKalSite::kFiltered) != 0 ?
+ // &(TKalTrackState&) ((TVKalSite&)site).GetState( TVKalSite::kFiltered ) : 0 ) ;
+ // if( tsFiltered == &trkState )
+ // streamlog_out(DEBUG2) << "************ MarlinDDKalTestTrack::intersectionWithLayer : using filtered TrackState !!!!! " << std::endl ;
+ // //------------------------------
+
+
+ THelicalTrack helix = trkState.GetHelix() ;
+
+ TVector3 xto; // reference point at destination to be returned by CalcXinPointWith
+
+ int crossing_exist = meas_module.getIntersectionAndCellID(helix, xto, dphi, detElementID, mode);
+ // int crossing_exist = surf->CalcXingPointWith(helix, xto, dphi, mode) ;
+ std::cout << "MarlinDDKalTestTrack::intersectionWithLayer crossing_exist = " << crossing_exist << " dphi " << dphi << " with detElementIDs: " << detElementID << std::endl;
+ streamlog_out(DEBUG1) << "MarlinDDKalTestTrack::intersectionWithLayer crossing_exist = " << crossing_exist << " dphi " << dphi << " with detElementIDs: " << detElementID ;
+
+ streamlog_out(DEBUG1) << std::endl ;
+
+ if( crossing_exist == 0 ) {
+ return no_intersection ;
+ }
+ else {
+ point.x = xto.X();
+ point.y = xto.Y();
+ point.z = xto.Z();
+ }
+
+ return success;
+ }
+
+
+ int MarlinDDKalTestTrack::findIntersection( std::vector<DDVMeasLayer const*>& meas_modules, const TKalTrackSite& site, edm4hep::Vector3d& point, int& detElementID, const DDVMeasLayer*& ml, int mode ) {
+
+ unsigned int n_modules = meas_modules.size() ;
+
+ double dphi_min = DBL_MAX; // use to store the min deflection angle found so that can avoid the crossing on the far side of the layer
+ bool surf_found(false);
+
+ for( unsigned int i = 0 ; i < n_modules ; ++i ){
+
+ double dphi = 0;
+ // need to send a temporary point as we may get the crossing point with the layer on the oposite side of the layer
+ edm4hep::Vector3d point_temp ;
+
+ int temp_detElementID;
+
+ int error_code = findIntersection( *meas_modules[i], site, point_temp, dphi, temp_detElementID, mode ) ;
+
+ if( error_code == success ) {
+
+ // make sure we get the next crossing
+ if( fabs(dphi) < dphi_min ) {
+
+ dphi_min = fabs(dphi) ;
+ surf_found = true ;
+ ml = meas_modules[i];
+ detElementID = temp_detElementID;
+ point = point_temp ;
+ }
+
+ }
+ else if( error_code != no_intersection ) { // in which case error_code is an error rather than simply a lack of intersection, so return
+
+ return error_code ;
+
+ }
+
+ }
+
+ // check if the surface was found and return accordingly
+ if ( surf_found ) {
+ return success ;
+ }
+ else {
+ return no_intersection ;
+ }
+
+
+ }
+
+
+
+ void MarlinDDKalTestTrack::ToLCIOTrackState( const THelicalTrack& helix, const TMatrixD& cov, edm4hep::TrackState& ts, double& chi2, int& ndf) const {
+
+ chi2 = _kaltrack->GetChi2();
+ ndf = _kaltrack->GetNDF();
+
+ //============== convert parameters to LCIO convention ====
+
+ // fill 5x5 covariance matrix from the 6x6 covariance matrix above
+ TMatrixD covK(5,5) ; for(int i=0;i<5;++i) for(int j=0;j<5;++j) covK[i][j] = cov[i][j] ;
+
+ // this is for incomming tracks ...
+ double phi = toBaseRange( helix.GetPhi0() + M_PI/2. ) ;
+ double omega = 1. /helix.GetRho() ;
+ double d0 = - helix.GetDrho() ;
+ double z0 = helix.GetDz() ;
+ double tanLambda = helix.GetTanLambda() ;
+
+ ts.D0 = d0;
+ ts.phi = phi; // fi0 - M_PI/2. ) ;
+ ts.omega = omega;
+ ts.Z0 = z0;
+ ts.tanLambda = tanLambda;
+
+ Double_t cpa = helix.GetKappa();
+ double alpha = omega / cpa ; // conversion factor for omega (1/R) to kappa (1/Pt)
+
+ decltype(ts.covMatrix) covLCIO{};
+ covLCIO[ 0] = covK( 0 , 0 ) ; // d0, d0
+
+ covLCIO[ 1] = - covK( 1 , 0 ) ; // phi0, d0
+ covLCIO[ 2] = covK( 1 , 1 ) ; // phi0, phi
+
+ covLCIO[ 3] = - covK( 2 , 0 ) * alpha ; // omega, d0
+ covLCIO[ 4] = covK( 2 , 1 ) * alpha ; // omega, phi
+ covLCIO[ 5] = covK( 2 , 2 ) * alpha * alpha ; // omega, omega
+
+ covLCIO[ 6] = - covK( 3 , 0 ) ; // z0 , d0
+ covLCIO[ 7] = covK( 3 , 1 ) ; // z0 , phi
+ covLCIO[ 8] = covK( 3 , 2 ) * alpha ; // z0 , omega
+ covLCIO[ 9] = covK( 3 , 3 ) ; // z0 , z0
+
+ covLCIO[10] = - covK( 4 , 0 ) ; // tanl, d0
+ covLCIO[11] = covK( 4 , 1 ) ; // tanl, phi
+ covLCIO[12] = covK( 4 , 2 ) * alpha ; // tanl, omega
+ covLCIO[13] = covK( 4 , 3 ) ; // tanl, z0
+ covLCIO[14] = covK( 4 , 4 ) ; // tanl, tanl
+
+
+ ts.covMatrix = covLCIO;
+
+
+ float pivot[3] ;
+
+ pivot[0] = helix.GetPivot().X() ;
+ pivot[1] = helix.GetPivot().Y() ;
+ pivot[2] = helix.GetPivot().Z() ;
+
+ ts.referencePoint = edm4hep::Vector3f(pivot);
+
+ streamlog_out( DEBUG2 ) << " kaltest track parameters: "
+ << " chi2/ndf " << chi2 / ndf
+ << " chi2 " << chi2
+ << " ndf " << ndf
+ << " prob " << TMath::Prob(chi2, ndf)
+ << std::endl
+
+ << "\t D0 " << d0 << "[+/-" << sqrt( covLCIO[0] ) << "] "
+ << "\t Phi :" << phi << "[+/-" << sqrt( covLCIO[2] ) << "] "
+ << "\t Omega " << omega << "[+/-" << sqrt( covLCIO[5] ) << "] "
+ << "\t Z0 " << z0 << "[+/-" << sqrt( covLCIO[9] ) << "] "
+ << "\t tan(Lambda) " << tanLambda << "[+/-" << sqrt( covLCIO[14]) << "] "
+
+ << "\t pivot : [" << pivot[0] << ", " << pivot[1] << ", " << pivot[2]
+ << " - r: " << std::sqrt( pivot[0]*pivot[0]+pivot[1]*pivot[1] ) << "]"
+ << std::endl ;
+
+
+ }
+
+
+ void MarlinDDKalTestTrack::ToLCIOTrackState( const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf ) const {
+
+ TKalTrackState& trkState = (TKalTrackState&) site.GetCurState(); // GetCutState will return the last added state to this site
+ // Assuming everything has proceeded as expected
+ // this will be Predicted -> Filtered -> Smoothed
+
+
+ // streamlog_out( DEBUG3 ) << " MarlinDDKalTestTrack::ToLCIOTrackState : " << std::endl ;
+ // trkState.DebugPrint() ;
+
+ THelicalTrack helix = trkState.GetHelix() ;
+
+ TMatrixD c0(trkState.GetCovMat());
+
+ this->ToLCIOTrackState( helix, c0, ts, chi2, ndf );
+
+ }
+
+
+ std::string MarlinDDKalTestTrack::toString() {
+
+ std::stringstream str ;
+
+ str << IMarlinTrack::toString() ;
+
+ str << _kaltrack->toString() ;
+
+ str << " --------------------- " << std::endl ;
+
+ return str.str() ;
+ }
+
+ int MarlinDDKalTestTrack::getSiteFromLCIOHit( edm4hep::TrackerHit& trkhit, TKalTrackSite*& site ) const {
+
+ std::map<edm4hep::TrackerHit,TKalTrackSite*>::const_iterator it;
+
+ it = _hit_used_for_sites.find(trkhit) ;
+
+ if( it == _hit_used_for_sites.end() ) { // hit not associated with any site
+
+ bool found = false;
+
+ for( unsigned int i = 0; i < _hit_not_used_for_sites.size(); ++i) {
+ if( trkhit == _hit_not_used_for_sites[i] ) found = true ;
+ }
+
+ if( found ) {
+ streamlog_out( DEBUG2 ) << "MarlinDDKalTestTrack::getSiteFromLCIOHit: hit was rejected during filtering" << std::endl ;
+ return site_discarded ;
+ }
+ else {
+ streamlog_out( DEBUG2 ) << "MarlinDDKalTestTrack::getSiteFromLCIOHit: hit " << trkhit << " not in list of supplied hits" << std::endl ;
+ return bad_intputs ;
+ }
+ }
+
+ site = it->second;
+
+
+ streamlog_out( DEBUG1 ) << "MarlinDDKalTestTrack::getSiteFromLCIOHit: site " << site << " found for hit " << trkhit << std::endl ;
+ return success ;
+
+ }
+
+} // end of namespace MarlinTrk
diff --git a/Service/TrackSystemSvc/src/MarlinDDKalTestTrack.h b/Service/TrackSystemSvc/src/MarlinDDKalTestTrack.h
new file mode 100644
index 0000000000000000000000000000000000000000..ccdf009f33a1c138c53db8e304deff92b483f2a5
--- /dev/null
+++ b/Service/TrackSystemSvc/src/MarlinDDKalTestTrack.h
@@ -0,0 +1,389 @@
+#ifndef MarlinDDKalTestTrack_h
+#define MarlinDDKalTestTrack_h
+
+#include "TrackSystemSvc/IMarlinTrack.h"
+#include "TrackSystemSvc/IMarlinTrkSystem.h"
+
+#include <TObjArray.h>
+
+#include <cmath>
+
+#include "TMatrixD.h"
+
+
+class TKalTrack ;
+class THelicalTrack ;
+class TKalTrackSite ;
+class DDVTrackHit ;
+class DDVMeasLayer ;
+
+namespace MarlinTrk {
+ class MarlinDDKalTest;
+}
+
+namespace edm4hep{
+ class TrackerHit ;
+}
+
+
+
+/** Implementation of the IMarlinTrack interface, using KalTest and KalDet to provide
+ * the needed functionality for a Kalman Filter.
+ *
+ * @version $Id: MarlinDDKalTestTrack.h 3641 2012-06-13 13:04:36Z aplin $
+ * @author S.Aplin, F. Gaede DESY
+ */
+
+namespace MarlinTrk{
+
+class MarlinDDKalTestTrack : public MarlinTrk::IMarlinTrack {
+
+public:
+
+#ifdef MARLINTRK_DIAGNOSTICS_ON
+ friend class DiagnosticsController;
+#endif
+
+ MarlinDDKalTestTrack(MarlinDDKalTest* ktest) ;
+
+ ~MarlinDDKalTestTrack() ;
+
+protected:
+
+private:
+
+ MarlinDDKalTestTrack(const MarlinDDKalTestTrack&) ; // Prevent copy-construction
+ MarlinDDKalTestTrack& operator=(const MarlinDDKalTestTrack&) ; // Prevent assignment
+
+ // make member functions private to force use through interface
+
+ /** set the mass of the charged particle (GeV) that is used for energy loss and multiple scattering -
+ * default value if this method is not called is the pion mass.
+ */
+ void setMass(double mass) ;
+
+ /** return the of the charged particle (GeV) that is used for energy loss and multiple scattering.
+ */
+ double getMass() ;
+
+ /** add hit to track - the hits have to be added ordered in time ( i.e. typically outgoing )
+ * this order will define the direction of the energy loss used in the fit
+ */
+ int addHit(edm4hep::TrackerHit& hit) ;
+
+ /** add hit to track - the hits have to be added ordered in time ( i.e. typically outgoing )
+ * this order will define the direction of the energy loss used in the fit
+ */
+ int addHit(edm4hep::TrackerHit& trkhit, const DDVMeasLayer* ml) ;
+
+ /** add hit to track - the hits have to be added ordered in time ( i.e. typically outgoing )
+ * this order will define the direction of the energy loss used in the fit
+ */
+ int addHit(edm4hep::TrackerHit& trkhit, DDVTrackHit* kalhit, const DDVMeasLayer* ml) ;
+
+ /** initialise the fit using the hits added up to this point -
+ * the fit direction has to be specified using IMarlinTrack::backward or IMarlinTrack::forward.
+ * this is the order wrt the order used in addHit() that will be used in the fit()
+ */
+ int initialise( bool fitDirection );
+
+ /** initialise the fit with a track state
+ * the fit direction has to be specified using IMarlinTrack::backward or IMarlinTrack::forward.
+ * this is the order that will be used in the fit().
+ * it is the users responsibility that the track state is consistent with the order
+ * of the hits used in addHit() ( i.e. the direction of energy loss )
+ * Note: the bfield_z is not taken from the argument but from the first hit
+ * should consider changing the interface ...
+ */
+ int initialise( const edm4hep::TrackState& ts, double /*bfield_z*/, bool fitDirection ) ;
+
+
+ /** perform the fit of all current hits, returns error code ( IMarlinTrack::success if no error ) .
+ * the fit will be performed in the order specified at initialise() wrt the order used in addHit(), i.e.
+ * IMarlinTrack::backward implies fitting from the outside to the inside for tracks comming from the IP.
+ */
+ int fit( double maxChi2Increment=DBL_MAX ) ;
+
+
+ /** smooth all track states
+ */
+ int smooth() ;
+
+
+ /** smooth track states from the last filtered hit back to the measurement site associated with the given hit
+ */
+ int smooth( edm4hep::TrackerHit& hit ) ;
+
+
+ /** update the current fit using the supplied hit, return code via int. Provides the Chi2 increment to the fit from adding the hit via reference.
+ * the given hit will not be added if chi2increment > maxChi2Increment.
+ */
+ int addAndFit( edm4hep::TrackerHit& hit, double& chi2increment, double maxChi2Increment=DBL_MAX ) ;
+
+ /** update the current fit using the supplied hit, return code via int. Provides the Chi2 increment to the fit from adding the hit via reference.
+ * the given hit will not be added if chi2increment > maxChi2Increment.
+ */
+ int addAndFit( DDVTrackHit* kalhit, double& chi2increment, TKalTrackSite*& site, double maxChi2Increment=DBL_MAX ) ;
+
+
+ /** obtain the chi2 increment which would result in adding the hit to the fit. This method will not alter the current fit, and the hit will not be stored in the list of hits or outliers
+ */
+ int testChi2Increment( edm4hep::TrackerHit& hit, double& chi2increment ) ;
+
+
+ // Track State Accessesors
+
+ /** get track state, returning TrackState, chi2 and ndf via reference
+ */
+ int getTrackState( edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+
+
+ /** get track state at measurement associated with the given hit, returning TrackState, chi2 and ndf via reference
+ */
+ int getTrackState( edm4hep::TrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+
+
+ /** get the list of hits included in the fit, together with the chi2 contributions of the hits.
+ * Pointers to the hits together with their chi2 contribution will be filled into a vector of
+ * pairs consitining of the pointer as the first part of the pair and the chi2 contribution as
+ * the second.
+ */
+ int getHitsInFit( std::vector<std::pair<edm4hep::TrackerHit, double> >& hits ) ;
+
+ /** get the list of hits which have been rejected by from the fit due to the a chi2 increment greater than threshold,
+ * Pointers to the hits together with their chi2 contribution will be filled into a vector of
+ * pairs consitining of the pointer as the first part of the pair and the chi2 contribution as
+ * the second.
+ */
+ int getOutliers( std::vector<std::pair<edm4hep::TrackerHit, double> >& hits ) ;
+
+
+ /** get the current number of degrees of freedom for the fit.
+ */
+ int getNDF( int& ndf ) ;
+
+ /** get TrackeHit at which fit became constrained, i.e. ndf >= 0
+ */
+ int getTrackerHitAtPositiveNDF( edm4hep::TrackerHit& trkhit ) ;
+
+ // PROPAGATORS
+
+ /** propagate the fit to the point of closest approach to the given point, returning TrackState, chi2 and ndf via reference
+ */
+ int propagate( const edm4hep::Vector3d& point, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+
+
+ /** propagate the fit at the measurement site associated with the given hit, to the point of closest approach to the given point,
+ * returning TrackState, chi2 and ndf via reference
+ */
+ int propagate( const edm4hep::Vector3d& point, edm4hep::TrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+
+
+ /** propagate the fit at the provided measurement site, to the point of closest approach to the given point,
+ * returning TrackState, chi2 and ndf via reference
+ */
+ int propagate( const edm4hep::Vector3d& point, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, const DDVMeasLayer* ml = 0 ) ;
+
+
+ /** propagate the fit to the numbered sensitive layer, returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
+ */
+ int propagateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
+
+ /** propagate the fit at the measurement site associated with the given hit, to numbered sensitive layer,
+ * returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
+ */
+ int propagateToLayer( int layerID, edm4hep::TrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
+
+ /** propagate the fit at the measurement site, to numbered sensitive layer,
+ * returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
+ */
+ int propagateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
+
+ /** propagate the fit to sensitive detector element, returning TrackState, chi2 and ndf via reference
+ */
+ int propagateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
+
+ /** propagate the fit at the measurement site associated with the given hit, to sensitive detector element,
+ * returning TrackState, chi2 and ndf via reference
+ */
+ int propagateToDetElement( int detEementID, edm4hep::TrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
+
+ /** propagate the fit at the measurement site, to sensitive detector element,
+ * returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
+ */
+ int propagateToDetElement( int detEementID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
+
+
+
+ // EXTRAPOLATORS
+
+ /** extrapolate the fit to the point of closest approach to the given point, returning TrackState, chi2 and ndf via reference
+ */
+ int extrapolate( const edm4hep::Vector3d& point, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+
+ /** extrapolate the fit at the measurement site associated with the given hit, to the point of closest approach to the given point,
+ * returning TrackState, chi2 and ndf via reference
+ */
+ int extrapolate( const edm4hep::Vector3d& point, edm4hep::TrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+
+ /** extrapolate the fit at the measurement site, to the point of closest approach to the given point,
+ * returning TrackState, chi2 and ndf via reference
+ */
+ int extrapolate( const edm4hep::Vector3d& point, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+
+ /** extrapolate the fit to numbered sensitive layer, returning TrackState via provided reference
+ */
+ int extrapolateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site associated with the given hit, to numbered sensitive layer,
+ * returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
+ */
+ int extrapolateToLayer( int layerID, edm4hep::TrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site, to numbered sensitive layer,
+ * returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
+ */
+ int extrapolateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
+
+ /** extrapolate the fit to sensitive detector element, returning TrackState, chi2 and ndf via reference
+ */
+ int extrapolateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site associated with the given hit, to sensitive detector element,
+ * returning TrackState, chi2 and ndf via reference
+ */
+ int extrapolateToDetElement( int detEementID, edm4hep::TrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site, to sensitive detector element,
+ * returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
+ */
+ int extrapolateToDetElement( int detEementID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
+
+
+ /** Dump this track to a string for debugging.
+ */
+ std::string toString() ;
+
+ // INTERSECTORS
+
+
+ /** extrapolate the fit to numbered sensitive layer, returning intersection point in global coordinates and integer ID of the
+ * intersected sensitive detector element via reference
+ */
+ int intersectionWithLayer( int layerID, edm4hep::Vector3d& point, int& detElementID, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site associated with the given hit, to numbered sensitive layer,
+ * returning intersection point in global coordinates and integer ID of the intersected sensitive detector element via reference
+ */
+ int intersectionWithLayer( int layerID, edm4hep::TrackerHit& hit, edm4hep::Vector3d& point, int& detElementID, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site, to numbered sensitive layer,
+ * returning intersection point in global coordinates and integer ID of the intersected sensitive detector element via reference
+ */
+ int intersectionWithLayer( int layerID, const TKalTrackSite& site, edm4hep::Vector3d& point, int& detElementID, const DDVMeasLayer*& ml, int mode=modeClosest ) ;
+
+
+ /** extrapolate the fit to numbered sensitive detector element, returning intersection point in global coordinates via reference
+ */
+ int intersectionWithDetElement( int detElementID, edm4hep::Vector3d& point, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site associated with the given hit, to sensitive detector element,
+ * returning intersection point in global coordinates via reference
+ */
+ int intersectionWithDetElement( int detElementID, edm4hep::TrackerHit& hit, edm4hep::Vector3d& point, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site, to sensitive detector element,
+ * returning intersection point in global coordinates via reference
+ */
+ int intersectionWithDetElement( int detElementID, const TKalTrackSite& site, edm4hep::Vector3d& point, const DDVMeasLayer*& ml, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site, to sensitive detector elements contained in the std::vector,
+ * and return intersection point in global coordinates via reference
+ */
+ int findIntersection( std::vector<DDVMeasLayer const*>& meas_modules, const TKalTrackSite& site, edm4hep::Vector3d& point, int& detElementID, const DDVMeasLayer*& ml, int mode=modeClosest ) ;
+
+ /** extrapolate the fit at the measurement site, to the DDVMeasLayer,
+ * and return intersection point in global coordinates via reference
+ */
+ int findIntersection( const DDVMeasLayer& meas_module, const TKalTrackSite& site, edm4hep::Vector3d& point, double& dphi, int& detElementIDconst, int mode=modeClosest ) ;
+
+
+
+
+ //** end of memeber functions from IMarlinTrack interface
+
+ /** fill LCIO Track State with parameters from helix and cov matrix
+ */
+ void ToLCIOTrackState( const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf ) const ;
+
+ /** fill LCIO Track State with parameters from helix and cov matrix
+ */
+ void ToLCIOTrackState( const THelicalTrack& helix, const TMatrixD& cov, edm4hep::TrackState& ts, double& chi2, int& ndf ) const ;
+
+ /** get the measurement site associated with the given lcio TrackerHit trkhit
+ */
+ int getSiteFromLCIOHit( edm4hep::TrackerHit& trkhit, TKalTrackSite*& site ) const ;
+
+
+
+ /** helper function to restrict the range of the azimuthal angle to ]-pi,pi]*/
+ inline double toBaseRange( double phi) const {
+ while( phi <= -M_PI ){ phi += 2. * M_PI ; }
+ while( phi > M_PI ){ phi -= 2. * M_PI ; }
+ return phi ;
+ }
+
+
+
+ // memeber variables
+
+ TKalTrack* _kaltrack=nullptr;
+
+ std::vector<edm4hep::TrackerHit> _lcioHits{};
+
+ TObjArray* _kalhits=nullptr;
+
+ MarlinDDKalTest* _ktest=nullptr;
+
+ edm4hep::TrackerHit _trackHitAtPositiveNDF=nullptr;
+ int _hitIndexAtPositiveNDF=-1;
+
+ /** used to store whether initial track state has been supplied or created
+ */
+ bool _initialised=false;
+
+ /** used to store the fit direction supplied to intialise
+ */
+ bool _fitDirection=false;
+
+
+ /** used to store whether smoothing has been performed
+ */
+ bool _smoothed=false;
+
+ /** map to store relation between lcio hits and measurement sites
+ */
+ std::map<edm4hep::TrackerHit,TKalTrackSite*> _hit_used_for_sites{};
+
+ /** map to store relation between lcio hits kaltest hits
+ */
+ std::map<edm4hep::TrackerHit,DDVTrackHit*> _lcio_hits_to_kaltest_hits{};
+
+ /** vector to store lcio hits rejected for measurement sites
+ */
+ std::vector<edm4hep::TrackerHit> _hit_not_used_for_sites{};
+
+ /** vector to store the chi-sqaure increment for measurement sites
+ */
+ std::vector< std::pair<edm4hep::TrackerHit, double> > _hit_chi2_values{};
+
+ /** vector to store the chi-sqaure increment for measurement sites
+ */
+ std::vector< std::pair<edm4hep::TrackerHit, double> > _outlier_chi2_values{};
+
+
+} ;
+
+} // end of namespace MarlinTrk
+
+#endif
diff --git a/Service/TrackSystemSvc/src/TrackSystemSvc.cpp b/Service/TrackSystemSvc/src/TrackSystemSvc.cpp
index 31d5a5b7757e94d38036317d3c44422d763ee5b5..723a39551cf0df3e264980d3c78110c362ec6215 100644
--- a/Service/TrackSystemSvc/src/TrackSystemSvc.cpp
+++ b/Service/TrackSystemSvc/src/TrackSystemSvc.cpp
@@ -2,6 +2,7 @@
#include "gear/GearMgr.h"
#include "MarlinKalTest.h"
+#include "MarlinDDKalTest.h"
#include "TrackSystemSvc.h"
@@ -14,31 +15,44 @@ TrackSystemSvc::TrackSystemSvc(const std::string& name, ISvcLocator* svc)
TrackSystemSvc::~TrackSystemSvc(){
}
-MarlinTrk::IMarlinTrkSystem* TrackSystemSvc::getTrackSystem(void* address){
+MarlinTrk::IMarlinTrkSystem* TrackSystemSvc::getTrackSystem(void* address, std::string type) {
std::map<void*, MarlinTrk::IMarlinTrkSystem*>::iterator it=m_trackSystems.find(address);
if(it==m_trackSystems.end()){
- gear::GearMgr* mgr=0;
- auto _gear = service<IGearSvc>("GearSvc");
- if ( !_gear ) {
- info() << "Failed to find GearSvc ..." << endmsg;
+ MarlinTrk::IMarlinTrkSystem* sys = nullptr;
+
+ if (type=="KalTest") {
+ gear::GearMgr* mgr=0;
+ auto _gear = service<IGearSvc>("GearSvc");
+ if ( !_gear ) {
+ info() << "Failed to find GearSvc ..." << endmsg;
+ }
+ else{
+ mgr = _gear->getGearMgr();
+ }
+
+ auto _geoSvc = service<IGeomSvc>("GeomSvc");
+ if ( !_geoSvc ) {
+ info() << "Failed to find GeomSvc ..." << endmsg;
+ }
+ if(mgr==0&&_geoSvc==0){
+ fatal() << "Both GearSvc and GeomSvc invalid!" << endmsg;
+ return 0;
+ }
+ debug() << "GearMgr=" << mgr << " GeomSvc=" << _geoSvc << endmsg;
+ //sys = new MarlinTrk::MarlinKalTest( *mgr, _geoSvc );
+ sys = new MarlinTrk::MarlinKalTest(*mgr);
}
- else{
- mgr = _gear->getGearMgr();
+ else if (type=="DDKalTest") {
+ sys = new MarlinTrk::MarlinDDKalTest();
}
-
- auto _geoSvc = service<IGeomSvc>("GeomSvc");
- if ( !_geoSvc ) {
- info() << "Failed to find GeomSvc ..." << endmsg;
+ else {
+ error() << type << " not support, check type" << endmsg;
}
- if(mgr==0&&_geoSvc==0){
- fatal() << "Both GearSvc and GeomSvc invalid!" << endmsg;
- return 0;
+
+ if (sys) {
+ m_trackSystems[address] = sys;
+ debug() << "Track system created successfully for " << address << endmsg;
}
- debug() << "GearMgr=" << mgr << " GeomSvc=" << _geoSvc << endmsg;
- //MarlinTrk::IMarlinTrkSystem* sys = new MarlinTrk::MarlinKalTest( *mgr, _geoSvc );
- MarlinTrk::IMarlinTrkSystem* sys = new MarlinTrk::MarlinKalTest(*mgr);
- m_trackSystems[address] = sys;
- debug() << "Track system created successfully for " << address << endmsg;
return sys;
}
return it->second;
diff --git a/Service/TrackSystemSvc/src/TrackSystemSvc.h b/Service/TrackSystemSvc/src/TrackSystemSvc.h
index b03dbd9b26d74f16c954206867ab976a26e39398..9d4928399fcf700f40ab2a0e27c9475130520d59 100644
--- a/Service/TrackSystemSvc/src/TrackSystemSvc.h
+++ b/Service/TrackSystemSvc/src/TrackSystemSvc.h
@@ -9,7 +9,7 @@ class TrackSystemSvc : public extends<Service, ITrackSystemSvc>{
TrackSystemSvc(const std::string& name, ISvcLocator* svc);
~TrackSystemSvc();
- MarlinTrk::IMarlinTrkSystem* getTrackSystem(void* address=0) override;
+ MarlinTrk::IMarlinTrkSystem* getTrackSystem(void* address=0, std::string type="KalTest") override;
void removeTrackSystem(void* address=0) override;
StatusCode initialize() override;
diff --git a/Utilities/CMakeLists.txt b/Utilities/CMakeLists.txt
index fd5eb932fd7ba2bdce7640d104ca0c8ea8eac86d..856908482e98e6f30cdbac2d0469d385ed26c046 100644
--- a/Utilities/CMakeLists.txt
+++ b/Utilities/CMakeLists.txt
@@ -1,5 +1,6 @@
add_subdirectory(DataHelper)
-add_subdirectory(KalTest)
+#obsolete, import through FetchContent now
+#add_subdirectory(KalTest)
add_subdirectory(KalDet)
add_subdirectory(KiTrack)
add_subdirectory(DecoderHelper)
diff --git a/Utilities/KalDet/CMakeLists.txt b/Utilities/KalDet/CMakeLists.txt
index 8511f2378df384b1241510f672208cc56e625967..61d4387f3e14b3a70523db59d3f92eedaf0d76d7 100644
--- a/Utilities/KalDet/CMakeLists.txt
+++ b/Utilities/KalDet/CMakeLists.txt
@@ -2,13 +2,13 @@
# Package: KalDet
# Desc: import from ILCSoft
##############################################################################
-
-get_target_property(to_incl KalTestLib SOURCE_DIR)
-if (to_incl)
- LIST( APPEND DICT_INCLUDE_DIRS ${to_incl}/include)
-else()
- message(FATAL_ERROR "Failed to get the source dir for package KalTestLib")
-endif()
+#get_target_property(to_incl KalTestLib SOURCE_DIR)
+#if (to_incl)
+# LIST( APPEND DICT_INCLUDE_DIRS ${to_incl}/include)
+#else()
+# message(FATAL_ERROR "Failed to get the source dir for package KalTestLib")
+#endif()
+LIST( APPEND DICT_INCLUDE_DIRS ${KalTest_INCLUDE_DIRS})
set( DICT_CINT_DEFINITIONS "HANDLE_DICT_EXCEPTIONS=IGNORED_FOR_CINT" )
set( DICT_INPUT_DIRS gen kern lctpc/gearTPC )
@@ -65,7 +65,7 @@ gaudi_add_library(KalDetLib
SOURCES ${KalDetLib_srcs}
LINK DetInterface
DetIdentifier
- KalTestLib
+ ${KalTest_LIBRARIES}
Gaudi::GaudiKernel
${ROOT_LIBRARIES}
${CLHEP_LIBRARIES}
@@ -82,6 +82,8 @@ target_include_directories(KalDetLib
PUBLIC ${LCIO_INCLUDE_DIRS}
)
+add_dependencies(KalDetLib KalTest)
+
install(TARGETS KalDetLib
EXPORT CEPCSWTargets
RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
diff --git a/cmake/CEPCSWDependencies.cmake b/cmake/CEPCSWDependencies.cmake
index 2778a89e7591e711ba3f697767adcf819c6338d4..dec2acd8d7eeea35d62c0db95bec42f29ffbc44a 100644
--- a/cmake/CEPCSWDependencies.cmake
+++ b/cmake/CEPCSWDependencies.cmake
@@ -18,6 +18,10 @@ Find all the dependencies here, so in each package user don't need to find the p
- podio
- ROOT
- CKF
+- ILCUTIL
+- KalTest
+- aidaTT
+- DDKalTest
- RDAnalysis
#]]
@@ -56,6 +60,37 @@ else()
include("${CMAKE_CURRENT_LIST_DIR}/internal_ckf.cmake")
endif()
+if (CEPCSW_USE_SYSTEM_ILCUTIL)
+ message("Try to use an existing installation of ILCUTIL")
+ find_package(ILCUTIL)
+else()
+ message("Try to use an internal installation of ILCUTIL")
+ include("${CMAKE_CURRENT_LIST_DIR}/internal_ilcutil.cmake")
+endif()
+
+if (CEPCSW_USE_SYSTEM_KALTEST)
+ message("Try to use an existing installation of KalTest")
+ find_package(KalTest)
+else()
+ message("Try to use an internal installation of KalTest")
+ include("${CMAKE_CURRENT_LIST_DIR}/internal_kaltest.cmake")
+endif()
+
+if (CEPCSW_USE_SYSTEM_AIDATT)
+ message("Try to use an existing installation of aidaTT")
+ find_package(aidaTT)
+else()
+ message("Try to use an internal installation of aidaTT")
+ include("${CMAKE_CURRENT_LIST_DIR}/internal_aidatt.cmake")
+endif()
+
+if (CEPCSW_USE_SYSTEM_DDKALTEST)
+ message("Try to use an existing installation of DDKalTest")
+ find_package(DDKalTest)
+else()
+ message("Try to use an internal installation of DDKalTest")
+ include("${CMAKE_CURRENT_LIST_DIR}/internal_ddkaltest.cmake")
+endif()
if (CEPCSW_USE_SYSTEM_EDM4CEPC)
message("Try to use an existing installation of EDM4CEPC")
diff --git a/cmake/CEPCSWOptions.cmake b/cmake/CEPCSWOptions.cmake
index a853022fe2e1b1c98683eb56177849d632074a1f..5dca0ca413d83b3f5520004947e2c5f68e521df3 100644
--- a/cmake/CEPCSWOptions.cmake
+++ b/cmake/CEPCSWOptions.cmake
@@ -2,6 +2,17 @@
##############################################################################
## CKF Belle
+## EDM4CEPC
+## ILCUTIL
+## KalTest
+## aidaTT
+## DDKalTest
+##
+## CEPCSW_EXTERNAL_AS_COMPONENT
+## TRUE install external project as component of CEPCSW in directory same
+## as CEPCSW
+## FALSE install external project as external library in their speical
+## directory, optional through CEPCSW_EXTERNAL_ROOT
##############################################################################
option(CEPCSW_USE_SYSTEM_CKF_BELLE
@@ -11,3 +22,29 @@ option(CEPCSW_USE_SYSTEM_CKF_BELLE
option(CEPCSW_USER_SYSTEM_EDM4CEPC
"Use the existing installation of EDM4CEPC. Otherwise the internal version will be used."
FALSE)
+
+option(CEPCSW_USER_SYSTEM_ILCUTIL
+ "Use the existing installation of ILCUTIL. Otherwise the internal version will be used."
+ FALSE)
+
+option(CEPCSW_USER_SYSTEM_AIDATT
+ "Use the existing installation of aidaTT. Otherwise the internal version will be used."
+ FALSE)
+
+option(CEPCSW_USER_SYSTEM_KALTEST
+ "Use the existing installation of KalTest. Otherwise the internal version will be used."
+ FALSE)
+
+option(CEPCSW_USER_SYSTEM_DDKALTEST
+ "Use the existing installation of DDKalTest. Otherwise the internal version will be used."
+ FALSE)
+
+option(CEPCSW_EXTERNAL_AS_COMPONENT
+ "Use special option to install external into CEPCSW install directory."
+ TRUE)
+
+if (NOT CEPCSW_EXTERNAL_AS_COMPONENT)
+ if (NOT DEFINED CEPCSW_EXTERNAL_ROOT)
+ set(CEPCSW_EXTERNAL_ROOT ${CMAKE_INSTALL_PREFIX}/External)
+ endif()
+endif()
diff --git a/cmake/internal_aidatt.cmake b/cmake/internal_aidatt.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..c1faf6c5d7d3287ba0e3eca8081f9b5f48d83b82
--- /dev/null
+++ b/cmake/internal_aidatt.cmake
@@ -0,0 +1,26 @@
+include(ExternalProject)
+
+if (CEPCSW_EXTERNAL_AS_COMPONENT)
+ set( aidaTT_DIR ${CMAKE_INSTALL_PREFIX} )
+else()
+ set( aidaTT_DIR ${CEPCSW_EXTERNAL_ROOT}/aidaTT )
+endif()
+
+ExternalProject_Add(
+ aidaTT
+ GIT_REPOSITORY https://code.ihep.ac.cn/cepc/externals/mirroring/aidaTT.git
+ GIT_TAG v00-10-cepcsw
+ PREFIX ${CMAKE_BINARY_DIR}/_deps
+ SOURCE_DIR ${CMAKE_BINARY_DIR}/_deps/aidatt-src
+ BINARY_DIR ${CMAKE_BINARY_DIR}/_deps/aidatt-build
+ STAMP_DIR ${CMAKE_BINARY_DIR}/_deps/aidatt-stamp
+ CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${aidaTT_DIR}
+ -DCMAKE_LIBRARY_OUTPUT_DIRECTORY=${CMAKE_BINARY_DIR}/lib
+ -DCMAKE_PREFIX_PATH=${ILCUTIL_DIR}/lib/cmake/ILCUTIL
+ -DCONFIGURE_FILE_INSTALL_DIR=lib/cmake/aidaTT
+ -DUSE_CXX11=FALSE
+ UPDATE_COMMAND ""
+ DEPENDS ILCUTIL
+)
+
+message(STATUS "will install into ${aidaTT_DIR}")
diff --git a/cmake/internal_ddkaltest.cmake b/cmake/internal_ddkaltest.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..ee442f443b636b889b57e23ca51d0ee339151d81
--- /dev/null
+++ b/cmake/internal_ddkaltest.cmake
@@ -0,0 +1,33 @@
+include(ExternalProject)
+
+if (CEPCSW_EXTERNAL_AS_COMPONENT)
+ set(DDKalTest_DIR ${CMAKE_INSTALL_PREFIX})
+else()
+ set(DDKalTest_DIR ${CEPCSW_EXTERNAL_ROOT}/DDKalTest)
+endif()
+
+ExternalProject_Add(
+ DDKalTest
+ GIT_REPOSITORY https://code.ihep.ac.cn/cepc/externals/DDKalTest_iLCSoft.git
+ GIT_TAG v01-07-cepcsw
+ PREFIX ${CMAKE_BINARY_DIR}/_deps
+ SOURCE_DIR ${CMAKE_BINARY_DIR}/_deps/ddkaltest-src
+ BINARY_DIR ${CMAKE_BINARY_DIR}/_deps/ddkaltest-build
+ STAMP_DIR ${CMAKE_BINARY_DIR}/_deps/ddkaltest-stamp
+ UPDATE_COMMAND ""
+ CONFIGURE_COMMAND ${CMAKE_COMMAND}
+ -DCMAKE_INSTALL_PREFIX=${DDKalTest_DIR}
+ -DCMAKE_LIBRARY_OUTPUT_DIRECTORY=${CMAKE_BINARY_DIR}/lib
+ -DCONFIGURE_FILE_INSTALL_DIR=lib/cmake/DDKalTest
+ -DKalTest_DIR=${KalTest_DIR}/lib/cmake/KalTest
+ -DILCUTIL_DIR=${ILCUTIL_DIR}/lib/cmake/ILCUTIL
+ -DaidaTT_DIR=${aidaTT_DIR}/lib/cmake/aidaTT
+ <SOURCE_DIR>
+ BUILD_COMMAND ${CMAKE_COMMAND} --build . --target install
+ DEPENDS ILCUTIL KalTest aidaTT
+)
+
+message(STATUS "will install into ${DDKalTest_DIR}")
+
+set(DDKalTest_INCLUDE_DIRS ${DDKalTest_DIR}/include)
+set(DDKalTest_LIBRARIES ${DDKalTest_DIR}/lib/libDDKalTest.so)
diff --git a/cmake/internal_ilcutil.cmake b/cmake/internal_ilcutil.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..fe2d512cb221b2ed5aa15ac7037d77ab868845f5
--- /dev/null
+++ b/cmake/internal_ilcutil.cmake
@@ -0,0 +1,27 @@
+include(ExternalProject)
+
+if (CEPCSW_EXTERNAL_AS_COMPONENT)
+ set( ILCUTIL_DIR ${CMAKE_INSTALL_PREFIX} )
+else()
+ set( ILCUTIL_DIR ${CEPCSW_EXTERNAL_ROOT}/ILCUTIL )
+endif()
+
+ExternalProject_Add(
+ ILCUTIL
+ GIT_REPOSITORY https://code.ihep.ac.cn/cepc/externals/mirroring/iLCUtil.git
+ GIT_TAG v01-07-03-cepcsw
+ PREFIX ${CMAKE_BINARY_DIR}/_deps
+ SOURCE_DIR ${CMAKE_BINARY_DIR}/_deps/ilcutil-src
+ BINARY_DIR ${CMAKE_BINARY_DIR}/_deps/ilcutil-build
+ STAMP_DIR ${CMAKE_BINARY_DIR}/_deps/ilcutil-stamp
+ CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${ILCUTIL_DIR}
+ -DCMAKE_LIBRARY_OUTPUT_DIRECTORY=${CMAKE_BINARY_DIR}/lib
+ -DCONFIGURE_FILE_INSTALL_DIR=lib/cmake/ILCUTIL
+ -DINSTALL_DOC=OFF
+ UPDATE_COMMAND ""
+)
+
+message(STATUS "will install into ${ILCUTIL_DIR}")
+
+set(streamlog_INCLUDE_DIRS ${ILCUTIL_DIR}/include)
+set(streamlog_LIBRARIES ${ILCUTIL_DIR}/lib/libstreamlog.so)
diff --git a/cmake/internal_kaltest.cmake b/cmake/internal_kaltest.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..c8d6764de4a32e89b5d00a15c738675ec5d9ceea
--- /dev/null
+++ b/cmake/internal_kaltest.cmake
@@ -0,0 +1,34 @@
+include(ExternalProject)
+
+OPTION( BUILD_WITH_T0_FIT "Set to ON to build with t0 fit (kSdim=6)" ON )
+
+if (CEPCSW_EXTERNAL_AS_COMPONENT)
+ set( KalTest_DIR ${CMAKE_INSTALL_PREFIX} )
+else()
+ set( KalTest_DIR ${CEPCSW_EXTERNAL_ROOT}/KalTest )
+endif()
+
+ExternalProject_Add(
+ KalTest
+ GIT_REPOSITORY https://code.ihep.ac.cn/cepc/externals/KalTest_iLCSoft.git
+ GIT_TAG v02-05-cepcsw
+ PREFIX ${CMAKE_BINARY_DIR}/_deps
+ SOURCE_DIR ${CMAKE_BINARY_DIR}/_deps/kaltest-src
+ BINARY_DIR ${CMAKE_BINARY_DIR}/_deps/kaltest-build
+ STAMP_DIR ${CMAKE_BINARY_DIR}/_deps/kaltest-stamp
+ UPDATE_COMMAND ""
+ CONFIGURE_COMMAND ${CMAKE_COMMAND}
+ -DCMAKE_INSTALL_PREFIX=${KalTest_DIR}
+ -DCMAKE_LIBRARY_OUTPUT_DIRECTORY=${CMAKE_BINARY_DIR}/lib
+ -DCMAKE_PREFIX_PATH=${ILCUTIL_DIR}/lib/cmake/ILCUTIL
+ -DCONFIGURE_FILE_INSTALL_DIR=lib/cmake/KalTest
+ -DBUILD_WITH_T0_FIT=${BUILD_WITH_T0_FIT}
+ <SOURCE_DIR>
+ BUILD_COMMAND ${CMAKE_COMMAND} --build . --target install
+ DEPENDS ILCUTIL
+)
+
+message(STATUS "will install into ${KalTest_DIR}")
+
+set(KalTest_INCLUDE_DIRS ${KalTest_DIR}/include ${KalTest_DIR}/include/kaltest)
+set(KalTest_LIBRARIES ${KalTest_DIR}/lib/libKalTest.so)