From 98613701ea819cdb5786b9168833f3c8c5d91f39 Mon Sep 17 00:00:00 2001
From: FU Chengdong <fucd@ihep.ac.cn>
Date: Wed, 14 Aug 2024 03:36:59 +0000
Subject: [PATCH] =?UTF-8?q?Tracking=EF=BC=9Alittle=20update=20pack=20for?=
=?UTF-8?q?=20geometry=20and=20tracking?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
.../SiTrackerStaggeredLadder_v01_geo.cpp | 2 +-
Digitisers/SimpleDigi/src/TPCDigiAlg.cpp | 21 ++++++++-------
.../Tracking/src/Clupatra/ClupatraAlg.cpp | 26 +++++++++++++++++--
.../Tracking/src/Clupatra/ClupatraAlg.h | 7 ++++-
Service/TrackSystemSvc/src/MarlinTrkUtils.cc | 17 ++++++++----
5 files changed, 55 insertions(+), 18 deletions(-)
diff --git a/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp
index b8245d66..64fe00bb 100644
--- a/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp
+++ b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp
@@ -101,7 +101,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h
std::string deadwireVis = x_display.attr<string>(_Unicode(deadwire));
//fetch the shell parameters
- if (x_det.hasAttr(_Unicode(shell))) {
+ if (x_det.hasChild(_Unicode(shell))) {
xml_comp_t x_shell(x_det.child(_Unicode(shell)));
double rmin_shell = x_shell.rmin();
double rmax_shell = x_shell.rmax();
diff --git a/Digitisers/SimpleDigi/src/TPCDigiAlg.cpp b/Digitisers/SimpleDigi/src/TPCDigiAlg.cpp
index 11931a84..6b0d042a 100644
--- a/Digitisers/SimpleDigi/src/TPCDigiAlg.cpp
+++ b/Digitisers/SimpleDigi/src/TPCDigiAlg.cpp
@@ -850,7 +850,8 @@ StatusCode TPCDigiAlg::execute()
// create a tpc voxel hit and store it for this row
Voxel_tpc * atpcVoxel = new Voxel_tpc(iRowHit,iPhiHit,iZHit, thisPoint, edep, tpcRPhiRes, tpcZRes);
-
+ debug() << "to seed hit position: " << atpcVoxel->getX() << "," << atpcVoxel->getY() << "," << atpcVoxel->getZ()
+ << " iRow=" << iRowHit << " iPhi=" << iPhiHit << " padIndex=" << padIndex << endmsg;
_tpcRowHits.at(iRowHit).push_back(atpcVoxel);
++numberOfVoxelsCreated;
@@ -1007,7 +1008,8 @@ StatusCode TPCDigiAlg::execute()
if( momentum_set ){
const edm4hep::Vector3f Momentum1 = Hit1.getMomentum() ;
- const edm4hep::Vector3f Momentum2 = Hit1.getMomentum() ;
+ //fucd: Hit1.getMomentum() -> Hit2.getMomentum()
+ const edm4hep::Vector3f Momentum2 = Hit2.getMomentum() ;
CLHEP::Hep3Vector mom1(Momentum1[0],Momentum1[1],Momentum1[2]);
CLHEP::Hep3Vector mom2(Momentum2[0],Momentum2[1],Momentum2[2]);
@@ -1185,7 +1187,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
Voxel_tpc* seed_hit = aVoxel;
// if( seed_hit->getRowIndex() > 5 ) return ;
- debug() << "==============" << endmsg;
+
//store hit variables
edm4hep::MutableTrackerHit trkHit;// = _trkhitVec->create();
//now the hit pos has to be smeared
@@ -1205,7 +1207,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
// make sure the hit is not smeared beyond the TPC Max DriftLength
if( fabs(point.z()) > gearTPC.getMaxDriftLength() ) point.setZ( (fabs(point.z()) / point.z() ) * gearTPC.getMaxDriftLength() );
- debug() << "==============" << endmsg;
+ debug() << seed_hit->getX() << "," << seed_hit->getY() << "," << seed_hit->getZ() << " -> " << point << endmsg;
edm4hep::Vector3d pos(point.x(),point.y(),point.z());
trkHit.setPosition(pos);
trkHit.setEDep(seed_hit->getEDep());
@@ -1256,7 +1258,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
<< "\n" ;
throw errorMsg.str();
}
- debug() << "==============" << endmsg;
+
// For no error in R
std::array<float,TRKHITNCOVMATRIX> covMat={sin(unsmearedPhi)*sin(unsmearedPhi)*tpcRPhiRes*tpcRPhiRes,
-cos(unsmearedPhi)*sin(unsmearedPhi)*tpcRPhiRes*tpcRPhiRes,
@@ -1274,7 +1276,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
<< "\n" ;
throw errorMsg.str();
}
- debug() << "==============" << endmsg;
+
if(pos[0]*pos[0]+pos[1]*pos[1]>0.0){
// push back the SimTHit for this TrackerHit
@@ -1291,7 +1293,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
_NRecTPCHits++;
}
- debug() << "==============" << endmsg;
+
#ifdef DIGIPLOTS
// edm4hep::SimTrackerHit* theSimHit = _tpcHitMap[seed_hit];
// double rSimSqrd = theSimHit->getPosition()[0]*theSimHit->getPosition()[0] + theSimHit->getPosition()[1]*theSimHit->getPosition()[1];
@@ -1333,7 +1335,7 @@ void TPCDigiAlg::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerge){
unsigned number_of_hits_to_merge = hitsToMerge->size();
-
+ debug() << "number_of_hits_to_merge = " << number_of_hits_to_merge << endmsg;
for(unsigned int ihitCluster = 0; ihitCluster < number_of_hits_to_merge; ++ihitCluster){
sumZ += hitsToMerge->at(ihitCluster)->getZ();
@@ -1345,6 +1347,7 @@ void TPCDigiAlg::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerge){
if (_use_raw_hits_to_store_simhit_pointer) {
trkHit.addToRawHits(_tpcHitMap[hitsToMerge->at(ihitCluster)].getObjectID());
}
+ debug() << "raw hit: " << _tpcHitMap[hitsToMerge->at(ihitCluster)].getPosition() << endmsg;
auto rel = _relCol->create();
rel.setRec (trkHit);
@@ -1384,7 +1387,7 @@ void TPCDigiAlg::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerge){
if( fabs(point.z()) > gearTPC.getMaxDriftLength() ) point.setZ( (fabs(point.z()) / point.z() ) * gearTPC.getMaxDriftLength() );
double pos[3] = {point.x(),point.y(),point.z()};
-
+ debug() << "to hit: " << pos[0] << "," << pos[1] << "," << pos[2] << endmsg;
//---------------------------------------------------------------------------------
trkHit.setPosition(pos);
trkHit.setEDep(sumEDep);
diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
index 8889bc19..678beddd 100644
--- a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
+++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
@@ -46,7 +46,7 @@
#include "edm4hep/TrackerHitCollection.h"
#include "edm4hep/TrackCollection.h"
// #include "edm4hep/TrackerHitPlane.h"
-
+#include <TStopwatch.h>
using namespace edm4hep ;
using namespace MarlinTrk ;
@@ -227,6 +227,23 @@ StatusCode ClupatraAlg::initialize() {
//tree->Branch("omega", &omega, "omega/D");
//tree->Branch("totalCandidates", &totalCandidates, "totalCandidates/I");
//tree->Branch("eventNumber", &_nEvt, "eventNumber/I");
+ if(_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;
+ return StatusCode::FAILURE;
+ }
+ }
+ else{
+ m_tuple = nt1;
+ }
+ }
return GaudiAlgorithm::initialize();
@@ -237,7 +254,7 @@ StatusCode ClupatraAlg::execute() {
info() << "Clupatra Algorithm started" << endmsg;
-
+ auto stopwatch = TStopwatch();
// clock_t start = clock() ;
Timer timer ;
unsigned t_init = timer.registerTimer(" initialization " ) ;
@@ -1264,6 +1281,11 @@ StatusCode ClupatraAlg::execute() {
delete t;
}
+ if(m_tuple){
+ m_timeTotal = stopwatch.RealTime()*1000;
+ m_tuple->write();
+ }
+
_nEvt++ ;
TrackInfo_of_edm4hepTrack.clear();
diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
index fb451483..2e81b9e6 100644
--- a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
+++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h
@@ -24,6 +24,8 @@
#include "TrackSystemSvc/IMarlinTrack.h"
#include "Tracking/ITrackFitterTool.h"
+#include "GaudiKernel/NTuple.h"
+
namespace MarlinTrk{
class IMarlinTrkSystem ;
}
@@ -141,7 +143,7 @@ class ClupatraAlg : public GaudiAlgorithm {
Gaudi::Property<bool> _ElossOn{this, "EnergyLossOn", true};
Gaudi::Property<bool> _SmoothOn{this, "SmoothOn", false};
Gaudi::Property<std::string> m_fitToolName{this, "FitterTool", "KalTestTool/KalTest010"};
-
+ Gaudi::Property<bool> _DumpTime{this, "DumpTime", false};
DataHandle<edm4hep::TrackerHitCollection> _TPCHitCollectionHandle{"TPCTrackerHits", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::TrackerHitCollection> _SITHitCollectionHandle{"SIDTrackerHits", Gaudi::DataHandle::Reader, this};
@@ -165,6 +167,9 @@ class ClupatraAlg : public GaudiAlgorithm {
double pt;
int totalCandidates;
+ NTuple::Tuple* m_tuple = nullptr;
+ NTuple::Item<float> m_timeTotal;
+ NTuple::Item<float> m_timeKalman;
// NNClusterer* _clusterer ;
} ;
diff --git a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
index bf2ab05d..2a4c9d8a 100644
--- a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
+++ b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
@@ -212,7 +212,16 @@ namespace MarlinTrk {
///////////////////////////////////////////////////////
// add hits to IMarlinTrk
///////////////////////////////////////////////////////
-
+ std::vector<edm4hep::TrackerHit> hit_list_copy;
+ if (fit_backwards) {
+ std::reverse_copy(hit_list.begin(), hit_list.end(), std::back_inserter(hit_list_copy));
+ }
+ else {
+ hit_list_copy.reserve(hit_list.size());
+ hit_list_copy.insert(hit_list_copy.end(), hit_list.begin(), hit_list.end());
+ }
+
+ hit_list_copy.swap(hit_list);
std::vector<edm4hep::TrackerHit>::iterator it = hit_list.begin();
// start by trying to add the hits to the track we want to finally use.
@@ -262,9 +271,8 @@ namespace MarlinTrk {
//std::cout << "ERROR<<<<<MarlinTrk::createFit : Cannot fit less with less than " << MIN_NDF << " degrees of freedom. Number of hits = " << added_hits.size() << " ndof = " << ndof_added << std::endl;
return IMarlinTrack::bad_intputs;
}
-
-
-
+
+ hit_list_copy.swap(hit_list);
///////////////////////////////////////////////////////
// set the initial track parameters
///////////////////////////////////////////////////////
@@ -347,7 +355,6 @@ namespace MarlinTrk {
helixTrack.moveRefPoint(hit_list.front().getPosition()[0], hit_list.front().getPosition()[1], hit_list.front().getPosition()[2]);
}
-
const float referencePoint[3] = { helixTrack.getRefPointX() , helixTrack.getRefPointY() , helixTrack.getRefPointZ() };
pre_fit->D0 = helixTrack.getD0();
--
GitLab