From a8ca1ffae80c79de0089b594bd3e5921b527059c Mon Sep 17 00:00:00 2001
From: Chengdong Fu <fucd@ihep.ac.cn>
Date: Mon, 1 Apr 2024 14:15:56 +0800
Subject: [PATCH] change to ITrackFitterTool
---
Reconstruction/SiliconTracking/CMakeLists.txt | 1 +
.../src/ForwardTrackingAlg.cpp | 17 +-
.../SiliconTracking/src/ForwardTrackingAlg.h | 4 +
.../src/SiliconTrackingAlg.cpp | 237 ++++++++++++------
.../SiliconTracking/src/SiliconTrackingAlg.h | 16 +-
.../SiliconTracking/src/TrackSubsetAlg.cpp | 29 ++-
.../SiliconTracking/src/TrackSubsetAlg.h | 4 +
.../Tracking/src/Clupatra/ClupatraAlg.cpp | 3 +
.../Tracking/src/Clupatra/ClupatraAlg.h | 3 +
.../FullLDCTracking/FullLDCTrackingAlg.cpp | 37 ++-
.../src/FullLDCTracking/FullLDCTrackingAlg.h | 3 +
Service/GearSvc/src/GearSvc.h | 3 +-
12 files changed, 256 insertions(+), 101 deletions(-)
mode change 100644 => 100755 Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
diff --git a/Reconstruction/SiliconTracking/CMakeLists.txt b/Reconstruction/SiliconTracking/CMakeLists.txt
index 9eb87350..8fd9be94 100644
--- a/Reconstruction/SiliconTracking/CMakeLists.txt
+++ b/Reconstruction/SiliconTracking/CMakeLists.txt
@@ -7,6 +7,7 @@ gaudi_add_module(SiliconTracking
src/TrackSubsetAlg.cpp
LINK GearSvc
EventSeeder
+ TrackingLib
TrackSystemSvcLib
DataHelperLib
KiTrackLib
diff --git a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp
index 3eb06e29..3184563f 100644
--- a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp
+++ b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp
@@ -4,7 +4,7 @@
#include "DataHelper/Navigation.h"
#include "edm4hep/TrackerHit.h"
-#include "edm4hep/TrackerHit.h"
+#include "edm4hep/TrackerHitPlane.h"
#include "edm4hep/Track.h"
#if __has_include("edm4hep/EDM4hepVersion.h")
#include "edm4hep/EDM4hepVersion.h"
@@ -74,12 +74,16 @@ StatusCode ForwardTrackingAlg::initialize(){
// can be printed. As this is mainly used for debugging it is not a steerable parameter.
//if( _useCED )MarlinCED::init(this) ; //CED
+ // Set up the track fit tool
+ m_fitTool = ToolHandle<ITrackFitterTool>(m_fitToolName.value());
+
if(m_dumpTime){
NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name());
if ( !nt1 ) {
m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time");
if ( 0 != m_tuple ) {
m_tuple->addItem ("timeTotal", m_timeTotal ).ignore();
+ m_tuple->addItem ("timeKalman", m_timeKalman ).ignore();
}
else {
fatal() << "Cannot book MyTuples/Time"+name() <<endmsg;
@@ -306,7 +310,9 @@ StatusCode ForwardTrackingAlg::execute(){
_map_sector_hits[ ftdHit->getSector() ].push_back( ftdHit );
}
}
-
+
+ if(m_dumpTime&&m_tuple) m_timeKalman = 0;
+
if( !_map_sector_hits.empty() ){
/**********************************************************************************************/
/* Check if no sector is overflowing with hits */
@@ -661,7 +667,7 @@ StatusCode ForwardTrackingAlg::execute(){
debug() << "\t\t---Save Tracks---" << endmsg ;
//auto trkCol = _outColHdl.createAndPut();
-
+
for (unsigned int i=0; i < tracks.size(); i++){
FTDTrack* myTrack = dynamic_cast< FTDTrack* >( tracks[i] );
@@ -942,7 +948,7 @@ bool ForwardTrackingAlg::setCriteria( unsigned round ){
}
void ForwardTrackingAlg::finaliseTrack( edm4hep::MutableTrack* trackImpl ){
-
+ auto stopwatch = TStopwatch();
Fitter fitter( trackImpl , _trkSystem );
//trackImpl->trackStates().clear();
@@ -1024,6 +1030,9 @@ void ForwardTrackingAlg::finaliseTrack( edm4hep::MutableTrack* trackImpl ){
trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SET]);
trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::ETD]);
#endif
+
+ if(m_dumpTime&&m_tuple) m_timeKalman += stopwatch.RealTime()*1000;
+
return;
}
diff --git a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h
index bd7a5173..fee28994 100644
--- a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h
+++ b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h
@@ -14,6 +14,7 @@
#include "edm4hep/TrackerHit.h"
#include "TrackSystemSvc/IMarlinTrkSystem.h"
+#include "Tracking/ITrackFitterTool.h"
//#include "gear/BField.h"
#include "KiTrack/Segment.h"
@@ -239,6 +240,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm {
NTuple::Tuple* m_tuple;
NTuple::Item<float> m_timeTotal;
+ NTuple::Item<float> m_timeKalman;
/** B field in z direction */
double _Bz;
@@ -261,6 +263,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm {
Gaudi::Property<std::vector<float> > _critMinimaInit{this, "CriteriaMin", {} };
Gaudi::Property<std::vector<float> > _critMaximaInit{this, "CriteriaMax", {} };
Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true};
+ Gaudi::Property<std::string> m_fitToolName{this, "FitterTool", "KalTestTool/KalTest110"};
std::map<std::string, std::vector<float> > _critMinima;
std::map<std::string, std::vector<float> > _critMaxima;
@@ -291,6 +294,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm {
unsigned _nTrackCandidatesPlus;
MarlinTrk::IMarlinTrkSystem* _trkSystem;
+ ToolHandle<ITrackFitterTool> m_fitTool;
/** The quality of the output track collection */
int _output_track_col_quality ;
diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp
index 37110a92..c1068912 100644
--- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp
+++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp
@@ -113,13 +113,17 @@ StatusCode SiliconTrackingAlg::initialize() {
_nRun = -1 ;
_nEvt = 0 ;
//printParameters() ;
- if(m_dumpTime){
+
+ // Set up the track fit tool
+ m_fitTool = ToolHandle<ITrackFitterTool>(m_fitToolName.value());
+
+ if (m_dumpTime) {
NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name());
if ( !nt1 ) {
m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time");
if ( 0 != m_tuple ) {
- m_tuple->addItem ("timeTotal", m_timeTotal ).ignore();
- m_tuple->addItem ("timeKalman", m_timeKalman ).ignore();
+ m_tuple->addItem("timeTotal", m_timeTotal).ignore();
+ m_tuple->addItem("timeKalman", m_timeKalman).ignore();
}
else {
fatal() << "Cannot book MyTuples/Time"+name() <<endmsg;
@@ -130,6 +134,27 @@ StatusCode SiliconTrackingAlg::initialize() {
m_tuple = nt1;
}
}
+
+ if (m_debug) {
+ NTuplePtr nt2(ntupleSvc(), "MyTuples/"+name());
+ if ( !nt2 ) {
+ m_tupleDebug = ntupleSvc()->book("MyTuples/"+name(),CLID_ColumnWiseTuple,"Tracking time");
+ if ( 0 != m_tuple ) {
+ m_tupleDebug->addItem("nhit", m_nhit, 0, 100).ignore();
+ m_tupleDebug->addIndexedItem("layer", m_nhit, m_layer).ignore();
+ m_tupleDebug->addIndexedItem("chi2", m_nhit, m_chi2).ignore();
+ m_tupleDebug->addIndexedItem("used", m_nhit, m_used).ignore();
+ }
+ else {
+ fatal() << "Cannot book MyTuples/"+name() <<endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ else{
+ m_tupleDebug = nt2;
+ }
+ }
+
// set up the geometery needed by KalTest
//FIXME: for now do KalTest only - make this a steering parameter to use other fitters
auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc");
@@ -192,7 +217,7 @@ StatusCode SiliconTrackingAlg::execute(){
// zero triplet counters
_ntriplets = _ntriplets_good = _ntriplets_2MCP = _ntriplets_3MCP = _ntriplets_1MCP_Bad = _ntriplets_bad = 0;
-
+
// Clearing the working containers from the previous event
// FIXME: partly done at the end of the event, in CleanUp. Make it consistent.
//_tracksWithNHitsContainer.clear();
@@ -205,14 +230,14 @@ StatusCode SiliconTrackingAlg::execute(){
//int runNo = header.getRunNumber();
//debug() << "Processing Run[" << runNo << "]::Event[" << evtNo << "]" << endmsg;
- _trackImplVec.reserve(100);
+ //_trackImplVec.reserve(100);
int successVTX = InitialiseVTX();
//int successFTD = 0;
int successFTD = InitialiseFTD();
if (successVTX == 1) {
- debug() << " phi theta layer nh o : m : i :: o*m*i " << endmsg;
+ debug() << " --phi-- --theta-- --layer-- nh o : m : i :: o*m*i " << endmsg;
for (int iPhi=0; iPhi<_nDivisionsInPhi; ++iPhi) {
for (int iTheta=0; iTheta<_nDivisionsInTheta;++iTheta) {
@@ -225,7 +250,7 @@ StatusCode SiliconTrackingAlg::execute(){
}
if (successFTD == 1) {
- debug() << " phi side layer nh o : m : i :: o*m*i " << endmsg;
+ debug() << " --phi-- --side-- --layer-- nh o : m : i :: o*m*i " << endmsg;
TrackingInFTD(); // Perform tracking in the FTD
debug() << "End of Processing FTD sectors" << endmsg;
}
@@ -250,7 +275,17 @@ StatusCode SiliconTrackingAlg::execute(){
trackIter < tracksWithNHits.end(); trackIter++) {
CreateTrack( *trackIter );
}
- debug() << "End of creating "<< nHits << " hits tracks " << endmsg;
+ debug() << "End of creating "<< nHits << " hits tracks, total " << _trackImplVec.size() << " tracks now" << endmsg;
+
+ if (msgLevel(MSG::DEBUG)) {
+ for (auto trackAR : _trackImplVec) {
+ debug() << "track " << trackAR << " has hits:" << endmsg;
+ TrackerHitExtendedVec& hitVec = trackAR->getTrackerHitExtendedVec();
+ for (auto hit : hitVec) {
+ debug() << " " << hit->getTrackerHit().id().collectionID << "-" << hit->getTrackerHit().id().index << endmsg;
+ }
+ }
+ }
}
if (_attachFast == 0) {
@@ -539,7 +574,7 @@ int SiliconTrackingAlg::InitialiseFTD() {
int iCode = iSemiSphere + 2*layer + 2*_nlayersFTD*iPhi;
_sectorsFTD[iCode].push_back( hitExt );
- debug() << " FTD Pixel Hit added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iSemiSphere " << iSemiSphere << " iCode = " << iCode << " layer = " << layer << " cellID = " << hit.getCellID() << endmsg;
+ debug() << " FTD Pixel Hit " << hit.id().index << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iSemiSphere " << iSemiSphere << " iCode = " << iCode << " layer = " << layer << " cellID = " << hit.getCellID() << endmsg;
}
}
// Reading out FTD SpacePoint Collection
@@ -603,7 +638,7 @@ int SiliconTrackingAlg::InitialiseFTD() {
double Phi = atan2(pos[1],pos[0]);
if (Phi < 0.) Phi = Phi + TWOPI;
-
+
// get the layer number
unsigned int layer = static_cast<unsigned int>(getLayerID(hit));
unsigned int petalIndex = static_cast<unsigned int>(getModuleID(hit));
@@ -636,7 +671,7 @@ int SiliconTrackingAlg::InitialiseFTD() {
int iCode = iSemiSphere + 2*layer + 2*_nlayersFTD*iPhi;
_sectorsFTD[iCode].push_back( hitExt );
- debug() << " FTD SpacePoint Hit added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iSemiSphere " << iSemiSphere << " iCode = " << iCode << " layer = " << layer << " cellID = " << hit.getCellID() << endmsg;
+ debug() << " FTD SpacePoint Hit added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iSemiSphere " << iSemiSphere << " iCode = " << iCode << " layer = " << layer << " cellID = " << hit.getCellID() << endmsg;
}
}
@@ -713,7 +748,7 @@ int SiliconTrackingAlg::InitialiseVTX() {
// SJA:FIXME: just use planar res for now
hitExt->setResolutionRPhi(hit.getCovMatrix()[2]);
hitExt->setResolutionZ(hit.getCovMatrix()[5]);
-
+
// set type is now only used in one place where it is set to 0 to reject hits from a fit, set to INT_MAX to try and catch any missuse
hitExt->setType(int(INT_MAX));
// det is no longer used set to INT_MAX to try and catch any missuse
@@ -733,7 +768,7 @@ int SiliconTrackingAlg::InitialiseVTX() {
double Phi = atan2(pos[1],pos[0]);
if (Phi < 0.) Phi = Phi + TWOPI;
-
+
// get the layer number
int layer = getLayerID(hit);
@@ -747,7 +782,7 @@ int SiliconTrackingAlg::InitialiseVTX() {
int iCode = layer + _nLayers*iPhi + _nLayers*_nDivisionsInPhi*iTheta;
_sectors[iCode].push_back( hitExt );
- debug() << " VXD Hit " << hit.id() << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iTheta " << iTheta << " iCode = " << iCode << " layer = " << layer << endmsg;
+ debug() << " VXD Hit " << hit.id().index << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iTheta " << iTheta << " iCode = " << iCode << " layer = " << layer << endmsg;
}
}
@@ -906,7 +941,7 @@ int SiliconTrackingAlg::InitialiseVTX() {
int iCode = layer + _nLayers*iPhi + _nLayers*_nDivisionsInPhi*iTheta;
_sectors[iCode].push_back( hitExt );
- debug() << " SIT Hit " << trkhit.id() << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iTheta " << iTheta << " iCode = " << iCode << " layer = " << layer << endmsg;
+ debug() << " SIT Hit " << trkhit.id().index << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iTheta " << iTheta << " iCode = " << iCode << " layer = " << layer << endmsg;
}
}
@@ -1063,13 +1098,13 @@ void SiliconTrackingAlg::ProcessOneSector(int iPhi, int iTheta) {
if (nHitsInner > 0) {
- debug() << " "
- << std::setw(3) << iPhi << " " << std::setw(3) << ipMiddle << " " << std::setw(3) << ipInner << " "
- << std::setw(3) << iTheta << " " << std::setw(3) << itMiddle << " " << std::setw(3) << itInner << " "
- << std::setw(3) << nLR[0] << " " << std::setw(3) << nLR[1] << " " << std::setw(3) << nLR[2] << " "
- << std::setw(3) << nHitsOuter << " : " << std::setw(3) << nHitsMiddle << " : " << std::setw(3) << nHitsInner << " :: "
- << std::setw(3) << nHitsOuter*nHitsMiddle* nHitsInner << endmsg;
-
+ debug() << "pattern "
+ << std::setw(3) << iPhi << " " << std::setw(3) << ipMiddle << " " << std::setw(3) << ipInner << " "
+ << std::setw(3) << iTheta << " " << std::setw(3) << itMiddle << " " << std::setw(3) << itInner << " "
+ << std::setw(3) << nLR[0] << " " << std::setw(3) << nLR[1] << " " << std::setw(3) << nLR[2] << " "
+ << std::setw(3) << nHitsOuter << " : " << std::setw(3) << nHitsMiddle << " : " << std::setw(3) << nHitsInner << " :: "
+ << std::setw(3) << nHitsOuter*nHitsMiddle*nHitsInner << endmsg;
+
// test all triplets
for (int iOuter=0; iOuter<nHitsOuter; ++iOuter) { // loop over hits in the outer sector
@@ -1219,6 +1254,8 @@ TrackExtended * SiliconTrackingAlg::TestTriplet(TrackerHitExtended * outerHit,
float d0 = par[3];
float z0 = par[4];
+ debug() << "after fastHelixFit: chi2RPhi = " << chi2RPhi << " chi2Z = " << chi2Z
+ << " tanlambda = " << tanlambda << " coslambda = " << cos(atan(tanlambda)) << endmsg;
// chi2 is weighted here by a factor for both rphi and z
float Chi2 = chi2RPhi*_chi2WRPhiTriplet+chi2Z*_chi2WZTriplet;
int ndf = 2*NPT-5;
@@ -1476,6 +1513,7 @@ int SiliconTrackingAlg::BuildTrack(TrackerHitExtended * outerHit,
*/
debug() << " BuildTrack starting " << endmsg;
+ debug() << outerHit->getTrackerHit().id().index << " " << middleHit->getTrackerHit().id().index << " " << innerHit->getTrackerHit().id().index << " innerLayer=" << innerLayer << endmsg;
for (int layer = innerLayer-1; layer>=0; layer--) { // loop over remaining layers
float distMin = 1.0e+20;
@@ -1571,7 +1609,7 @@ int SiliconTrackingAlg::BuildTrack(TrackerHitExtended * outerHit,
float chi2RPhi;
float chi2Z;
- //debug() << "######## number of hits to fit with _fastfitter = " << NPT << endmsg;
+ //debug() << "######## number of hits to fit with _fastfitter = " << NPT << endmsg;
_fastfitter->fastHelixFit(NPT, xh, yh, rh, ph, wrh, zh, wzh,iopt, par, epar, chi2RPhi, chi2Z);
par[3] = par[3]*par[0]/fabs(par[0]);
@@ -1623,7 +1661,8 @@ int SiliconTrackingAlg::BuildTrack(TrackerHitExtended * outerHit,
} // endloop over remaining layers
TrackerHitExtendedVec& hvec = trackAR->getTrackerHitExtendedVec();
int nTotalHits = int(hvec.size());
- debug() << "######## number of hits to return = " << nTotalHits << endmsg;
+ debug() << "######## number of hits to return = " << nTotalHits << endmsg;
+
return nTotalHits;
}
@@ -1654,15 +1693,20 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) {
Method which creates Track out of TrackExtended objects. Checks for possible
track splitting (separate track segments in VXD and FTD).
*/
-
+ debug() << "CreateTrack " << trackAR << endmsg;
TrackerHitExtendedVec& hitVec = trackAR->getTrackerHitExtendedVec();
int nHits = int(hitVec.size());
for (int i=0; i<nHits; ++i) {
TrackExtendedVec& trackVec = hitVec[i]->getTrackExtendedVec();
- if (trackVec.size() != 0)
+ //edm4hep::TrackerHit trkHit = hitVec[i]->getTrackerHit();
+ //debug() << "Hit " << trkHit.id().collectionID << "-" << trkHit.id().index << endmsg;
+ if (trackVec.size() != 0) {
+ //debug() << " used in TrackExtended " << endmsg;
return ;
+ //continue;
+ }
}
// First check if the current track is piece of the split one
@@ -1674,7 +1718,6 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) {
for (int itrk=0; itrk<nTrk; ++itrk) {
TrackExtended * trackOld = _trackImplVec[itrk];
- // fucd: TrackerHitExtendedVec& will change after merge split tracks, so must TrackerHitExtendedVec
TrackerHitExtendedVec hitVecOld = trackOld->getTrackerHitExtendedVec();
float phiNew = trackAR->getPhi();
@@ -1714,6 +1757,7 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) {
}
for (int ih=0;ih<nHitsOld;++ih) {
edm4hep::TrackerHit trkHit = hitVecOld[ih]->getTrackerHit();
+ //debug() << "old track's hit " << trkHit.id().collectionID << " " << trkHit.id().index << endmsg;
xh[ih+nHits] = trkHit.getPosition()[0];
yh[ih+nHits] = trkHit.getPosition()[1];
zh[ih+nHits] = float(trkHit.getPosition()[2]);
@@ -1810,6 +1854,15 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) {
// Attach hits belonging to the current track segment to
// the track already created
if (found == 1) {
+ if (iBad!=-1) {
+ if (iBad>=nHits)
+ debug() << "Bad hit " << iBad << " " << hitVecOld[iBad-nHits]->getTrackerHit().id().collectionID << "-" << hitVecOld[iBad-nHits]->getTrackerHit().id().index << endmsg;
+ else
+ debug() << "Bad hit " << iBad << " " << hitVec[iBad]->getTrackerHit().id().collectionID << "-" << hitVec[iBad]->getTrackerHit().id().index << endmsg;
+ }
+ else
+ debug() << "Bad hit " << iBad << endmsg;
+
trackOld->ClearTrackerHitExtendedVec();
for (int i=0;i<nHits;++i) {
TrackerHitExtended * trkHit = hitVec[i];
@@ -1824,6 +1877,7 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) {
for (int i=0;i<nHitsOld;++i) {
int icur = i+nHits;
TrackerHitExtended * trkHit = hitVecOld[i];
+ debug() << "add old hit " << i << " " << trkHit->getTrackerHit().id().collectionID << "-" << trkHit->getTrackerHit().id().index << endmsg;
trkHit->clearTrackVec();
if (icur == iBad) {
}
@@ -1848,7 +1902,8 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) {
trackOld->setChi2(chi2Min*float(ndf));
trackOld->setNDF(ndf);
trackOld->setCovMatrix(eparmin);
- // trackOld->setReferencePoint(refPointMin);
+ // TODO: by fucd, removed why in ILCSoft
+ //trackOld->setReferencePoint(refPointMin);
}
delete[] xh;
@@ -2027,7 +2082,8 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsFast() {
}
void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() {
-
+ debug() << "AttachRemainingVTXHitsSlow()" << endmsg;
+
TrackerHitExtendedVec nonAttachedHits;
nonAttachedHits.clear();
@@ -2050,11 +2106,9 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() {
if(isize>maxTrackSize)maxTrackSize = isize;
}
if (maxTrackSize<=3) {
- debug() << " Add non attached hit to list: id = " << hit->getTrackerHit().id() << endmsg;
+ debug() << " Add non attached hit to list: id = " << hit->getTrackerHit().id().index << endmsg;
nonAttachedHits.push_back( hit );
}
-
-
}
}
}
@@ -2065,7 +2119,7 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() {
int nTrk = int(_trackImplVec.size());
for (int iHit=0; iHit<nNotAttached; ++iHit) {
TrackerHitExtended * hit = nonAttachedHits[iHit];
- debug() << " Try hit: id = " << hit->getTrackerHit().id() << endmsg;
+ debug() << " Try hit: id = " << hit->getTrackerHit().id().index << endmsg;
int layer = getLayerID( hit->getTrackerHit() );
if (layer > _minimalLayerToAttach) {
float pos[3];
@@ -2097,7 +2151,7 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() {
distance = sqrt(distance);
if (distance<_minDistToDelta) {
consider = false;
- debug() << " hit: id = " << hit->getTrackerHit().id() << " condsidered delta together with hit " << trkhit2.id() << endmsg;
+ debug() << " hit: id = " << hit->getTrackerHit().id().index << " condsidered delta together with hit " << trkhit2.id().index << endmsg;
break;
}
}
@@ -2110,12 +2164,13 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() {
float z0 = trackAR->getZ0();
float omega = trackAR->getOmega();
float tanlambda = trackAR->getTanLambda();
+ //debug() << "attach hit now, field = " << _bField << " d0 = " << d0 << " phi0 = " << phi0 << " z0 = " << z0 << " omega = " << omega << " tanlambda = " << tanlambda << endmsg;
helix.Initialize_Canonical(phi0,d0,z0,omega,tanlambda,_bField);
float distance[3];
float time = helix.getDistanceToPoint(pos,distance);
if (time < 1.0e+10) {
if (distance[2] < minDist) {
- minDist = distance[2];
+ minDist = distance[2];
trackToAttach = trackAR;
}
}
@@ -2123,16 +2178,18 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() {
}
if (minDist < _minDistCutAttach && trackToAttach != NULL) {
int iopt = 2;
- debug() << " Hit: id = " << hit->getTrackerHit().id() << " : try attachement"<< endmsg;
+ debug() << " Hit: id = " << hit->getTrackerHit().id().index << " : try attachement"<< endmsg;
AttachHitToTrack(trackToAttach,hit,iopt);
} else {
- debug() << " Hit: id = " << hit->getTrackerHit().id() << " rejected due to distance cut of " <<_minDistCutAttach<< " min distance = " << minDist << endmsg;
+ debug() << " Hit: id = " << hit->getTrackerHit().id().index << " rejected due to distance cut of " <<_minDistCutAttach<< " min distance = " << minDist << endmsg;
}
}
}
}
void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() {
+ debug() << "AttachRemainingFTDHitsSlow()" << endmsg;
+
TrackerHitExtendedVec nonAttachedHits;
nonAttachedHits.clear();
@@ -2146,6 +2203,7 @@ void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() {
TrackerHitExtended * hit = hitVec[iH];
TrackExtendedVec& trackVec = hit->getTrackExtendedVec();
if (trackVec.size()==0)
+ debug() << " Add non attached hit to list: id = " << hit->getTrackerHit().id().index << endmsg;
nonAttachedHits.push_back( hit );
}
}
@@ -2173,7 +2231,7 @@ void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() {
// SJA:FIXME: check to see if allowing no hits in the same sensor vs no hits in the same layer works
// if (hit->getTrackerHit()->getType() == hitVector[IHIT]->getTrackerHit()->getType()) {
if (hit->getTrackerHit().getCellID() == hitVector[IHIT]->getTrackerHit().getCellID()) {
-
+ debug() << " hit: id = " << hit->getTrackerHit().id().index << " condsidered delta together with hit " << hitVector[IHIT]->getTrackerHit().id().index << endmsg;
consider = false;
break;
}
@@ -2193,7 +2251,7 @@ void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() {
float time = helix.getDistanceToPoint(pos,distance);
if (time < 1.0e+10) {
if (distance[2] < minDist) {
- minDist = distance[2];
+ minDist = distance[2];
trackToAttach = trackAR;
}
}
@@ -2202,8 +2260,12 @@ void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() {
}
if (minDist < _minDistCutAttach && trackToAttach != NULL) {
int iopt = 2;
+ debug() << " Hit: id = " << hit->getTrackerHit().id().index << " : try attachement"<< endmsg;
AttachHitToTrack(trackToAttach,hit,iopt);
- }
+ }
+ else {
+ debug() << " Hit: id = " << hit->getTrackerHit().id().index << " rejected due to distance cut of " <<_minDistCutAttach<< " min distance = " << minDist << endmsg;
+ }
}
}
@@ -2321,6 +2383,7 @@ void SiliconTrackingAlg::TrackingInFTD() {
if( iCodeOuter >= _sectorsFTD.size()){
error() << "iCodeOuter index out of range: iCodeOuter = " << iCodeOuter << " _sectorsFTD.size() = " << _sectorsFTD.size() << " exit(1) called from file " << __FILE__ << " line " << __LINE__<< endmsg;
+ info() << "iS=" << iS << " nLS=" << nLS[0] << " _nlayersFTD=" << _nlayersFTD << " ipOuter=" << ipOuter << endmsg;
exit(1);
}
@@ -2371,12 +2434,12 @@ void SiliconTrackingAlg::TrackingInFTD() {
// << "\nMiddle Hit Type "<< hitMiddle->getTrackerHit()->getType() << " z = " << hitMiddle->getTrackerHit().getPosition()[2]
// << "\nInner Hit Type "<< hitInner->getTrackerHit()->getType() << " z = " << hitInner->getTrackerHit().getPosition()[2] << endmsg;
- debug() << " "
- << std::setw(3) << ipOuter << " " << std::setw(3) << ipMiddle << " " << std::setw(3) << ipInner << " "
- << std::setw(3) << iS << " "
- << std::setw(3) << nLS[0] << " " << std::setw(3) << nLS[1] << " " << std::setw(3) << nLS[2] << " "
- << std::setw(3) << nOuter << " : " << std::setw(3) << nMiddle << " : " << std::setw(3) << nInner << " :: "
- << std::setw(3) << nOuter*nMiddle* nInner << endmsg;
+ debug() << "pattern "
+ << std::setw(3) << ipOuter << " " << std::setw(3) << ipMiddle << " " << std::setw(3) << ipInner << " "
+ << ((iS==0)?" +z ":" -z ")
+ << std::setw(3) << nLS[0] << " " << std::setw(3) << nLS[1] << " " << std::setw(3) << nLS[2] << " "
+ << std::setw(3) << nOuter << " : " << std::setw(3) << nMiddle << " : " << std::setw(3) << nInner << " :: "
+ << std::setw(3) << nOuter*nMiddle*nInner << endmsg;
TrackExtended * trackAR = TestTriplet(hitOuter,hitMiddle,hitInner,helix);
@@ -2583,7 +2646,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
TrackerHitExtendedVec& hitVec = trackAR->getTrackerHitExtendedVec();
int nHits = int(hitVec.size());
- //debug() << " Track " << iTrk << " has hits: " << nHits << endmsg;
+ debug() << " Track " << iTrk << " has hits: " << nHits << endmsg;
if( nHits >= _minimalHits) {
// int * lh = new int[nHits];
std::vector<int> lh;
@@ -2701,10 +2764,12 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
int nFit = 0;
for (int i=0; i<nHits; ++i) {
+ edm4hep::TrackerHit trkHit = hitVec[i]->getTrackerHit();
+ debug() << "TrackerHit " << i << " id = " << trkHit.id().collectionID << "-" << trkHit.id().index << endmsg;
// check if the hit has been rejected as being on the same layer and further from the helix lh==0
if (lh[i] == 1) {
- edm4hep::TrackerHit trkHit = hitVec[i]->getTrackerHit();
- debug() << "TrackerHit " << i << " id = " << trkHit.id() << endmsg;
+ //edm4hep::TrackerHit trkHit = hitVec[i]->getTrackerHit();
+ debug() << " accept" << endmsg;
nFit++;
if(trkHit.isAvailable()) {
trkHits.push_back(trkHit);
@@ -2723,6 +2788,10 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
debug() << "SiliconTrackingAlg::FinalRefit: Cannot fit less than 3 hits. Number of hits = " << trkHits.size() << endmsg;
continue ;
}
+
+ // test w/o fit
+ //continue;
+
//TrackImpl* Track = new TrackImpl ;
//auto track = trk_col->create();
//fucd
@@ -2786,6 +2855,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
*/
bool fit_backwards = IMarlinTrack::backward;
+ /*
MarlinTrk::IMarlinTrack* marlinTrk = nullptr;
try{
marlinTrk = _trksystem->createTrack();
@@ -2794,11 +2864,13 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
error() << "Cannot create MarlinTrack ! " << endmsg;
return;
}
+ */
int status = 0;
debug() << "call createFinalisedLCIOTrack now" << endmsg;
try {
- status = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, covMatrix, _bField, _maxChi2PerHit);
+ //status = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, covMatrix, _bField, _maxChi2PerHit);
+ status = m_fitTool->Fit(track, trkHits, covMatrix, _maxChi2PerHit, fit_backwards);
} catch (...) {
// delete Track;
// delete marlinTrk;
@@ -2819,29 +2891,49 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
std::vector<std::pair<edm4hep::TrackerHit , double> > hits_in_fit ;
std::vector<std::pair<edm4hep::TrackerHit , double> > outliers ;
std::vector<edm4hep::TrackerHit> all_hits;
- all_hits.reserve(300);
-
- marlinTrk->getHitsInFit(hits_in_fit);
-
+ //all_hits.reserve(300);
+
+ UTIL::BitField64 cellID_encoder( UTIL::ILDCellID0::encoder_string ) ;
+
+ //marlinTrk->getHitsInFit(hits_in_fit);
+ hits_in_fit = m_fitTool->GetHitsInFit();
+
for ( unsigned ihit = 0; ihit < hits_in_fit.size(); ++ihit) {
- debug() << "Hit id =" << hits_in_fit[ihit].first.id() << endmsg;
- edm4hep::TrackerHit trk = hits_in_fit[ihit].first;
- all_hits.push_back(trk);//hits_in_fit[ihit].first);
+ edm4hep::TrackerHit hit = hits_in_fit[ihit].first;
+ //all_hits.push_back(hit);//hits_in_fit[ihit].first);
+ if (m_tupleDebug) {
+ cellID_encoder.setValue(hit.getCellID());
+ m_layer[m_nhit] = cellID_encoder[lcio::ILDCellID0::layer];
+ m_chi2[m_nhit] = hits_in_fit[ihit].second;
+ m_used[m_nhit] = true;
+ debug() << "Hit id =" << hit.id().index << " layer = " << m_layer[m_nhit] << " det = " << cellID_encoder[lcio::ILDCellID0::subdet] << endmsg;
+ m_nhit++;
+ }
+ all_hits.push_back(hit);
}
- UTIL::BitField64 cellID_encoder( UTIL::ILDCellID0::encoder_string ) ;
-
MarlinTrk::addHitNumbersToTrack(&track, all_hits, true, cellID_encoder);
- marlinTrk->getOutliers(outliers);
-
+ //marlinTrk->getOutliers(outliers);
+ outliers = m_fitTool->GetOutliers();
+
for ( unsigned ihit = 0; ihit < outliers.size(); ++ihit) {
all_hits.push_back(outliers[ihit].first);
+ if (m_tupleDebug) {
+ cellID_encoder.setValue(outliers[ihit].first.getCellID());
+ m_layer[m_nhit] = cellID_encoder[lcio::ILDCellID0::layer];
+ m_chi2[m_nhit] = outliers[ihit].second;
+ m_used[m_nhit] = false;
+ m_nhit++;
+ }
}
MarlinTrk::addHitNumbersToTrack(&track, all_hits, false, cellID_encoder);
-
- delete marlinTrk;
+ if (m_tupleDebug) m_tupleDebug->write();
+
+ //delete marlinTrk;
+ m_fitTool->Clear();
+
#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
int nhits_in_vxd = track.getSubdetectorHitNumbers(0);
int nhits_in_ftd = track.getSubdetectorHitNumbers(1);
@@ -2852,8 +2944,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
int nhits_in_sit = track.getSubDetectorHitNumbers(2);
#endif
- //debug() << " Hit numbers for Track "<< track.id() << ": "
- debug() << " Hit numbers for Track "<< iTrk <<": "
+ debug() << " Hit numbers for Track "<< iTrk <<" (id=" << track.id().index << "): "
<< " vxd hits = " << nhits_in_vxd
<< " ftd hits = " << nhits_in_ftd
<< " sit hits = " << nhits_in_sit
@@ -2865,13 +2956,13 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
if( status != IMarlinTrack::success ) {
//delete track;
- debug() << "FinalRefit: Track fit failed with error code " << status << " track dropped. Number of hits = "<< trkHits.size() << endmsg;
+ debug() << "FinalRefit: Track fit failed with error code " << status << " track dropped. Number of hits = "<< trkHits.size() << endmsg;
continue ;
}
if( track.getNdf() < 0) {
//delete track;
- debug() << "FinalRefit: Track fit returns " << track.getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg;
+ debug() << "FinalRefit: Track fit returns " << track.getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg;
//delete track;
continue ;
}
@@ -2885,7 +2976,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
if(trkStateIP.location !=1) continue;
/*
if (trkStateIP == 0) {
- debug() << "SiliconTrackingAlg::FinalRefit: Track fit returns " << track->getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg;
+ debug() << "SiliconTrackingAlg::FinalRefit: Track fit returns " << track->getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg;
throw EVENT::Exception( std::string("SiliconTrackingAlg::FinalRefit: trkStateIP pointer == NULL ") ) ;
}
*/
@@ -2898,8 +2989,8 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
float cov[15];
- for (int i = 0 ; i<15 ; ++i) {
- cov[i] = trkStateIP.covMatrix.operator[](i);
+ for (int icov = 0 ; icov < 15 ; ++icov) {
+ cov[icov] = trkStateIP.covMatrix.operator[](icov);
}
trackAR->setCovMatrix(cov);
@@ -2925,8 +3016,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
}
if(m_dumpTime&&m_tuple) m_timeKalman = stopwatch.RealTime()*1000;
- debug() << " -> run " << _nRun << " event " << _nEvt << endmsg;
- debug() << "Number of Si tracks = " << nSiSegments << endmsg;
+ info() << " -> run " << _nRun << " event " << _nEvt << " Number of Si tracks = " << nSiSegments << endmsg;
debug() << "Total 4-momentum of Si Tracks : E = " << std::setprecision(7) << eTot
<< " Px = " << pxTot
<< " Py = " << pyTot
@@ -3076,6 +3166,9 @@ StatusCode SiliconTrackingAlg::setupGearGeom(){
}
}
+
+ info() << "nvxd = " << _nLayersVTX << " nsit = " << _nLayersSIT << " nftd = " << _nlayersFTD << endmsg;
+
return StatusCode::SUCCESS;
}
diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h
index e8e915b5..c457d446 100644
--- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h
+++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h
@@ -23,6 +23,7 @@
#include "DataHelper/HelixClass.h"
#include "TrackSystemSvc/IMarlinTrack.h"
+#include "Tracking/ITrackFitterTool.h"
#include <UTIL/BitField64.h>
#include <UTIL/ILDConf.h>
@@ -211,6 +212,7 @@ class SiliconTrackingAlg : public GaudiAlgorithm {
/** pointer to the IMarlinTrkSystem instance
*/
MarlinTrk::IMarlinTrkSystem* _trksystem ;
+ ToolHandle<ITrackFitterTool> m_fitTool;
//bool _runMarlinTrkDiagnostics;
//std::string _MarlinTrkDiagnosticsName;
typedef std::vector<int> IntVec;
@@ -238,7 +240,7 @@ class SiliconTrackingAlg : public GaudiAlgorithm {
Gaudi::Property<float> _cutOnPt{this, "CutOnPt", 0.05};
Gaudi::Property<int> _minimalHits{this, "MinimalHits",3};
Gaudi::Property<int> _nHitsChi2{this, "NHitsChi2", 5};
- Gaudi::Property<int> _max_hits_per_sector{this, "MaxHitsPerSector", 100};
+ Gaudi::Property<int> _max_hits_per_sector{this, "MaxHitsPerSector", 200};
Gaudi::Property<int> _attachFast{this, "FastAttachment", 0};
Gaudi::Property<bool> _useSIT{this, "UseSIT", true};
Gaudi::Property<float> _initialTrackError_d0{this, "InitialTrackErrorD0",1e6};
@@ -253,8 +255,11 @@ class SiliconTrackingAlg : public GaudiAlgorithm {
Gaudi::Property<bool> _ElossOn{this, "EnergyLossOn", true};
Gaudi::Property<bool> _SmoothOn{this, "SmoothOn", true};
Gaudi::Property<float> _helix_max_r{this, "HelixMaxR", 2000.};
- Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true};
+ Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true};
+ Gaudi::Property<bool> m_debug{this, "Debug", false};
+
//std::vector<int> _colours;
+ Gaudi::Property<std::string> m_fitToolName{this, "FitterTool", "KalTestTool/KalTest111"};
/** helper function to get collection using try catch block */
//LCCollection* GetCollection( LCEvent * evt, std::string colName ) ;
@@ -302,10 +307,15 @@ class SiliconTrackingAlg : public GaudiAlgorithm {
std::vector<TrackerHitExtendedVec> _sectors;
std::vector<TrackerHitExtendedVec> _sectorsFTD;
- NTuple::Tuple* m_tuple;
+ NTuple::Tuple* m_tuple = nullptr;
NTuple::Item<float> m_timeTotal;
NTuple::Item<float> m_timeKalman;
+ NTuple::Tuple* m_tupleDebug = nullptr;
+ NTuple::Item<int> m_nhit;
+ NTuple::Array<int> m_layer;
+ NTuple::Array<float> m_chi2;
+ NTuple::Array<bool> m_used;
/**
* A helper class to allow good code readability by accessing tracks with N hits.
* As the smalest valid track contains three hits, but the first index in a vector is 0,
diff --git a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp
index d4741f2e..c8b6d1ea 100644
--- a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp
+++ b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp
@@ -41,12 +41,16 @@ StatusCode TrackSubsetAlg::initialize() {
_nRun = 0 ;
_nEvt = 0 ;
+ // Set up the track fit tool
+ m_fitTool = ToolHandle<ITrackFitterTool>(m_fitToolName.value());
+
if(m_dumpTime){
NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name());
if ( !nt1 ) {
m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time");
if ( 0 != m_tuple ) {
m_tuple->addItem ("timeTotal", m_timeTotal ).ignore();
+ m_tuple->addItem ("timeKalman", m_timeKalman ).ignore();
}
else {
fatal() << "Cannot book MyTuples/Time"+name() <<endmsg;
@@ -220,7 +224,7 @@ StatusCode TrackSubsetAlg::execute(){
debug() << "Fitting and saving of the tracks" << endmsg;
//auto trkCol = _outColHdl.createAndPut();
-
+ auto stopwatch_kalman = TStopwatch();
for( unsigned i=0; i < accepted.size(); i++ ){
edm4hep::MutableTrack trackImpl;
@@ -266,14 +270,15 @@ StatusCode TrackSubsetAlg::execute(){
bool fit_backwards = MarlinTrk::IMarlinTrack::backward;
- MarlinTrk::IMarlinTrack* marlinTrk = _trkSystem->createTrack();
+ //MarlinTrk::IMarlinTrack* marlinTrk = _trkSystem->createTrack();
int error = 0;
try {
- error = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trackerHits, &trackImpl, fit_backwards, covMatrix, _bField, _maxChi2PerHit);
-
+ //error = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trackerHits, &trackImpl, fit_backwards, covMatrix, _bField, _maxChi2PerHit);
+ error = m_fitTool->Fit(trackImpl, trackerHits, covMatrix, _maxChi2PerHit, fit_backwards);
+
} catch (...) {
// delete Track;
@@ -288,10 +293,11 @@ StatusCode TrackSubsetAlg::execute(){
std::vector<std::pair<edm4hep::TrackerHit , double> > hits_in_fit ;
std::vector<std::pair<edm4hep::TrackerHit , double> > outliers ;
std::vector<edm4hep::TrackerHit> all_hits;
- all_hits.reserve(300);
-
- marlinTrk->getHitsInFit(hits_in_fit);
+ //all_hits.reserve(300);
+ //marlinTrk->getHitsInFit(hits_in_fit);
+ hits_in_fit = m_fitTool->GetHitsInFit();
+
for ( unsigned ihit = 0; ihit < hits_in_fit.size(); ++ihit) {
all_hits.push_back(hits_in_fit[ihit].first);
}
@@ -300,15 +306,17 @@ StatusCode TrackSubsetAlg::execute(){
MarlinTrk::addHitNumbersToTrack(&trackImpl, all_hits, true, cellID_encoder);
- marlinTrk->getOutliers(outliers);
-
+ //marlinTrk->getOutliers(outliers);
+ outliers = m_fitTool->GetOutliers();
+
for ( unsigned ihit = 0; ihit < outliers.size(); ++ihit) {
all_hits.push_back(outliers[ihit].first);
}
MarlinTrk::addHitNumbersToTrack(&trackImpl, all_hits, false, cellID_encoder);
- delete marlinTrk;
+ //delete marlinTrk;
+ m_fitTool->Clear();
if( error != MarlinTrk::IMarlinTrack::success ) {
//delete trackImpl;
@@ -355,6 +363,7 @@ StatusCode TrackSubsetAlg::execute(){
if(m_dumpTime&&m_tuple){
m_timeTotal = stopwatch.RealTime()*1000;
+ m_timeKalman = stopwatch_kalman.RealTime()*1000;
m_tuple->write();
}
diff --git a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h
index 1023eba9..534524f9 100644
--- a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h
+++ b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h
@@ -8,6 +8,7 @@
#include "edm4hep/TrackerHitCollection.h"
//#include "edm4hep/Track.h"
#include "TrackSystemSvc/IMarlinTrkSystem.h"
+#include "Tracking/ITrackFitterTool.h"
#include "Math/ProbFunc.h"
@@ -56,6 +57,7 @@ class TrackSubsetAlg : public GaudiAlgorithm {
protected:
MarlinTrk::IMarlinTrkSystem* _trkSystem;
+ ToolHandle<ITrackFitterTool> m_fitTool;
/* Input collection */
std::vector<DataHandle<edm4hep::TrackCollection>* > _inTrackColHdls;
std::vector<DataHandle<edm4hep::TrackerHitCollection>* > _inTrackerHitColHdls;
@@ -76,6 +78,7 @@ class TrackSubsetAlg : public GaudiAlgorithm {
Gaudi::Property<double> _maxChi2PerHit{this, "MaxChi2PerHit", 1e2};
Gaudi::Property<double> _omega{this, "Omega", 0.75};
Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true};
+ Gaudi::Property<std::string> m_fitToolName{this, "FitterTool", "KalTestTool/KalTest111"};
float _bField;
@@ -84,6 +87,7 @@ class TrackSubsetAlg : public GaudiAlgorithm {
NTuple::Tuple* m_tuple;
NTuple::Item<float> m_timeTotal;
+ NTuple::Item<float> m_timeKalman;
};
/** A functor to return whether two tracks are compatible: The criterion is if the share a TrackerHit or more */
diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
index d218e041..937adba6 100644
--- a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
+++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
@@ -201,6 +201,9 @@ StatusCode ClupatraAlg::initialize() {
// don't need, since Gaudi Algorithm will print all Property
//printParameters() ;
+ // Set up the track fit tool
+ m_fitTool = ToolHandle<ITrackFitterTool>(m_fitToolName.value());
+
auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc");
if ( !_trackSystemSvc ) {
error() << "Fail to find TrackSystemSvc ..." << endmsg;
diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
index 962b1c54..fb451483 100644
--- a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
+++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
@@ -22,6 +22,7 @@
#include "TrackSystemSvc/HelixTrack.h"
#include "TrackSystemSvc/HelixFit.h"
#include "TrackSystemSvc/IMarlinTrack.h"
+#include "Tracking/ITrackFitterTool.h"
namespace MarlinTrk{
class IMarlinTrkSystem ;
@@ -139,6 +140,7 @@ class ClupatraAlg : public GaudiAlgorithm {
Gaudi::Property<bool> _MSOn{this, "MultipleScatteringOn", false};
Gaudi::Property<bool> _ElossOn{this, "EnergyLossOn", true};
Gaudi::Property<bool> _SmoothOn{this, "SmoothOn", false};
+ Gaudi::Property<std::string> m_fitToolName{this, "FitterTool", "KalTestTool/KalTest010"};
DataHandle<edm4hep::TrackerHitCollection> _TPCHitCollectionHandle{"TPCTrackerHits", Gaudi::DataHandle::Reader, this};
@@ -154,6 +156,7 @@ class ClupatraAlg : public GaudiAlgorithm {
int _nEvt ;
MarlinTrk::IMarlinTrkSystem* _trksystem ;
+ ToolHandle<ITrackFitterTool> m_fitTool;
const gear::TPCParameters* _gearTPC ;
diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
old mode 100644
new mode 100755
index ee5f1feb..95265a4c
--- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
+++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
@@ -6,7 +6,7 @@
#include <GearSvc/IGearSvc.h>
#include <edm4hep/TrackerHit.h>
-#include <edm4hep/TrackerHit.h>
+#include <edm4hep/TrackerHitPlane.h>
#include <edm4hep/Track.h>
#if __has_include("edm4hep/EDM4hepVersion.h")
#include "edm4hep/EDM4hepVersion.h"
@@ -94,7 +94,11 @@ std::string toString( int iTrk, edm4hep::Track tpcTrack, float bField=3.5 ) {
float pzTPC = helixTPC.getMomentum()[2];
const float ptot = sqrt(pxTPC*pxTPC+pyTPC*pyTPC+pzTPC*pzTPC);
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ sprintf(strg,"%3i %5i %9.3f %9.3f %9.3f %7.2f %7.2f %7.2f %4i %4i %8.3f %8i",iTrk,tpcTrack.id().index,
+#else
sprintf(strg,"%3i %5i %9.3f %9.3f %9.3f %7.2f %7.2f %7.2f %4i %4i %8.3f %8i",iTrk,tpcTrack.id(),
+#endif
ptot, d0TPC,z0TPC,pxTPC,pyTPC,pzTPC,nHits,ndfTPC,Chi2TPC,nlinkedTracks);
return std::string( strg ) ;
@@ -154,6 +158,9 @@ StatusCode FullLDCTrackingAlg::initialize() {
PIOVER2 = 0.5*PI;
TWOPI = 2*PI;
+ // Set up the track fit tool
+ m_fitTool = ToolHandle<ITrackFitterTool>(m_fitToolName.value());
+
if(m_dumpTime){
NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name());
if ( !nt1 ) {
@@ -414,13 +421,14 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
debug() << "trkHits ready" << endmsg;
bool fit_backwards = IMarlinTrack::backward;
- MarlinTrk::IMarlinTrack* marlinTrk = _trksystem->createTrack();
- debug() << "marlinTrk created " << endmsg;
+ //MarlinTrk::IMarlinTrack* marlinTrk = _trksystem->createTrack();
+ //debug() << "marlinTrk created " << endmsg;
int error_code = 0;
try {
- error_code = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, &ts_initial, _bField, _maxChi2PerHit);
+ //error_code = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, &ts_initial, _bField, _maxChi2PerHit);
+ error_code = m_fitTool->Fit(track, trkHits, ts_initial, _maxChi2PerHit, fit_backwards);
} catch (...) {
// delete track;
@@ -445,8 +453,9 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
std::vector<edm4hep::TrackerHit> all_hits;
all_hits.reserve(hits_in_fit.size());
- marlinTrk->getHitsInFit(hits_in_fit);
-
+ //marlinTrk->getHitsInFit(hits_in_fit);
+ hits_in_fit = m_fitTool->GetHitsInFit();
+
for ( unsigned ihit = 0; ihit < hits_in_fit.size(); ++ihit) {
all_hits.push_back(hits_in_fit[ihit].first);
}
@@ -455,17 +464,18 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
MarlinTrk::addHitNumbersToTrack(&track, all_hits, true, cellID_encoder);
- marlinTrk->getOutliers(outliers);
-
+ //marlinTrk->getOutliers(outliers);
+ outliers = m_fitTool->GetOutliers();
+
for ( unsigned ihit = 0; ihit < outliers.size(); ++ihit) {
all_hits.push_back(outliers[ihit].first);
}
MarlinTrk::addHitNumbersToTrack(&track, all_hits, false, cellID_encoder);
-
- delete marlinTrk;
-
+ //delete marlinTrk;
+ m_fitTool->Clear();
+
if( error_code != IMarlinTrack::success ) {
debug() << "FullLDCTrackingAlg::AddTrackColToEvt: Track fit failed with error code " << error_code << " track dropped. Number of hits = "<< trkHits.size() << endmsg;
//delete Track;
@@ -534,6 +544,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
int nhits_in_tpc = track.getSubDetectorHitNumbers(3);
int nhits_in_set = track.getSubDetectorHitNumbers(4);
#endif
+
//int nhits_in_vxd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::VXD - 2 ];
//int nhits_in_ftd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - 2 ];
//int nhits_in_sit = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::SIT - 2 ];
@@ -1264,7 +1275,11 @@ void FullLDCTrackingAlg::prepareVectors() {
const float pTot = sqrt(pxSi*pxSi+pySi*pySi+pzSi*pzSi);
const int ndfSi = trackExt->getNDF();
float Chi2Si = trackExt->getChi2()/float(trackExt->getNDF());
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ sprintf(strg,"%3i %5i %9.3f %9.3f %9.3f %7.2f %7.2f %7.2f %4i %4i %8.3f",iTrk, siTrack.id().index, pTot, d0Si,z0Si,pxSi,pySi,pzSi,nHits, ndfSi, Chi2Si);
+#else
sprintf(strg,"%3i %5i %9.3f %9.3f %9.3f %7.2f %7.2f %7.2f %4i %4i %8.3f",iTrk, siTrack.id(), pTot, d0Si,z0Si,pxSi,pySi,pzSi,nHits, ndfSi, Chi2Si);
+#endif
debug() << strg << endmsg;
if(nHits>0){
diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
index cfc80487..2e9b9586 100755
--- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
+++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
@@ -19,6 +19,7 @@
#include "TrackSystemSvc/MarlinTrkUtils.h"
#include "TrackSystemSvc/IMarlinTrack.h"
+#include "Tracking/ITrackFitterTool.h"
#include "edm4hep/SimTrackerHitCollection.h"
#include "edm4hep/TrackerHitCollection.h"
@@ -319,6 +320,7 @@ protected:
/** pointer to the IMarlinTrkSystem instance
*/
MarlinTrk::IMarlinTrkSystem* _trksystem ;
+ ToolHandle<ITrackFitterTool> m_fitTool;
Gaudi::Property<bool> _MSOn{this, "MultipleScatteringOn", true};
Gaudi::Property<bool> _ElossOn{this, "EnergyLossOn", true};
@@ -397,6 +399,7 @@ protected:
Gaudi::Property<int> _maxAllowedSiHitRejectionsForTrackCombination{this, "MaxAllowedSiHitRejectionsForTrackCombination", 2};
Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true};
//float _dPCutForForcedMerging;
+ Gaudi::Property<std::string> m_fitToolName{this, "FitterTool", "KalTestTool/KalTest111"};
bool _runMarlinTrkDiagnostics = false;
std::string _MarlinTrkDiagnosticsName;
diff --git a/Service/GearSvc/src/GearSvc.h b/Service/GearSvc/src/GearSvc.h
index 3048526c..312491eb 100644
--- a/Service/GearSvc/src/GearSvc.h
+++ b/Service/GearSvc/src/GearSvc.h
@@ -10,7 +10,7 @@ class GearSvc : public extends<Service, IGearSvc>
{
public:
GearSvc(const std::string& name, ISvcLocator* svc);
- ~GearSvc();
+ virtual ~GearSvc();
gear::GearMgr* getGearMgr() override;
@@ -28,6 +28,7 @@ class GearSvc : public extends<Service, IGearSvc>
TGeoNode* FindNode(TGeoNode* mother, char* name);
Gaudi::Property<std::string> m_gearFile{this, "GearXMLFile", ""};
+ Gaudi::Property<float> m_field{this, "MagneticField", 0};
gear::GearMgr* m_gearMgr;
};
--
GitLab