diff --git a/Detector/DetDriftChamber/compact/det.xml b/Detector/DetDriftChamber/compact/det.xml
index 3c86627d3bde171386776eb74904f81129f0fe51..fcc55f72795c20cc355c0fb277cd50b6a09e2e38 100644
--- a/Detector/DetDriftChamber/compact/det.xml
+++ b/Detector/DetDriftChamber/compact/det.xml
@@ -32,7 +32,7 @@
<constant name="DC_rbegin" value="800*mm"/>
<constant name="DC_rend" value="1800*mm"/>
- <constant name="DC_layer_number" value="100"/>
+ <constant name="DC_layer_number" value="55"/>
<constant name="Alpha" value="0*deg"/>
<constant name="Gas_radius_min" value="DC_rbegin+DC_inner_wall_thickness+DC_safe_distance"/>
@@ -44,7 +44,7 @@
<constant name="Gas_half_length" value="DC_half_length-DC_Endcap_dz-DC_safe_distance"/>
<constant name="Gas_length" value="Gas_half_length*2"/>
- <constant name="DC_cell_width" value="10*mm"/>
+ <constant name="DC_cell_width" value="18*mm"/>
<constant name="DC_inner_wall_radius_min" value="DC_rbegin"/>
<constant name="DC_inner_wall_radius_max" value="DC_rbegin+DC_inner_wall_thickness"/>
diff --git a/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc
index d97a82d6e86aa8b2d97b52b6504011d571c6f172..6a1e0cc20da795c05f1e7bc328060f5b232f91a8 100644
--- a/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc
+++ b/Reconstruction/PFA/Arbor/src/ArborToolLCIO.cc
@@ -35,6 +35,8 @@
#include <DDRec/CellIDPositionConverter.h>
#include "DetInterface/IGeomSvc.h"
+#include "podio/podioVersion.h"
+
using namespace std;
/*
void ClusterBuilding( LCEvent * evtPP, std::string Name, std::vector<CalorimeterHit*> Hits, std::vector< std::vector<int> > BranchOrder, int DHCALFlag )
@@ -859,7 +861,12 @@ edm4hep::ClusterCollection* ArborToolLCIO::ClusterVecMerge( std::vector<edm4hep:
edm4hep::Cluster Mergebranch_A;
edm4hep::Cluster Mergebranch_B;
edm4hep::Cluster tmpMergebranch;
+#if PODIO_BUILD_VERSION < PODIO_VERSION(0, 17, 4)
edm4hep::Cluster Mainbranch (0);
+#else
+ auto Mainbranch = edm4hep::Cluster::makeEmpty();
+#endif
+
TVector3 tmpClusterSeedPos, MBSeedPos;
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/MCParticleCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/MCParticleCreator.cpp
index 1912eee92592f0b8071fcabed38259a61c049f87..48908a6439a27cc352020a83572640803034f777 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/MCParticleCreator.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/MCParticleCreator.cpp
@@ -53,7 +53,7 @@ pandora::StatusCode MCParticleCreator::CreateMCParticles(const CollectionMaps& c
mcParticleParameters.m_particleId = pMcParticle.getPDG();
mcParticleParameters.m_mcParticleType = pandora::MC_3D;
mcParticleParameters.m_pParentAddress = &pMcParticle;
- unsigned int p_id = pMcParticle.id();
+ unsigned int p_id = pMcParticle.id().index;
//auto p_mc = const_cast<edm4hep::MCParticle*>(&pMcParticle);
auto p_mc = &pMcParticle;
(*m_id_pMC_map) [p_id] = p_mc;
@@ -131,8 +131,8 @@ pandora::StatusCode MCParticleCreator::CreateCaloHitToMCParticleRelationships(co
auto conb = pSimHit.getContributions(iCont);
auto ipa = conb.getParticle();
float ien = conb.getEnergy();
- if( m_id_pMC_map->find(ipa.id()) == m_id_pMC_map->end() ) continue;
- auto p_tmp = (*m_id_pMC_map)[ipa.id()];
+ if( m_id_pMC_map->find(ipa.id().index) == m_id_pMC_map->end() ) continue;
+ auto p_tmp = (*m_id_pMC_map)[ipa.id().index];
mcParticleToEnergyWeightMap[p_tmp] += ien;
}
@@ -190,13 +190,13 @@ pandora::StatusCode MCParticleCreator::CreateTrackToMCParticleRelationships(cons
if( pMCRecoTrackerAssociationCollection.at(ic).getRec().id() != pTrack->getTrackerHits(ith).id() ) continue;
auto pSimHit = pMCRecoTrackerAssociationCollection.at(ic).getSim();
auto ipa = pSimHit.getMCParticle();
- if( m_id_pMC_map->find(ipa.id()) == m_id_pMC_map->end() ) continue;
+ if( m_id_pMC_map->find(ipa.id().index) == m_id_pMC_map->end() ) continue;
const float trueMomentum(pandora::CartesianVector(ipa.getMomentum()[0], ipa.getMomentum()[1], ipa.getMomentum()[2]).GetMagnitude());
const float deltaMomentum(std::fabs(recoMomentum - trueMomentum));
if (deltaMomentum < bestDeltaMomentum)
{
//pBestMCParticle =((*m_id_pMC_map)[ipa.id()]);
- best_mc_id = ipa.id() ;
+ best_mc_id = ipa.id().index ;
bestDeltaMomentum = deltaMomentum;
}
}
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp
index cfe8f2f25f1906723527534f6ec3154616c752c7..a95ea9c70b872d0731245f911af378d8e8a6dc30 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp
@@ -669,9 +669,9 @@ StatusCode PandoraPFAlg::CreateMCRecoParticleAssociation()
if(it->getRec().id() != hit.id()) continue;
for(auto itc = it->getSim().contributions_begin(); itc != it->getSim().contributions_end(); itc++)
{
- if(mc_map.find(itc->getParticle().id()) == mc_map.end()) mc_map[itc->getParticle().id()] = itc->getParticle() ;
- if(id_edep_map.find(itc->getParticle().id()) != id_edep_map.end()) id_edep_map[itc->getParticle().id()] = id_edep_map[itc->getParticle().id()] + itc->getEnergy() ;
- else id_edep_map[itc->getParticle().id()] = itc->getEnergy() ;
+ if(mc_map.find(itc->getParticle().id().index) == mc_map.end()) mc_map[itc->getParticle().id().index] = itc->getParticle() ;
+ if(id_edep_map.find(itc->getParticle().id().index) != id_edep_map.end()) id_edep_map[itc->getParticle().id().index] = id_edep_map[itc->getParticle().id().index] + itc->getEnergy() ;
+ else id_edep_map[itc->getParticle().id().index] = itc->getEnergy() ;
tot_en += itc->getEnergy() ;
}
}
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp
index 940d116f14d5682f6adfdd530deba7777a0c9d0d..8ebe8e239504f3475fa07c64db0199a72ea2e72b 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp
@@ -10,6 +10,14 @@
#include "edm4hep/Vertex.h"
#include "edm4hep/Vector3f.h"
#include "edm4hep/ReconstructedParticle.h"
+#if __has_include("edm4hep/EDM4hepVersion.h")
+#include "edm4hep/EDM4hepVersion.h"
+#else
+// Copy the necessary parts from the header above to make whatever we need to work here
+#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch)))
+// v00-09 is the last version without the capitalization change of the track vector members
+#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0)
+#endif
#include "gear/BField.h"
#include "gear/CalorimeterParameters.h"
@@ -313,7 +321,7 @@ pandora::StatusCode TrackCreator::ExtractKinks(const CollectionMaps& collectionM
for (unsigned int iTrack = 0, nTracks = pReconstructedParticle.tracks_size(); iTrack < nTracks; ++iTrack)
{
auto pTrack = pReconstructedParticle.getTracks(iTrack);
- (0 == iTrack) ? m_parentTrackList.insert(pTrack.id()) : m_daughterTrackList.insert(pTrack.id());
+ (0 == iTrack) ? m_parentTrackList.insert(pTrack.id().index) : m_daughterTrackList.insert(pTrack.id().index);
int trackPdgCode = pandora::UNKNOWN_PARTICLE_TYPE;
@@ -413,7 +421,7 @@ pandora::StatusCode TrackCreator::ExtractProngsAndSplits(const CollectionMaps& c
for (unsigned int iTrack = 0, nTracks = pReconstructedParticle.tracks_size(); iTrack < nTracks; ++iTrack)
{
edm4hep::Track pTrack = pReconstructedParticle.getTracks(iTrack);
- (0 == iTrack) ? m_parentTrackList.insert(pTrack.id()) : m_daughterTrackList.insert(pTrack.id());
+ (0 == iTrack) ? m_parentTrackList.insert(pTrack.id().index) : m_daughterTrackList.insert(pTrack.id().index);
if (0 == m_settings.m_shouldFormTrackRelationships) continue;
@@ -482,7 +490,7 @@ pandora::StatusCode TrackCreator::ExtractV0s(const CollectionMaps& collectionMap
for (unsigned int iTrack = 0, nTracks = pReconstructedParticle.tracks_size(); iTrack < nTracks; ++iTrack)
{
edm4hep::Track pTrack = pReconstructedParticle.getTracks(iTrack);
- m_v0TrackList.insert(pTrack.id());
+ m_v0TrackList.insert(pTrack.id().index);
int trackPdgCode = pandora::UNKNOWN_PARTICLE_TYPE;
@@ -538,7 +546,7 @@ bool TrackCreator::IsConflictingRelationship(const edm4hep::ReconstructedParticl
for (unsigned int iTrack = 0, nTracks = Particle.tracks_size(); iTrack < nTracks; ++iTrack)
{
edm4hep::Track pTrack = Particle.getTracks(iTrack) ;
- unsigned int pTrack_id = pTrack.id() ;
+ unsigned int pTrack_id = pTrack.id().index ;
if (this->IsDaughter(pTrack_id) || this->IsParent(pTrack_id) || this->IsV0(pTrack_id))
return true;
@@ -851,7 +859,7 @@ void TrackCreator::DefineTrackPfoUsage(const edm4hep::Track *const pTrack, Pando
bool canFormPfo(false);
bool canFormClusterlessPfo(false);
- if (trackParameters.m_reachesCalorimeter.Get() && !this->IsParent(pTrack->id()))
+ if (trackParameters.m_reachesCalorimeter.Get() && !this->IsParent(pTrack->id().index))
{
const float d0(std::fabs(pTrack->getTrackStates(0).D0)), z0(std::fabs(pTrack->getTrackStates(0).Z0));
@@ -880,8 +888,8 @@ void TrackCreator::DefineTrackPfoUsage(const edm4hep::Track *const pTrack, Pando
const float zCutForNonVertexTracks(m_tpcInnerR * std::fabs(pZ / pT) + m_settings.m_zCutForNonVertexTracks);
const bool passRzQualityCuts((zMin < zCutForNonVertexTracks) && (rInner < m_tpcInnerR + m_settings.m_maxTpcInnerRDistance));
- const bool isV0(this->IsV0(pTrack->id()));
- const bool isDaughter(this->IsDaughter(pTrack->id()));
+ const bool isV0(this->IsV0(pTrack->id().index));
+ const bool isDaughter(this->IsDaughter(pTrack->id().index));
// Decide whether track can be associated with a pandora cluster and used to form a charged PFO
if ((d0 < m_settings.m_d0TrackCut) && (z0 < m_settings.m_z0TrackCut) && (rInner < m_tpcInnerR + m_settings.m_maxTpcInnerRDistance))
@@ -918,7 +926,7 @@ void TrackCreator::DefineTrackPfoUsage(const edm4hep::Track *const pTrack, Pando
}
}
}
- else if (this->IsDaughter(pTrack->id()) || this->IsV0(pTrack->id()))
+ else if (this->IsDaughter(pTrack->id().index) || this->IsV0(pTrack->id().index))
{
std::cout<<"WARNING Recovering daughter or v0 track " << trackParameters.m_momentumAtDca.Get().GetMagnitude() << std::endl;
canFormPfo = true;
@@ -1020,9 +1028,15 @@ int TrackCreator::GetNTpcHits(const edm4hep::Track *const pTrack) const
//According to FG: [ 2 * lcio::ILDDetID::TPC - 2 ] is the first number and it is supposed to
//be the number of hits in the fit and this is what should be used !
// at least for DD4hep/DDSim
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ return pTrack->getSubdetectorHitNumbers(3);//FIXME https://github.com/wenxingfang/CEPCSW/blob/master/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp#L483
+ }
+ else return pTrack->getSubdetectorHitNumbers(2 * 4 - 1);// lcio::ILDDetID::TPC=4, still use LCIO code now
+#else
return pTrack->getSubDetectorHitNumbers(3);//FIXME https://github.com/wenxingfang/CEPCSW/blob/master/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp#L483
}
else return pTrack->getSubDetectorHitNumbers(2 * 4 - 1);// lcio::ILDDetID::TPC=4, still use LCIO code now
+#endif
}
//------------------------------------------------------------------------------------------------------------------------------------------
@@ -1036,9 +1050,16 @@ int TrackCreator::GetNFtdHits(const edm4hep::Track *const pTrack) const
// ---- use hitsInFit :
//return pTrack->getSubdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - 1 ];
if(m_settings.m_use_dd4hep_geo){
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ return pTrack->getSubdetectorHitNumbers(1);//FIXME https://github.com/wenxingfang/CEPCSW/blob/master/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp#L481
+ }
+ else return pTrack->getSubdetectorHitNumbers( 2 * 3 - 1 );// lcio::ILDDetID::FTD=3
+#else
return pTrack->getSubDetectorHitNumbers(1);//FIXME https://github.com/wenxingfang/CEPCSW/blob/master/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp#L481
}
else return pTrack->getSubDetectorHitNumbers( 2 * 3 - 1 );// lcio::ILDDetID::FTD=3
+
+#endif
}
//------------------------------------------------------------------------------------------------------------------------------------------
diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/MCParticleCreator.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/MCParticleCreator.cpp
index 40f44d1a4a5a3f76cf8f527c704d855eeb709ddb..2da13168523268d37056960708474d8891c12ad2 100644
--- a/Reconstruction/PFA/Pandora/MatrixPandora/src/MCParticleCreator.cpp
+++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/MCParticleCreator.cpp
@@ -57,7 +57,7 @@ pandora::StatusCode MCParticleCreator::CreateMCParticles(const CollectionMaps& c
mcParticleParameters.m_particleId = pMcParticle.getPDG();
mcParticleParameters.m_mcParticleType = pandora::MC_3D;
mcParticleParameters.m_pParentAddress = &pMcParticle;
- unsigned int p_id = pMcParticle.id();
+ unsigned int p_id = pMcParticle.id().index;
auto p_mc = &pMcParticle;
(*m_id_pMC_map) [p_id] = p_mc;
mcParticleParameters.m_momentum = pandora::CartesianVector(pMcParticle.getMomentum()[0], pMcParticle.getMomentum()[1],
@@ -271,8 +271,8 @@ pandora::StatusCode MCParticleCreator::CreateCaloHitToMCParticleRelationships(co
edm4hep::CaloHitContribution conb = pSimHit.getContributions(iCont);
auto ipa = conb.getParticle();
float ien = conb.getEnergy();
- if( m_id_pMC_map->find(ipa.id()) == m_id_pMC_map->end() ) continue;
- auto p_tmp = (*m_id_pMC_map)[ipa.id()];
+ if( m_id_pMC_map->find(ipa.id().index) == m_id_pMC_map->end() ) continue;
+ auto p_tmp = (*m_id_pMC_map)[ipa.id().index];
mcParticleToEnergyWeightMap[p_tmp] += ien;
}
@@ -331,12 +331,12 @@ pandora::StatusCode MCParticleCreator::CreateTrackToMCParticleRelationships(cons
if( pMCRecoTrackerAssociationCollection.at(ic).getRec().id() != pTrack->getTrackerHits(ith).id() ) continue;
auto pSimHit = pMCRecoTrackerAssociationCollection.at(ic).getSim();
auto ipa = pSimHit.getMCParticle();
- if( m_id_pMC_map->find(ipa.id()) == m_id_pMC_map->end() ) continue;
+ if( m_id_pMC_map->find(ipa.id().index) == m_id_pMC_map->end() ) continue;
const float trueMomentum(pandora::CartesianVector(ipa.getMomentum()[0], ipa.getMomentum()[1], ipa.getMomentum()[2]).GetMagnitude());
const float deltaMomentum(std::fabs(recoMomentum - trueMomentum));
if (deltaMomentum < bestDeltaMomentum)
{
- pBestMCParticle = ((*m_id_pMC_map)[ipa.id()]);
+ pBestMCParticle = ((*m_id_pMC_map)[ipa.id().index]);
bestDeltaMomentum = deltaMomentum;
}
}
diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp
index b0a1d7ca4fd4c57807960672cc588f5f3576cf7e..16996f6ca34d503c97f67350c1262ac7343bf45d 100644
--- a/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp
+++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp
@@ -591,7 +591,7 @@ StatusCode PandoraMatrixAlg::Ana()
{
if(reco_associa_col->at(j).getRec().id() != pReco.id() ) continue;
std::cout<<"MC pid ="<<reco_associa_col->at(j).getSim().getPDG()<<",weight="<<reco_associa_col->at(j).getWeight()<<", px="<<reco_associa_col->at(j).getSim().getMomentum()[0]<<", py="<<reco_associa_col->at(j).getSim().getMomentum()[1]<<",pz="<<reco_associa_col->at(j).getSim().getMomentum()[2]<<std::endl;
- tmp_mc_id .push_back(reco_associa_col->at(j).getSim().id());
+ tmp_mc_id .push_back(reco_associa_col->at(j).getSim().id().index);
tmp_mc_weight.push_back(reco_associa_col->at(j).getWeight());
}
m_pReco_mc_id .push_back(tmp_mc_id);
@@ -604,7 +604,7 @@ StatusCode PandoraMatrixAlg::Ana()
{
for(unsigned int i=0 ; i< MCParticle->size(); i++)
{
- m_mc_id .push_back(MCParticle->at(i).id());
+ m_mc_id .push_back(MCParticle->at(i).id().index);
m_mc_p_size.push_back(MCParticle->at(i).parents_size());
m_mc_pid .push_back(MCParticle->at(i).getPDG());
m_mc_mass .push_back(MCParticle->at(i).getMass());
@@ -655,9 +655,9 @@ StatusCode PandoraMatrixAlg::CreateMCRecoParticleAssociation()
if(it->getRec().id() != hit.id()) continue;
for(auto itc = it->getSim().contributions_begin(); itc != it->getSim().contributions_end(); itc++)
{
- if(mc_map.find(itc->getParticle().id()) == mc_map.end()) mc_map[itc->getParticle().id()] = itc->getParticle() ;
- if(id_edep_map.find(itc->getParticle().id()) != id_edep_map.end()) id_edep_map[itc->getParticle().id()] = id_edep_map[itc->getParticle().id()] + itc->getEnergy() ;
- else id_edep_map[itc->getParticle().id()] = itc->getEnergy() ;
+ if(mc_map.find(itc->getParticle().id().index) == mc_map.end()) mc_map[itc->getParticle().id().index] = itc->getParticle() ;
+ if(id_edep_map.find(itc->getParticle().id().index) != id_edep_map.end()) id_edep_map[itc->getParticle().id().index] = id_edep_map[itc->getParticle().id().index] + itc->getEnergy() ;
+ else id_edep_map[itc->getParticle().id().index] = itc->getEnergy() ;
tot_en += itc->getEnergy() ;
}
}
diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/TrackCreator.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/TrackCreator.cpp
index dd633f076d472d602dcf1406e37e68df69a50d32..c290e069311e26d5294033c7eaddf8bebdde34d9 100644
--- a/Reconstruction/PFA/Pandora/MatrixPandora/src/TrackCreator.cpp
+++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/TrackCreator.cpp
@@ -10,6 +10,14 @@
#include "edm4hep/Vertex.h"
#include "edm4hep/ReconstructedParticle.h"
+#if __has_include("edm4hep/EDM4hepVersion.h")
+#include "edm4hep/EDM4hepVersion.h"
+#else
+// Copy the necessary parts from the header above to make whatever we need to work here
+#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch)))
+// v00-09 is the last version without the capitalization change of the track vector members
+#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0)
+#endif
#include "gear/BField.h"
#include "gear/CalorimeterParameters.h"
@@ -177,7 +185,7 @@ pandora::StatusCode TrackCreator::ExtractKinks(const CollectionMaps& collectionM
for (unsigned int iTrack = 0, nTracks = pReconstructedParticle.tracks_size(); iTrack < nTracks; ++iTrack)
{
auto pTrack = pReconstructedParticle.getTracks(iTrack);
- (0 == iTrack) ? m_parentTrackList.insert(pTrack.id()) : m_daughterTrackList.insert(pTrack.id());
+ (0 == iTrack) ? m_parentTrackList.insert(pTrack.id().index) : m_daughterTrackList.insert(pTrack.id().index);
int trackPdgCode = pandora::UNKNOWN_PARTICLE_TYPE;
@@ -277,7 +285,7 @@ pandora::StatusCode TrackCreator::ExtractProngsAndSplits(const CollectionMaps& c
for (unsigned int iTrack = 0, nTracks = pReconstructedParticle.tracks_size(); iTrack < nTracks; ++iTrack)
{
auto pTrack = pReconstructedParticle.getTracks(iTrack);
- (0 == iTrack) ? m_parentTrackList.insert(pTrack.id()) : m_daughterTrackList.insert(pTrack.id());
+ (0 == iTrack) ? m_parentTrackList.insert(pTrack.id().index) : m_daughterTrackList.insert(pTrack.id().index);
if (0 == m_settings.m_shouldFormTrackRelationships) continue;
@@ -346,7 +354,7 @@ pandora::StatusCode TrackCreator::ExtractV0s(const CollectionMaps& collectionMap
for (unsigned int iTrack = 0, nTracks = pReconstructedParticle.tracks_size(); iTrack < nTracks; ++iTrack)
{
auto pTrack = pReconstructedParticle.getTracks(iTrack);
- m_v0TrackList.insert(pTrack.id());
+ m_v0TrackList.insert(pTrack.id().index);
int trackPdgCode = pandora::UNKNOWN_PARTICLE_TYPE;
@@ -401,7 +409,7 @@ bool TrackCreator::IsConflictingRelationship(const edm4hep::ReconstructedParticl
for (unsigned int iTrack = 0, nTracks = Particle.tracks_size(); iTrack < nTracks; ++iTrack)
{
edm4hep::Track pTrack = Particle.getTracks(iTrack) ;
- unsigned int pTrack_id = pTrack.id() ;
+ unsigned int pTrack_id = pTrack.id().index ;
if (this->IsDaughter(pTrack_id) || this->IsParent(pTrack_id) || this->IsV0(pTrack_id))
return true;
@@ -712,7 +720,7 @@ void TrackCreator::DefineTrackPfoUsage(edm4hep::Track *const pTrack, PandoraApi:
bool canFormPfo(false);
bool canFormClusterlessPfo(false);
- if (trackParameters.m_reachesCalorimeter.Get() && !this->IsParent(pTrack->id()))
+ if (trackParameters.m_reachesCalorimeter.Get() && !this->IsParent(pTrack->id().index))
{
const float d0(std::fabs(pTrack->getTrackStates(0).D0)), z0(std::fabs(pTrack->getTrackStates(0).Z0));
@@ -740,8 +748,8 @@ void TrackCreator::DefineTrackPfoUsage(edm4hep::Track *const pTrack, PandoraApi:
const float zCutForNonVertexTracks(m_tpcInnerR * std::fabs(pZ / pT) + m_settings.m_zCutForNonVertexTracks);
const bool passRzQualityCuts((zMin < zCutForNonVertexTracks) && (rInner < m_tpcInnerR + m_settings.m_maxTpcInnerRDistance));
- const bool isV0(this->IsV0(pTrack->id()));
- const bool isDaughter(this->IsDaughter(pTrack->id()));
+ const bool isV0(this->IsV0(pTrack->id().index));
+ const bool isDaughter(this->IsDaughter(pTrack->id().index));
// Decide whether track can be associated with a pandora cluster and used to form a charged PFO
if ((d0 < m_settings.m_d0TrackCut) && (z0 < m_settings.m_z0TrackCut) && (rInner < m_tpcInnerR + m_settings.m_maxTpcInnerRDistance))
@@ -778,7 +786,7 @@ void TrackCreator::DefineTrackPfoUsage(edm4hep::Track *const pTrack, PandoraApi:
}
}
}
- else if (this->IsDaughter(pTrack->id()) || this->IsV0(pTrack->id()))
+ else if (this->IsDaughter(pTrack->id().index) || this->IsV0(pTrack->id().index))
{
std::cout<<"WARNING Recovering daughter or v0 track " << trackParameters.m_momentumAtDca.Get().GetMagnitude() << std::endl;
canFormPfo = true;
@@ -870,14 +878,22 @@ bool TrackCreator::PassesQualityCuts(edm4hep::Track *const pTrack, const Pandora
int TrackCreator::GetNTpcHits(edm4hep::Track *const pTrack) const
{
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ return pTrack->getSubdetectorHitNumbers(2 * lcio::ILDDetID::TPC - 1);// still use LCIO code now
+#else
return pTrack->getSubDetectorHitNumbers(2 * lcio::ILDDetID::TPC - 1);// still use LCIO code now
+#endif
}
//------------------------------------------------------------------------------------------------------------------------------------------
int TrackCreator::GetNFtdHits(edm4hep::Track *const pTrack) const
{
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ return pTrack->getSubdetectorHitNumbers( 2 * lcio::ILDDetID::FTD - 1 );
+#else
return pTrack->getSubDetectorHitNumbers( 2 * lcio::ILDDetID::FTD - 1 );
+#endif
}
//------------------------------------------------------------------------------------------------------------------------------------------
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
index 1b154b957586939bffcfb4492eb1fb18fc5da86e..5e865cb162bd3b81bd9a7c90d7e24915d84d9050 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
@@ -67,6 +67,7 @@ const int GenfitTrack::s_PDG[2][5]
bool
sortDCSimHit(edm4hep::SimTrackerHit hit1,edm4hep::SimTrackerHit hit2)
+//sortDCSimHit(edm4hep::ConstSimTrackerHit hit1,edm4hep::ConstSimTrackerHit hit2)
{
//std::cout<<"hit1"<<hit1<<std::endl;
//std::cout<<"hit2"<<hit2<<std::endl;
@@ -75,11 +76,13 @@ sortDCSimHit(edm4hep::SimTrackerHit hit1,edm4hep::SimTrackerHit hit2)
}
bool
sortDCDigi(std::pair<double,edm4hep::TrackerHit*> hitPair1,std::pair<double,edm4hep::TrackerHit*> hitPair2)
+//sortDCDigi(std::pair<double,edm4hep::ConstTrackerHit*> hitPair1,std::pair<double,edm4hep::ConstTrackerHit*> hitPair2)
{
bool isEarly=hitPair1.first<hitPair2.first;
return isEarly;
}
//bool sortDCDigiLayer(std::pair<int,edm4hep::TrackerHit> hitPair1,std::pair<int,edm4hep::TrackerHit> hitPair2)
+//bool sortDCDigiLayer(std::pair<int,edm4hep::ConstTrackerHit> hitPair1,std::pair<int,edm4hep::ConstTrackerHit> hitPair2)
//{
// bool isEarly=hitPair1.first<hitPair2.first;
// return isEarly;
@@ -162,7 +165,7 @@ bool GenfitTrack::createGenfitTrackFromEDM4HepTrack(int pidType,
if(track.trackerHits_size()<=0) {
std::cout << " track.trackerHits_size = " << track.trackerHits_size() << std::endl;
if(m_debug){
- std::cout<<"createGenfitTrackFromEDM4HepTrack skip track w/o hit"<<std::endl;
+ std::cout<<"createGenfitTrackFromEDM4HepTrack skip track w/o hit"<<std::endl;
}
return false;
}
@@ -234,6 +237,7 @@ bool GenfitTrack::addSpacePointMeasurement(const TVector3& pos,
/// Return isurface of a silicon hit
const dd4hep::rec::ISurface*
GenfitTrack::getISurface(edm4hep::TrackerHit* hit){
+//GenfitTrack::getISurface(edm4hep::ConstTrackerHit* hit){
dd4hep::rec::SurfaceManager surfaceManager(*m_geomSvc->lcdd());
std::string detectorName;
@@ -280,6 +284,7 @@ GenfitTrack::getISurface(edm4hep::TrackerHit* hit){
/// Add a 1d strip or 2d pixel smeared by sigma
bool
GenfitTrack::addSiliconMeasurement(edm4hep::TrackerHit* hit,
+//GenfitTrack::addSiliconMeasurement(edm4hep::ConstTrackerHit* hit,
float sigmaU,float sigmaV,int cellID,int hitID)
{
if(m_debug>0) std::cout<<"addSiliconMeasurement "<<*hit<<std::endl;
@@ -362,6 +367,47 @@ int GenfitTrack::addSiliconMeasurements(edm4hep::Track& track,
return nHitAdd;
}
+//Add wire measurement on wire, unit conversion here
+int GenfitTrack::addWireMeasurementsFromListTrF(const edm4hep::TrackerHitCollection* trkHits,
+ float sigma,int sortMethod)
+{
+ double driftVelocity=40.;
+ std::vector<edm4hep::TrackerHit*> hits;
+ for(auto trkHit:*trkHits){
+ hits.push_back(&trkHit);
+ }
+
+ std::vector<edm4hep::TrackerHit*> sortedTrackerHits;
+ getSortedTrackerHitsTrF(hits,sortedTrackerHits,sortMethod);
+
+ if(m_debug>0){
+ std::cout<<"n sortedTrackerHits "<<sortedTrackerHits.size()<<std::endl;
+ }
+ int nHitAdd=0;
+ for(auto trackerHit : sortedTrackerHits){
+
+ TVector3 end0(0,0,0);
+ TVector3 end1(0,0,0);
+ getEndPointsOfWire(trackerHit->getCellID(),end0,end1);
+
+ ///New a TrackPoint,create connection between meas. and trackPoint
+ WireMeasurementDC* wireMeasurementDC=
+ new WireMeasurementDC(1e-3*driftVelocity*trackerHit->getTime(),sigma,
+ end0,end1,getDetTypeID(trackerHit->getCellID()),nHitAdd,nullptr);
+
+ int layer = m_decoderDC->get(trackerHit->getCellID(),"layer");
+ int cellID = m_decoderDC->get(trackerHit->getCellID(),"cellID");
+ const edm4hep::TrackerHit trackerHit_;
+ wireMeasurementDC->setTrackerHit(trackerHit_,layer,cellID);
+
+ m_track->insertPoint(new genfit::TrackPoint(wireMeasurementDC,m_track));
+ if(m_debug>=2){ std::cout<<nHitAdd-1; wireMeasurementDC->Print(); }
+ nHitAdd++;
+ }//end of loop over sotred hits
+ return nHitAdd;
+}//end of addWireMeasurementsFromList
+
+
//Add wire measurement on wire, unit conversion here
//int GenfitTrack::addWireMeasurementsFromList(podio::RelationRange<edm4hep::TrackerHit> hits,float sigma,
int GenfitTrack::addWireMeasurementsFromList(std::vector<edm4hep::TrackerHit*>& hits,float sigma,
@@ -371,11 +417,11 @@ int GenfitTrack::addWireMeasurementsFromList(std::vector<edm4hep::TrackerHit*>&
if(m_debug>0){ std::cout<<"addWireMeasurementsFromList"<<std::endl; }
//podio::RelationRange<edm4hep::TrackerHit> hits_t=track.getTrackerHits();
std::vector<edm4hep::TrackerHit*> sortedTrackerHits;
-// sortedTrackerHits.reserve(100);
+ // sortedTrackerHits.reserve(100);
getSortedTrackerHits(hits,assoHits,sortedTrackerHits,sortMethod);
if(m_debug>0){
- std::cout<<"n sortedTrackerHits "<<sortedTrackerHits.size()<<std::endl;
+ std::cout<<"n sortedTrackerHits "<<sortedTrackerHits.size()<<std::endl;
}
int nHitAdd=0;
for(auto trackerHit : sortedTrackerHits){
@@ -385,6 +431,7 @@ int GenfitTrack::addWireMeasurementsFromList(std::vector<edm4hep::TrackerHit*>&
GenfitHit* genfitHit=
makeAGenfitHit(trackerHit,&simTrackerHitAsso,sigma,truthAmbig,
skipCorner,skipNear);
+
if(nullptr==genfitHit) continue;
m_genfitHitVec.push_back(genfitHit);
@@ -397,6 +444,7 @@ int GenfitTrack::addWireMeasurementsFromList(std::vector<edm4hep::TrackerHit*>&
return nHitAdd;
}//end of addWireMeasurementsFromList
+
//Add wire measurement on wire, unit conversion here
int GenfitTrack::addWireMeasurementsOnTrack(edm4hep::Track& track,float sigma,
const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
@@ -405,7 +453,7 @@ int GenfitTrack::addWireMeasurementsOnTrack(edm4hep::Track& track,float sigma,
if(m_debug>0){ std::cout<<"addWireMeasurementsOnTrack"<<std::endl; }
if(m_debug>0){
- std::cout<<"n trackerHits size"<<track.trackerHits_size()<<std::endl;
+ std::cout<<"n trackerHits size"<<track.trackerHits_size()<<std::endl;
}
int nHitAdd=0;
if(0==track.trackerHits_size()){
@@ -486,8 +534,8 @@ unsigned int GenfitTrack::getNumPoints() const
}
GenfitHit* GenfitTrack::GetHit(long unsigned int i) const {
- if(i>=m_genfitHitVec.size())return nullptr;
- return m_genfitHitVec[i];
+ if(i>=m_genfitHitVec.size())return nullptr;
+ return m_genfitHitVec[i];
}
unsigned int GenfitTrack::getNumPointsDet(int detTypeID) const
@@ -819,8 +867,8 @@ int GenfitTrack::addSpacePointsDC(const edm4hep::Track& track,
std::vector<float> sigmaU,std::vector<float> sigmaV)
{
if(m_debug>=2){
- std::cout<<m_name<<" 5. addSpacePointsDC nTrackerHit="
- <<track.trackerHits_size()<<std::endl;
+ std::cout<<m_name<<" 5. addSpacePointsDC nTrackerHit="
+ <<track.trackerHits_size()<<std::endl;
}
///Get TrackerHit with min. time in Cell
@@ -1079,7 +1127,9 @@ bool GenfitTrack::storeTrack(edm4hep::MutableReconstructedParticle& recParticle,
const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
std::vector<double>& driftDis,
std::vector<double>& FittedDoca,
- std::vector<double>& Res)
+ std::vector<double>& truthDoca,
+ std::vector<double>& Res,
+ std::vector<double>& truthRes)
{
int id = 0;
@@ -1092,28 +1142,46 @@ bool GenfitTrack::storeTrack(edm4hep::MutableReconstructedParticle& recParticle,
// getNumRawMeasurements
unsigned int nPoints = m_track->getNumPoints();
std::cout << " nPoints = " << nPoints << std::endl;
+ unsigned int nPointsWithMea = m_track->getNumPointsWithMeasurement();
+ std::cout << " nPointsWithMea = " << nPointsWithMea << std::endl;
std::vector<double> hitMomMag;
+ //std::cout << __FILE__ << " " << __LINE__ << std::endl;
while(1){
genfit::TrackPoint* point = m_track->getPointWithFitterInfo(id);
+ //std::cout << __FILE__ << " " << __LINE__ << std::endl;
if(!point)break;
+ //std::cout << __FILE__ << " " << __LINE__ << " id = " << id << std::endl;
id++;
genfit::AbsMeasurement* absMea = point->getRawMeasurement();
+ //std::cout << __FILE__ << " " << __LINE__ << " absMea: " << std::endl;
+ //absMea->Print();
// getPoint FitterInfo
genfit::AbsFitterInfo* FitterInfo = point->getFitterInfo();
+ //std::cout << __FILE__ << " " << __LINE__ << " FitterInfo: " << std::endl;
+ //FitterInfo->Print();
genfit::KalmanFitterInfo *KalmanfitterInfo =
dynamic_cast<genfit::KalmanFitterInfo*>(FitterInfo);
+ //std::cout << __FILE__ << " " << __LINE__ << " KalmanfitterInfo: " << std::endl;
+ //KalmanfitterInfo->Print();
+
unsigned int nNumMea = KalmanfitterInfo->getNumMeasurements();
+ //std::cout << __FILE__ << " " << __LINE__ << " nNumMea = " << nNumMea << std::endl;
bool flag = false;
for(int i=0;i<nNumMea;i++)
{
genfit::MeasurementOnPlane* MeaOnPlane =
KalmanfitterInfo->getMeasurementOnPlane(i);
+ //std::cout << __FILE__ << " " << __LINE__ << " MeaOnPlane: " << std::endl;
+ //MeaOnPlane->Print();
+
double weight = MeaOnPlane->getWeight();
+ std::cout << __FILE__ << " " << __LINE__ << " weight = " << weight << std::endl;
+
if(weight>0.8) flag = true;
}
if(flag) fitid++;
@@ -1133,6 +1201,7 @@ bool GenfitTrack::storeTrack(edm4hep::MutableReconstructedParticle& recParticle,
dynamic_cast<WireMeasurementDC*>(absMea);
if(dcMea){
const edm4hep::TrackerHit* TrackerHit_ = dcMea->getTrackerHit();
+ if(nullptr==TrackerHit_) std::cout << __LINE__ << std::endl;
track.addToTrackerHits(*TrackerHit_);
@@ -1142,33 +1211,14 @@ bool GenfitTrack::storeTrack(edm4hep::MutableReconstructedParticle& recParticle,
truthMomEdep.push_back(TrackerHit_->getEDep());
double DriftDis,fittedDoca,res = 0;
- GetDocaRes(dcFit,DriftDis,fittedDoca,res);
- driftDis.push_back(DriftDis);
- FittedDoca.push_back(fittedDoca);
+ GetDocaRes((id-1),DriftDis,fittedDoca,res);
+ driftDis.push_back(10*DriftDis);
+ FittedDoca.push_back(10*fittedDoca);
+ truthDoca.push_back(40*1e-3*TrackerHit_->getTime());
+ truthRes.push_back((40*1e-3*TrackerHit_->getTime())-10*fittedDoca);
Res.push_back(res);
- // TMatrixDSym fittedHitCov(6);//cm, GeV
- // TLorentzVector fittedHitPos;
- // TVector3 fittedHitMom;
- // int fittedState=getFittedState(fittedHitPos,fittedHitMom,fittedHitCov,dcFit);
- // hitMomMag.push_back(fittedHitMom.Mag());
- //
- // for(int iSim=0;iSim<assoHits->size();iSim++)
- // {
- // //if(*TrackerHit_ == assoHits->at(iSim).getRec())
- // if(TrackerHit_->getCellID() == assoHits->at(iSim).getRec().getCellID())
- // {
- // // std::cout << " if TrackerHit_ = " << TrackerHit_->getCellID() << std::endl;
- // // std::cout << " if assoHits->at(iSim).getRec() = " << assoHits->at(iSim).getRec().getCellID() << std::endl;
- // // std::cout << " Sim Mom = " << assoHits->at(iSim).getSim().getMomentum() << std::endl;
- // // std::cout << " Sim MomMag = " << sqrt(assoHits->at(iSim).getSim().getMomentum()[0]*assoHits->at(iSim).getSim().getMomentum()[0]+assoHits->at(iSim).getSim().getMomentum()[1]*assoHits->at(iSim).getSim().getMomentum()[1]+assoHits->at(iSim).getSim().getMomentum()[2]*assoHits->at(iSim).getSim().getMomentum()[2]) << std::endl;
- //
- // break;
- // }
- // }
-
- //std::cout << " i = " << dcFit << "fittedMom = " << fittedHitMom.X() << " " << fittedHitMom.Y() << " " << fittedHitMom.Z() << std::endl;
- //std::cout << " i = " << dcFit << "fittedMomMag = " << fittedHitMom.Mag() << std::endl;
+
dcFit++;
}
}
@@ -1183,10 +1233,12 @@ bool GenfitTrack::storeTrack(edm4hep::MutableReconstructedParticle& recParticle,
std::cout<<"nDC: "<<nFittedDC<<", nSDT: "<<nFittedSDT<<std::endl;
std::cout<<"nFittedDC: "<<dcFit<<", nFittedSDT: "<<sdtFit<<std::endl;
+
+ if(fitid<1e-9) return false;
+
if(m_debug>0)std::cout<<m_name<<" store track ndfCut "<<ndfCut<<" chi2Cut "
<<chi2Cut<<std::endl;
-
/// Get fit status
const genfit::FitStatus* fitState = getFitStatus();
double ndf = fitState->getNdf();
@@ -1499,6 +1551,35 @@ void GenfitTrack::getTrackFromEDMTrack(const edm4hep::Track& edm4HepTrack,
}
}
+void GenfitTrack::getTrackFromEDMTrackFinding(const edm4hep::Track& edm4HepTrack,
+ double& charge, TVectorD& trackParam, TMatrixDSym& cov,TVector3& pos,
+ TVector3& mom){
+ //double Bz=m_genfitField->getBz(TVector3{0.,0.,0.})/GenfitUnit::tesla;
+ // FIXME
+ //double BZ=GenfitField::getBzFinding(TVector3{0.,0.,0.});
+ double Bz=3*GenfitUnit::tesla;
+ double charge_double;
+ CEPC::getPosMomFromTrackState(edm4HepTrack.getTrackStates(1),Bz,pos,mom,charge_double,cov);
+
+ //std::cout<<__LINE__<<" Bz "<<Bz<<" charge "<<charge_double<<std::endl;
+ //pos.Print();
+ //mom.Print();
+ charge=(int) charge_double;
+ trackParam[0]=pos[0]*GenfitUnit::mm;
+ trackParam[1]=pos[1]*GenfitUnit::mm;
+ trackParam[2]=pos[2]*GenfitUnit::mm;
+ trackParam[3]=mom[0]*GenfitUnit::GeV;
+ trackParam[4]=mom[1]*GenfitUnit::GeV;
+ trackParam[5]=mom[2]*GenfitUnit::GeV;
+
+ //cov unit conversion
+ for(int i=0;i<6;i++){
+ for(int j=0;j<6;j++){
+ cov(i,j) = cov(i,j)*GenfitUnit::mm;
+ }
+ }
+}
+
//to genfit unit
void GenfitTrack::getISurfaceOUV(const dd4hep::rec::ISurface* iSurface,TVector3& o,
TVector3& u,TVector3& v, double& lengthU,double& lengthV){
@@ -1570,6 +1651,51 @@ bool GenfitTrack::isCDCHit(edm4hep::TrackerHit* hit){
getDetTypeID(hit->getCellID());
}
+// TrackFinding
+void GenfitTrack::getSortedTrackerHitsTrF(
+ std::vector<edm4hep::TrackerHit*> trackerHits,
+ std::vector<edm4hep::TrackerHit*>& sortedDCTrackerHits,
+ int sortMethod){
+
+ std::vector<std::pair<double,edm4hep::TrackerHit*> > sortedDCTrackerHitPair;
+ for(auto trackerHit : trackerHits){
+ //edm4hep::TrackerHit* thisHit = trackerHit;
+ //if(!isCDCHit(thisHit))continue;//skip non-DC trackerHit
+
+ double time=trackerHit->getTime();
+ if(0==sortMethod){
+ //by time
+ sortedDCTrackerHitPair.push_back(std::make_pair(time,trackerHit));
+ if(m_debug>0){ std::cout<<"sorted DC digi by time"<<std::endl;}
+ }else if(1==sortMethod){
+ //by layer
+ sortedDCTrackerHitPair.push_back(std::make_pair(
+ m_decoderDC->get(trackerHit->getCellID(),"layer"),trackerHit));
+ if(m_debug>0){ std::cout<<"sorted DC digi by layer"<<std::endl;}
+ }else{
+ sortedDCTrackerHits.push_back(trackerHit);
+ }
+ }
+ if(0==sortMethod || 1==sortMethod){
+ std::sort(sortedDCTrackerHitPair.begin(),sortedDCTrackerHitPair.end(),
+ sortDCDigi);
+ for(auto trackerHit:sortedDCTrackerHitPair){
+ sortedDCTrackerHits.push_back(trackerHit.second);
+ }
+ }
+ if(m_debug>0){
+ std::cout<<"trackerHits on track after sort\n";
+ for(auto trackerHit:sortedDCTrackerHits){
+ std::cout<<trackerHit->getCellID()<<"("<<std::setw(2)
+ <<m_decoderDC->get(trackerHit->getCellID(),"layer")
+ <<","<<std::setw(3)
+ <<m_decoderDC->get(trackerHit->getCellID(),"cellID")<<")\n";
+ }
+ std::cout<<"\n"; std::cout.unsetf(std::ios_base::floatfield);
+ }
+ return;
+}//end of getSortedTrackerHits
+
void GenfitTrack::getSortedTrackerHits(
std::vector<edm4hep::TrackerHit*>& trackerHits,
const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
index 27ef200b6153e0354b80d67067e238d4ec664c2b..8a3e914a4971874d3bc2172684c2e41d14e19f2b 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
@@ -71,7 +71,7 @@ namespace dd4hep {
class GenfitTrack {
friend int GenfitFitter::processTrack(
GenfitTrack* track, bool resort);
-
+
friend int GenfitFitter::processTrackWithRep(
GenfitTrack* track, int repID, bool resort);
@@ -102,7 +102,6 @@ class GenfitTrack {
const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
std::vector<float> sigmaU,std::vector<float> sigmaV);
-
///Add WireMeasurements of hits on track
virtual int addWireMeasurementsOnTrack(edm4hep::Track& track,float sigma,
const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
@@ -113,6 +112,9 @@ class GenfitTrack {
const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
int sortMethod,bool truthAmbig,float skipCorner, float skipNear);
+ virtual int addWireMeasurementsFromListTrF(const edm4hep::TrackerHitCollection* trkHits,
+ float sigma,int sortMethod);
+
///Add one silicon hits
bool addSiliconMeasurement(edm4hep::TrackerHit* hit,
float sigmaU,float sigmaV,int cellID,int hitID);
@@ -141,7 +143,9 @@ class GenfitTrack {
const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
std::vector<double>& driftDis,
std::vector<double>& FittedDoca,
- std::vector<double>& Res);
+ std::vector<double>& truthDoca,
+ std::vector<double>& Res,
+ std::vector<double>& truthRes);
///A tool to convert track to the first layer of DC
void pivotToFirstLayer(const edm4hep::Vector3d& pos,
@@ -221,6 +225,10 @@ class GenfitTrack {
/// Get a hit according to index
GenfitHit* GetHit(long unsigned int i) const;
+
+ static void getTrackFromEDMTrackFinding(const edm4hep::Track& edm4HepTrack,
+ double& charge, TVectorD& trackParam, TMatrixDSym& cov,TVector3& pos,
+ TVector3& mom);
protected:
//genfit::Track* getTrack() {return m_track;}
@@ -265,6 +273,9 @@ class GenfitTrack {
const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
std::vector<edm4hep::TrackerHit*>& sortedDCTrackerHits,
int sortMethod);
+ void getSortedTrackerHitsTrF(std::vector<edm4hep::TrackerHit*> trackerHits,
+ std::vector<edm4hep::TrackerHit*>& sortedDCTrackerHits,
+ int sortMethod=1);
genfit::Track* m_track;/// track
int m_debug;/// debug level
diff --git a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp
index df95acc6cbe9cfa3fe53902b659792f367b5c558..3eb06e29713317bc809ef615a31e801259e8fd56 100644
--- a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp
+++ b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp
@@ -6,6 +6,14 @@
#include "edm4hep/TrackerHit.h"
#include "edm4hep/TrackerHit.h"
#include "edm4hep/Track.h"
+#if __has_include("edm4hep/EDM4hepVersion.h")
+#include "edm4hep/EDM4hepVersion.h"
+#else
+// Copy the necessary parts from the header above to make whatever we need to work here
+#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch)))
+// v00-09 is the last version without the capitalization change of the track vector members
+#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0)
+#endif
#include "UTIL/ILDConf.h"
@@ -1001,13 +1009,21 @@ void ForwardTrackingAlg::finaliseTrack( edm4hep::MutableTrack* trackImpl ){
//trackImpl->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::TPC - 1 ] = hitNumbers[lcio::ILDDetID::TPC];
//trackImpl->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::SET - 1 ] = hitNumbers[lcio::ILDDetID::SET];
//trackImpl->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::ETD - 1 ] = hitNumbers[lcio::ILDDetID::ETD];
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ trackImpl->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::VXD]);
+ trackImpl->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SIT]);
+ trackImpl->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::FTD]);
+ trackImpl->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::TPC]);
+ trackImpl->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SET]);
+ trackImpl->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::ETD]);
+#else
trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::VXD]);
trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SIT]);
trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::FTD]);
trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::TPC]);
trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SET]);
trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::ETD]);
-
+#endif
return;
}
diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp
index f9b6a5466aba7cf08527d77d2f3b6c5d429ba3d8..1ef104fa4b479dde6ce3be62b98ad6a6feb31b4a 100644
--- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp
+++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp
@@ -9,6 +9,14 @@
//#include "edm4hep/TrackerHitPlane.h"
#include "edm4hep/Track.h"
#include "edm4hep/TrackState.h"
+#if __has_include("edm4hep/EDM4hepVersion.h")
+#include "edm4hep/EDM4hepVersion.h"
+#else
+// Copy the necessary parts from the header above to make whatever we need to work here
+#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch)))
+// v00-09 is the last version without the capitalization change of the track vector members
+#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0)
+#endif
#include <iostream>
#include <algorithm>
@@ -2833,10 +2841,15 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
MarlinTrk::addHitNumbersToTrack(&track, all_hits, false, cellID_encoder);
delete marlinTrk;
-
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ int nhits_in_vxd = track.getSubdetectorHitNumbers(0);
+ int nhits_in_ftd = track.getSubdetectorHitNumbers(1);
+ int nhits_in_sit = track.getSubdetectorHitNumbers(2);
+#else
int nhits_in_vxd = track.getSubDetectorHitNumbers(0);
int nhits_in_ftd = track.getSubDetectorHitNumbers(1);
int nhits_in_sit = track.getSubDetectorHitNumbers(2);
+#endif
//debug() << " Hit numbers for Track "<< track.id() << ": "
debug() << " Hit numbers for Track "<< iTrk <<": "
diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
old mode 100755
new mode 100644
index 08e1fb022fdfc240dd127957085fb954ca2edf54..ee5f1feb25c5cda5c95f0f4d9dab9801e943db04
--- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
+++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
@@ -8,6 +8,14 @@
#include <edm4hep/TrackerHit.h>
#include <edm4hep/TrackerHit.h>
#include <edm4hep/Track.h>
+#if __has_include("edm4hep/EDM4hepVersion.h")
+#include "edm4hep/EDM4hepVersion.h"
+#else
+// Copy the necessary parts from the header above to make whatever we need to work here
+#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch)))
+// v00-09 is the last version without the capitalization change of the track vector members
+#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0)
+#endif
#include <iostream>
#include <algorithm>
@@ -513,11 +521,19 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
float z0TrkCand = trkCand->getZ0();
// float phi0TrkCand = trkCand->getPhi();
// FIXME, fucd
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ int nhits_in_vxd = track.getSubdetectorHitNumbers(0);
+ int nhits_in_ftd = track.getSubdetectorHitNumbers(1);
+ int nhits_in_sit = track.getSubdetectorHitNumbers(2);
+ int nhits_in_tpc = track.getSubdetectorHitNumbers(3);
+ int nhits_in_set = track.getSubdetectorHitNumbers(4);
+#else
int nhits_in_vxd = track.getSubDetectorHitNumbers(0);
int nhits_in_ftd = track.getSubDetectorHitNumbers(1);
int nhits_in_sit = track.getSubDetectorHitNumbers(2);
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 ];
@@ -1150,13 +1166,18 @@ void FullLDCTrackingAlg::prepareVectors() {
trackExt->setNDF(tpcTrack.getNdf());
trackExt->setChi2(tpcTrack.getChi2());
for (int iHit=0;iHit<nHits;++iHit) {
- edm4hep::TrackerHit hit = tpcTrack.getTrackerHits(iHit);//hitVec[iHit];
- if(!hit.isAvailable()) error() << "Tracker hit not available" << endmsg;
+ edm4hep::TrackerHit hit = tpcTrack.getTrackerHits(iHit);
+ if (!hit.isAvailable()) {
+ error() << "Tracker hit not available" << endmsg;
+ continue;
+ }
//info() << "hit " << hit.id() << " " << hit.getCellID() << " " << hit.getPosition()[0] << " " << hit.getPosition()[1] << " " << hit.getPosition()[2] << endmsg;
auto it = mapTrackerHits.find(hit);
- if(it==mapTrackerHits.end()) error() << "Cannot find hit " << hit.id() << " in map" << endmsg;
- else continue;
- TrackerHitExtended * hitExt = it->second;
+ if (it==mapTrackerHits.end()) {
+ error() << "Cannot find hit " << hit.id() << " in map" << endmsg;
+ continue;
+ }
+ TrackerHitExtended* hitExt = it->second;
//info() << hit.id() << " " << hitExt << endmsg;
hitExt->setTrackExtended( trackExt );
trackExt->addTrackerHitExtended( hitExt );
@@ -1215,8 +1236,17 @@ void FullLDCTrackingAlg::prepareVectors() {
char strg[200];
HelixClass helixSi;
for (int iHit=0;iHit<nHits;++iHit) {
- edm4hep::TrackerHit hit = siTrack.getTrackerHits(iHit);//hitVec[iHit];
- TrackerHitExtended * hitExt = mapTrackerHits[hit];
+ edm4hep::TrackerHit hit = siTrack.getTrackerHits(iHit);
+ if (!hit.isAvailable()) {
+ error() << "Tracker hit not available" << endmsg;
+ continue;
+ }
+ auto it = mapTrackerHits.find(hit);
+ if (it==mapTrackerHits.end()) {
+ error() << "Cannot find hit " << hit.id() << " in map" << endmsg;
+ continue;
+ }
+ TrackerHitExtended* hitExt = it->second;
hitExt->setTrackExtended( trackExt );
trackExt->addTrackerHitExtended( hitExt );
@@ -1518,8 +1548,8 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
int nTPCHits = int(tpcHitVec.size());
int nHits = nTPCHits + nSiHits;
- //std::cout << "FullLDCTrackingAlg::CombineTracks nSiHits = " << nSiHits << endmsg;
- //std::cout << "FullLDCTrackingAlg::CombineTracks nTPCHits = " << nTPCHits << endmsg;
+ //debug() << "FullLDCTrackingAlg::CombineTracks nSiHits = " << nSiHits << endmsg;
+ //debug() << "FullLDCTrackingAlg::CombineTracks nTPCHits = " << nTPCHits << endmsg;
TrackerHitVec trkHits;
trkHits.reserve(nHits);
@@ -1732,7 +1762,7 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
tpcHitInFit.push_back(tpcHitVec[i]);
}
}
-
+
debug() << "FullLDCTrackingAlg::CombineTracks: Check for Silicon Hit rejections ... " << endmsg;
if ( (int)siOutliers.size() > _maxAllowedSiHitRejectionsForTrackCombination ) {
diff --git a/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc b/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc
index ffedd64edb8558338e69a028ddcdcbf8021d42e1..3ca601461c16443ca6352740f3358c3a4eb499cb 100644
--- a/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc
+++ b/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc
@@ -26,6 +26,8 @@
#include "gear/GEAR.h"
#include "gear/BField.h"
+#include "podio/podioVersion.h"
+
//#include "streamlog/streamlog.h"
@@ -76,8 +78,12 @@ namespace MarlinTrk {
MarlinKalTestTrack::MarlinKalTestTrack(MarlinKalTest* ktest)
- : _ktest(ktest), _trackHitAtPositiveNDF(edm4hep::TrackerHit(0)) {
-
+ : _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() ;
diff --git a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
index ffa0f7f9395accd360f1efce633308447873d258..fcacb19cf8fead3e3204d70af7517f286372faef 100644
--- a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
+++ b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
@@ -16,6 +16,15 @@
#include "edm4hep/Track.h"
#include "edm4hep/MutableTrack.h"
+#if __has_include("edm4hep/EDM4hepVersion.h")
+#include "edm4hep/EDM4hepVersion.h"
+#else
+// Copy the necessary parts from the header above to make whatever we need to work here
+#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch)))
+// v00-09 is the last version without the capitalization change of the track vector members
+#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0)
+#endif
+
#include <UTIL/BitField64.h>
#include <UTIL/ILDConf.h>
#include <UTIL/BitSet32.h>
@@ -24,6 +33,8 @@
#include "TMatrixD.h"
+#include "podio/podioVersion.h"
+
#define MIN_NDF 6
namespace MarlinTrk {
@@ -258,7 +269,9 @@ namespace MarlinTrk {
// set the initial track parameters
///////////////////////////////////////////////////////
- return_error = marlinTrk->initialise( *pre_fit, bfield_z, IMarlinTrack::backward ) ;
+ //return_error = marlinTrk->initialise( *pre_fit, bfield_z, IMarlinTrack::backward ) ;
+ // fucd: previous fixed as IMarlinTrack::backward, can not change? using input value now
+ return_error = marlinTrk->initialise( *pre_fit, bfield_z, fit_backwards ) ;
if (return_error != IMarlinTrack::success) {
std::cout << "MarlinTrk::createFit Initialisation of track fit failed with error : " << return_error << std::endl;
@@ -388,10 +401,10 @@ namespace MarlinTrk {
return_error = marlintrk->getNDF(ndf);
if ( return_error != IMarlinTrack::success) {
- //streamlog_out(DEBUG3) << "MarlinTrk::finaliseLCIOTrack: getNDF returns " << return_error << std::endl;
+ //std::cout << "MarlinTrk::finaliseLCIOTrack: getNDF returns " << return_error << std::endl;
return return_error;
} else if( ndf < 0 ) {
- //streamlog_out(DEBUG2) << "MarlinTrk::finaliseLCIOTrack: number of degrees of freedom less than 0 track dropped : NDF = " << ndf << std::endl;
+ //std::cout << "MarlinTrk::finaliseLCIOTrack: number of degrees of freedom less than 0 track dropped : NDF = " << ndf << std::endl;
return IMarlinTrack::error;
} else {
//streamlog_out(DEBUG1) << "MarlinTrk::finaliseLCIOTrack: NDF = " << ndf << std::endl;
@@ -481,9 +494,14 @@ namespace MarlinTrk {
///////////////////////////////////////////////////////
edm4hep::TrackState* trkStateAtLastHit = new edm4hep::TrackState() ;
+
edm4hep::TrackerHit lastHit = hits_in_fit.front().first;
+#if PODIO_BUILD_VERSION < PODIO_VERSION(0, 17, 4)
edm4hep::TrackerHit last_constrained_hit(0);// = 0 ;
+#else
+ auto last_constrained_hit = edm4hep::TrackerHit::makeEmpty();
+#endif
marlintrk->getTrackerHitAtPositiveNDF(last_constrained_hit);
return_error = marlintrk->smooth(lastHit);
@@ -543,7 +561,6 @@ namespace MarlinTrk {
track->addToTrackStates(*trkStateAtFirstHit);
} else {
//streamlog_out( WARNING ) << " >>>>>>>>>>> MarlinTrk::finaliseLCIOTrack: could not get TrackState at First Hit " << firstHit << std::endl ;
- delete trkStateAtFirstHit;
}
double r_first = firstHit.getPosition()[0]*firstHit.getPosition()[0] + firstHit.getPosition()[1]*firstHit.getPosition()[1];
@@ -569,7 +586,7 @@ namespace MarlinTrk {
track->addToTrackStates(*trkStateAtLastHit);
} else {
std::cout << "ERROR>>>>>>>>>>> MarlinTrk::finaliseLCIOTrack: could not get TrackState at Last Hit " << last_constrained_hit.id() << std::endl ;
- delete trkStateAtLastHit;
+ //delete trkStateAtLastHit;
}
// const EVENT::FloatVec& ma = trkStateAtLastHit->getCovMatrix();
@@ -587,17 +604,22 @@ namespace MarlinTrk {
// return_error = createTrackStateAtCaloFace(marlintrk, trkStateCalo, lastHit, tanL_is_positive);
if ( return_error == 0 ) {
+ //std::cout << "fucdout referencePoint " << trkStateCalo.referencePoint << std::endl;
trkStateCalo.location = MarlinTrk::Location::AtCalorimeter;
track->addToTrackStates(trkStateCalo);
} else {
- //streamlog_out( WARNING ) << " >>>>>>>>>>> MarlinTrk::finaliseLCIOTrack: could not get TrackState at Calo Face " << std::endl ;
+ std::cout << " >>>>>>>>>>> MarlinTrk::finaliseLCIOTrack: could not get TrackState at Calo Face " << std::endl ;
//delete trkStateCalo;
}
} else {
track->addToTrackStates(*atLastHit);
track->addToTrackStates(*atCaloFace);
+ //delete trkStateAtLastHit;
}
-
+
+ if(trkStateAtFirstHit) delete trkStateAtFirstHit;
+ if(trkStateAtLastHit) delete trkStateAtLastHit;
+ if(trkStateIP) delete trkStateIP;
///////////////////////////////////////////////////////
// done
///////////////////////////////////////////////////////
@@ -688,12 +710,21 @@ namespace MarlinTrk {
if ( hits_in_fit == false ) { // all hit atributed by patrec
offset = 1 ;
}
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::VXD]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::FTD]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SIT]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::TPC]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SET]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::ETD]);
+#else
track->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::VXD]);
track->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::FTD]);
track->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SIT]);
track->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::TPC]);
track->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SET]);
track->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::ETD]);
+#endif
//track->subdetectorHitNumbers().resize(2 * lcio::ILDDetID::ETD);
//track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::VXD - offset ] = hitNumbers[lcio::ILDDetID::VXD];
//track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - offset ] = hitNumbers[lcio::ILDDetID::FTD];
@@ -732,12 +763,21 @@ namespace MarlinTrk {
if ( hits_in_fit == false ) { // all hit atributed by patrec
offset = 1 ;
}
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::VXD]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::FTD]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SIT]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::TPC]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SET]);
+ track->addToSubdetectorHitNumbers(hitNumbers[UTIL::ILDDetID::ETD]);
+#else
track->addToSubDetectorHitNumbers(hitNumbers[lcio::ILDDetID::VXD]);
track->addToSubDetectorHitNumbers(hitNumbers[lcio::ILDDetID::FTD]);
track->addToSubDetectorHitNumbers(hitNumbers[lcio::ILDDetID::SIT]);
track->addToSubDetectorHitNumbers(hitNumbers[lcio::ILDDetID::TPC]);
track->addToSubDetectorHitNumbers(hitNumbers[lcio::ILDDetID::SET]);
track->addToSubDetectorHitNumbers(hitNumbers[lcio::ILDDetID::ETD]);
+#endif
//track->subdetectorHitNumbers().resize(2 * lcio::ILDDetID::ETD);
//track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::VXD - offset ] = hitNumbers[lcio::ILDDetID::VXD];
//track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - offset ] = hitNumbers[lcio::ILDDetID::FTD];
diff --git a/Simulation/DetSimAna/CMakeLists.txt b/Simulation/DetSimAna/CMakeLists.txt
index 883dfd8215c9ccfd9e88b7d1713c442de1f5e151..b1bd935f32e1de1bb87d8c490f5967960deda351 100644
--- a/Simulation/DetSimAna/CMakeLists.txt
+++ b/Simulation/DetSimAna/CMakeLists.txt
@@ -5,12 +5,17 @@ include(${Geant4_USE_FILE})
gaudi_add_module(DetSimAna
SOURCES src/Edm4hepWriterAnaElemTool.cpp
LINK DetSimInterface
+ DetSimSDLib
${DD4hep_COMPONENT_LIBRARIES}
Gaudi::GaudiKernel
EDM4HEP::edm4hep EDM4HEP::edm4hepDict
k4FWCore::k4FWCore
)
+target_include_directories(DetSimAna PUBLIC
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/include
+ $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+
install(TARGETS DetSimAna
EXPORT CEPCSWTargets
RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp
index b05e7994b44cc54d2dbda14039a3700348fcc704..0c687f82701832b673c49e2a744566aa1d76b061 100644
--- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp
+++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp
@@ -1,5 +1,6 @@
#include "Edm4hepWriterAnaElemTool.h"
+#include "G4Step.hh"
#include "G4Event.hh"
#include "G4THitsCollection.hh"
#include "G4EventManager.hh"
@@ -12,13 +13,36 @@
#include "DDG4/Geant4Mapping.h"
#include "DDG4/Geant4HitCollection.h"
#include "DDG4/Geant4Data.h"
-#include "DDG4/Geant4Hits.h"
+#include "DetSimSD/Geant4Hits.h"
DECLARE_COMPONENT(Edm4hepWriterAnaElemTool)
void
Edm4hepWriterAnaElemTool::BeginOfRunAction(const G4Run*) {
+ auto msglevel = msgLevel();
+ if (msglevel == MSG::VERBOSE || msglevel == MSG::DEBUG) {
+ verboseOutput = true;
+ }
+
G4cout << "Begin Run of detector simultion..." << G4endl;
+
+ // access geometry service
+ m_geosvc = service<IGeomSvc>("GeomSvc");
+ if (m_geosvc) {
+ dd4hep::Detector* dd4hep_geo = m_geosvc->lcdd();
+ // try to get the constants
+ R = dd4hep_geo->constant<double>("tracker_region_rmax")/dd4hep::mm*CLHEP::mm;
+ Z = dd4hep_geo->constant<double>("tracker_region_zmax")/dd4hep::mm*CLHEP::mm;
+
+ info() << "Tracker Region "
+ << " R: " << R
+ << " Z: " << Z
+ << endmsg;
+
+ } else {
+ error() << "Failed to find GeomSvc." << endmsg;
+ }
+
}
void
@@ -30,6 +54,15 @@ void
Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) {
msg() << "Event " << anEvent->GetEventID() << endmsg;
+ // event info
+ m_userinfo = nullptr;
+ if (anEvent->GetUserInformation()) {
+ m_userinfo = dynamic_cast<CommonUserEventInfo*>(anEvent->GetUserInformation());
+ if (verboseOutput) {
+ m_userinfo->dumpIdxG4Track2Edm4hep();
+ }
+ }
+
auto mcGenCol = m_mcParGenCol.get();
mcCol = m_mcParCol.createAndPut();
@@ -50,17 +83,18 @@ Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) {
newparticle.setColorFlow (mcGenParticle.getColorFlow());
}
- msg() << "mcCol size: " << mcCol->size() << endmsg;
+ msg() << "mcCol size (original) : " << mcCol->size() << endmsg;
// reset
m_track2primary.clear();
-
+
}
void
Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) {
- // save all data
+ msg() << "mcCol size (after simulation) : " << mcCol->size() << endmsg;
+ // save all data
// create collections.
auto trackercols = m_trackerCol.createAndPut();
auto calorimetercols = m_calorimeterCol.createAndPut();
@@ -257,7 +291,12 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) {
pritrkid = 1;
}
- edm_trk_hit.setMCParticle(mcCol->at(pritrkid-1));
+ if (m_userinfo) {
+ auto idxedm4hep = m_userinfo->idxG4Track2Edm4hep(pritrkid);
+ if (idxedm4hep != -1) {
+ edm_trk_hit.setMCParticle(mcCol->at(idxedm4hep));
+ }
+ }
if (pritrkid != trackID) {
// If the track is a secondary, then the primary track id and current track id is different
@@ -297,7 +336,12 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) {
pritrkid = 1;
}
- edm_calo_contrib.setParticle(mcCol->at(pritrkid-1)); // todo
+ if (m_userinfo) {
+ auto idxedm4hep = m_userinfo->idxG4Track2Edm4hep(pritrkid);
+ if (idxedm4hep != -1) {
+ edm_calo_contrib.setParticle(mcCol->at(idxedm4hep)); // todo
+ }
+ }
edm_calo_hit.addToContributions(edm_calo_contrib);
}
}
@@ -345,17 +389,27 @@ void
Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) {
int curtrkid = track->GetTrackID(); // starts from 1
int curparid = track->GetParentID();
+ int idxedm4hep = -1;
- if (curparid == 0) {
+ if (m_userinfo) {
+ idxedm4hep = m_userinfo->idxG4Track2Edm4hep(curtrkid);
+ }
+
+ if (curparid == 0) { // Primary Track
// select the primary tracks (parentID == 0)
// auto mcCol = m_mcParCol.get();
- if (curtrkid-1>=mcCol->size()) {
+ if (idxedm4hep == -1) {
+ error () << "Failed to get idxedm4hep according to the g4track id " << curtrkid << endmsg;
+ return;
+ }
+
+ if (idxedm4hep>=mcCol->size()) {
error() << "out of range: curtrkid is " << curtrkid
<< " while the MCParticle size is " << mcCol->size() << endmsg;
return;
}
- auto primary_particle = mcCol->at(curtrkid-1);
+ auto primary_particle = mcCol->at(idxedm4hep);
const G4ThreeVector& stop_pos = track->GetPosition();
edm4hep::Vector3d endpoint(stop_pos.x()/CLHEP::mm,
@@ -407,14 +461,14 @@ Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) {
// select the necessary processes
if (creatorProcess==proc_decay) {
- info() << "Creator Process is Decay for secondary particle: "
- << " idx: " << i
- << " trkid: " << sectrk->GetTrackID() // not valid until track
- << " particle: " << secparticle->GetParticleName()
- << " pdg: " << secparticle->GetPDGEncoding()
- << " at position: " << sectrk->GetPosition() //
- << " time: " << sectrk->GetGlobalTime()
- << " momentum: " << sectrk->GetMomentum() //
+ debug() << "Creator Process is Decay for secondary particle: "
+ << " idx: " << i
+ << " trkid: " << sectrk->GetTrackID() // not valid until track
+ << " particle: " << secparticle->GetParticleName()
+ << " pdg: " << secparticle->GetPDGEncoding()
+ << " at position: " << sectrk->GetPosition() //
+ << " time: " << sectrk->GetGlobalTime()
+ << " momentum: " << sectrk->GetMomentum() //
<< endmsg;
is_decay = true;
@@ -445,6 +499,15 @@ Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) {
mcp.addToParents(primary_particle);
primary_particle.addToDaughters(mcp);
+
+ // store the edm4hep obj idx in track info.
+ // using this idx, the MCParticle object could be modified later.
+ auto trackinfo = new CommonUserTrackInfo();
+ trackinfo->setIdxEdm4hep(mcp.getObjectID().index);
+ sectrk->SetUserInformation(trackinfo);
+ debug() << " Appending MCParticle: (id: "
+ << mcp.getObjectID().index << ")"
+ << endmsg;
}
}
}
@@ -462,8 +525,36 @@ Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) {
}
void
-Edm4hepWriterAnaElemTool::UserSteppingAction(const G4Step*) {
-
+Edm4hepWriterAnaElemTool::UserSteppingAction(const G4Step* aStep) {
+ auto aTrack = aStep->GetTrack();
+ // try to get user track info
+ auto trackinfo = dynamic_cast<CommonUserTrackInfo*>(aTrack->GetUserInformation());
+
+ // ========================================================================
+ // Note:
+ // if there is no track info, then do nothing.
+ // ========================================================================
+ if (trackinfo) {
+ // back scattering is defined as following:
+ // - pre point is not in tracker
+ // - post point is in tracker
+ // For details, look at Mokka's source code:
+ // https://llrforge.in2p3.fr/trac/Mokka/browser/trunk/source/Kernel/src/SteppingAction.cc
+ const auto& pre_pos = aStep->GetPreStepPoint()->GetPosition();
+ const auto& post_pos = aStep->GetPostStepPoint()->GetPosition();
+
+ bool is_pre_in_tracker = pre_pos.perp() < R && std::fabs(pre_pos.z()) < Z;
+ bool is_post_in_tracker = post_pos.perp() < R && std::fabs(post_pos.z()) < Z;
+
+ if ((!is_pre_in_tracker) and is_post_in_tracker) {
+ // set back scattering
+ auto idxedm4hep = trackinfo->idxEdm4hep();
+ mcCol->at(idxedm4hep).setBackscatter(true);
+ debug() << " set back scatter for MCParticle "
+ << " (ID: " << idxedm4hep << ")"
+ << endmsg;
+ }
+ }
}
StatusCode
diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h
index 7b992d3076c0ddfa7b7bd9f5d8639274ba0538b1..e4415b17899012530817186c0428581b585bad62 100644
--- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h
+++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h
@@ -6,6 +6,10 @@
#include "GaudiKernel/AlgTool.h"
#include "k4FWCore/DataHandle.h"
#include "DetSimInterface/IAnaElemTool.h"
+#include "DetSimInterface/CommonUserEventInfo.hh"
+#include "DetSimInterface/CommonUserTrackInfo.hh"
+
+#include "DetInterface/IGeomSvc.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimTrackerHitCollection.h"
@@ -42,10 +46,10 @@ public:
private:
// In order to associate MCParticle with contribution, we need to access MC Particle.
// - collection MCParticle: the particles in Generator
- DataHandle<edm4hep::MCParticleCollection> m_mcParGenCol{"MCParticle",
+ DataHandle<edm4hep::MCParticleCollection> m_mcParGenCol{"MCParticleGen",
Gaudi::DataHandle::Writer, this};
// - collection MCParticleG4: the simulated particles in Geant4
- DataHandle<edm4hep::MCParticleCollection> m_mcParCol{"MCParticleG4",
+ DataHandle<edm4hep::MCParticleCollection> m_mcParCol{"MCParticle",
Gaudi::DataHandle::Writer, this};
edm4hep::MCParticleCollection* mcCol;
@@ -140,6 +144,14 @@ private:
std::map<int, int> m_track2primary;
+ CommonUserEventInfo* m_userinfo = nullptr;
+
+ // get the limitation of R/Z in tracker
+ SmartIF<IGeomSvc> m_geosvc;
+ double R = 0;
+ double Z = 0;
+
+ bool verboseOutput = false;
};
#endif
diff --git a/Simulation/DetSimAna/src/ExampleAnaElemTool.cpp b/Simulation/DetSimAna/src/ExampleAnaElemTool.cpp
index d912414d4108a8f45efdb394d26ddcfab5fc0d21..1e11da48b5c8872ff1553c95025597a2a7db251c 100644
--- a/Simulation/DetSimAna/src/ExampleAnaElemTool.cpp
+++ b/Simulation/DetSimAna/src/ExampleAnaElemTool.cpp
@@ -9,7 +9,7 @@
#include "DDG4/Geant4Mapping.h"
#include "DDG4/Geant4HitCollection.h"
#include "DDG4/Geant4Data.h"
-#include "DDG4/Geant4Hits.h"
+#include "DetSimSD/Geant4Hits.h"
DECLARE_COMPONENT(ExampleAnaElemTool)
diff --git a/Simulation/DetSimInterface/CMakeLists.txt b/Simulation/DetSimInterface/CMakeLists.txt
index ed99eead6c48943bbe639731e2705a8e70423d31..21963df1a8e0fd3289849181995d25908f69641e 100644
--- a/Simulation/DetSimInterface/CMakeLists.txt
+++ b/Simulation/DetSimInterface/CMakeLists.txt
@@ -4,7 +4,10 @@
gaudi_add_library(DetSimInterface
SOURCES src/IDetSimSvc.cpp
+ src/CommonUserEventInfo.cc
+ src/CommonUserTrackInfo.cc
LINK Gaudi::GaudiKernel
+ ${Geant4_LIBRARIES}
)
install(TARGETS DetSimInterface
diff --git a/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh b/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3319598a270d0d4cb38607d87b6c4fd3b2a85ad4
--- /dev/null
+++ b/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh
@@ -0,0 +1,39 @@
+#ifndef CommonUserEventInfo_hh
+#define CommonUserEventInfo_hh
+
+/*
+ * Description:
+ * This class is a part of simulation framework to allow users to extend the G4Event.
+ *
+ * For example, when G4 converts the EDM4hep/G4 primary vertex/particle to G4 track,
+ * the relationship between the EDM4hep track and G4 track is missing.
+ * So a map is used as a bookkeeping.
+ *
+ * Author:
+ * Tao Lin <lintao AT ihep.ac.cn>
+ */
+
+#include "G4VUserEventInformation.hh"
+#include <map>
+
+class CommonUserEventInfo: public G4VUserEventInformation {
+public:
+
+ CommonUserEventInfo();
+ virtual ~CommonUserEventInfo();
+
+public:
+ virtual void Print() const;
+
+ // set the relationship between idx in geant4 and idx in mc particle collection.
+ // idxG4: G4 track ID (starts from 1)
+ // idxEdm4hep: index in MC Particle collection (starts from 0)
+ bool setIdxG4Track2Edm4hep(int idxG4, int idxEdm4hep);
+ int idxG4Track2Edm4hep(int idxG4) const;
+ void dumpIdxG4Track2Edm4hep() const;
+
+private:
+ std::map<int, int> m_g4track_to_edm4hep;
+};
+
+#endif
diff --git a/Simulation/DetSimInterface/include/DetSimInterface/CommonUserTrackInfo.hh b/Simulation/DetSimInterface/include/DetSimInterface/CommonUserTrackInfo.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ce63e602d26267dd87465406d454c7d558eb83ad
--- /dev/null
+++ b/Simulation/DetSimInterface/include/DetSimInterface/CommonUserTrackInfo.hh
@@ -0,0 +1,35 @@
+#ifndef CommonUserTrackInfo_hh
+#define CommonUserTrackInfo_hh
+
+/* Description:
+ * This class is a part of simulation framework to extend the G4Track.
+ *
+ * Some secondaries are created due to decay. However, their G4 Track IDs are
+ * not valid until the tracks are tracking by geant4. In order to associate
+ * these tracks and their edm4hep MC particle, we use the track information
+ * to record the extra track information.
+ *
+ * Author:
+ * Tao Lin <lintao AT ihep.ac.cn>
+ */
+
+#include "G4VUserTrackInformation.hh"
+
+class CommonUserTrackInfo: public G4VUserTrackInformation {
+public:
+ CommonUserTrackInfo();
+ ~CommonUserTrackInfo();
+
+public:
+
+ virtual void Print() const;
+
+ // get the idx in the EDM4hep MC particle collection
+ bool setIdxEdm4hep(int idxEdm4hep);
+ int idxEdm4hep() const;
+
+private:
+ int m_idxEdm4hep = -1;
+};
+
+#endif
diff --git a/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h b/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h
index 8c804500a7b4322e1998a04a7a8a43f46ba3f499..4da01d532a69f4cba7078d82fb161fcc34a1634c 100644
--- a/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h
+++ b/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h
@@ -28,6 +28,7 @@ public:
virtual double dedx(const G4Step* aStep) = 0;
virtual double dedx(const edm4hep::MCParticle& mc) = 0;
virtual double dndx(double betagamma) = 0;
+ virtual void endOfEvent() {}
};
diff --git a/Simulation/DetSimInterface/src/CommonUserEventInfo.cc b/Simulation/DetSimInterface/src/CommonUserEventInfo.cc
new file mode 100644
index 0000000000000000000000000000000000000000..7992c69c897967e728040b3602a472663f604c6f
--- /dev/null
+++ b/Simulation/DetSimInterface/src/CommonUserEventInfo.cc
@@ -0,0 +1,53 @@
+#include "DetSimInterface/CommonUserEventInfo.hh"
+#include <iostream>
+
+CommonUserEventInfo::CommonUserEventInfo() {
+
+}
+
+CommonUserEventInfo::~CommonUserEventInfo() {
+
+}
+
+void
+CommonUserEventInfo::Print() const {
+
+}
+
+bool
+CommonUserEventInfo::setIdxG4Track2Edm4hep(int idxG4, int idxEdm4hep) {
+ auto it = m_g4track_to_edm4hep.find(idxG4);
+
+ // if already exists, return false
+ if (it != m_g4track_to_edm4hep.end()) {
+ return false;
+ }
+
+ m_g4track_to_edm4hep[idxG4] = idxEdm4hep;
+
+ return true;
+}
+
+int
+CommonUserEventInfo::idxG4Track2Edm4hep(int idxG4) const {
+ int ret = -1;
+
+ auto it = m_g4track_to_edm4hep.find(idxG4);
+
+ // if found
+ if (it != m_g4track_to_edm4hep.end()) {
+ ret = it->second;
+ }
+
+ return ret;
+}
+
+void
+CommonUserEventInfo::dumpIdxG4Track2Edm4hep() const {
+ std::cout << "---- Dumping IdxG4Track2Edm4hep: "
+ << " total size: " << m_g4track_to_edm4hep.size()
+ << std::endl;
+ for (auto [idxG4, idxE4]: m_g4track_to_edm4hep) {
+ std::cout << " - " << idxG4 << " -> " << idxE4 << std::endl;
+ }
+}
diff --git a/Simulation/DetSimInterface/src/CommonUserTrackInfo.cc b/Simulation/DetSimInterface/src/CommonUserTrackInfo.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9a923a35a22e7fafe496c827e066e591864f3c50
--- /dev/null
+++ b/Simulation/DetSimInterface/src/CommonUserTrackInfo.cc
@@ -0,0 +1,22 @@
+#include "DetSimInterface/CommonUserTrackInfo.hh"
+#include <iostream>
+
+CommonUserTrackInfo::CommonUserTrackInfo() {
+
+}
+
+CommonUserTrackInfo::~CommonUserTrackInfo() {
+
+}
+
+void CommonUserTrackInfo::Print() const {
+
+}
+
+bool CommonUserTrackInfo::setIdxEdm4hep(int idxEdm4hep) {
+ m_idxEdm4hep = idxEdm4hep;
+}
+
+int CommonUserTrackInfo::idxEdm4hep() const {
+ return m_idxEdm4hep;
+}
diff --git a/Simulation/DetSimSD/CMakeLists.txt b/Simulation/DetSimSD/CMakeLists.txt
index fa8ba9ef3e13d27da31d8c25763278ba795ef47c..f3a9618ac625fce03d1c277bcef6be28e355628e 100644
--- a/Simulation/DetSimSD/CMakeLists.txt
+++ b/Simulation/DetSimSD/CMakeLists.txt
@@ -2,20 +2,34 @@
find_package(Geant4 REQUIRED ui_all vis_all)
include(${Geant4_USE_FILE})
+gaudi_add_library(DetSimSDLib
+ SOURCES src/Geant4Hits.cpp
+ src/DDG4SensitiveDetector.cpp
+ src/CaloSensitiveDetector.cpp
+ src/DriftChamberSensitiveDetector.cpp
+ src/TimeProjectionChamberSensitiveDetector.cpp
+ src/GenericTrackerSensitiveDetector.cpp
+ src/TrackerCombineSensitiveDetector.cpp
+ LINK DetSimInterface
+ DetInterface
+ ${DD4hep_COMPONENT_LIBRARIES}
+)
+target_include_directories(DetSimSDLib PUBLIC
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/include
+ $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+install(TARGETS DetSimSDLib
+ EXPORT CEPCSWTargets
+ RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+ LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+ COMPONENT dev)
+
gaudi_add_module(DetSimSD
SOURCES src/CalorimeterSensDetTool.cpp
- src/DDG4SensitiveDetector.cpp
- src/CaloSensitiveDetector.cpp
-
src/DriftChamberSensDetTool.cpp
- src/DriftChamberSensitiveDetector.cpp
-
src/TimeProjectionChamberSensDetTool.cpp
- src/TimeProjectionChamberSensitiveDetector.cpp
-
src/GenericTrackerSensDetTool.cpp
- src/GenericTrackerSensitiveDetector.cpp
LINK DetSimInterface
+ DetSimSDLib
DetInterface
${DD4hep_COMPONENT_LIBRARIES}
Gaudi::GaudiKernel
diff --git a/Simulation/DetSimSD/include/DetSimSD/DDG4SensitiveDetector.h b/Simulation/DetSimSD/include/DetSimSD/DDG4SensitiveDetector.h
index 687771e344b9f12ae553e56888b1222288fefbaf..c4a29eb3983d0e1f240059eb1c767b8fc0773c15 100644
--- a/Simulation/DetSimSD/include/DetSimSD/DDG4SensitiveDetector.h
+++ b/Simulation/DetSimSD/include/DetSimSD/DDG4SensitiveDetector.h
@@ -15,7 +15,7 @@
#include "DD4hep/Detector.h"
-#include "DDG4/Geant4Hits.h"
+#include "DetSimSD/Geant4Hits.h"
#include "G4Step.hh"
#include "G4HCofThisEvent.hh"
diff --git a/Simulation/DetSimSD/include/DetSimSD/Geant4Hits.h b/Simulation/DetSimSD/include/DetSimSD/Geant4Hits.h
new file mode 100644
index 0000000000000000000000000000000000000000..60fc18d4792876eac5b8ebe6c558c74447a4f6c3
--- /dev/null
+++ b/Simulation/DetSimSD/include/DetSimSD/Geant4Hits.h
@@ -0,0 +1,201 @@
+//==========================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------------
+// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+// All rights reserved.
+//
+// For the licensing terms see $DD4hepINSTALL/LICENSE.
+// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+//
+// Author : M.Frank
+//
+//==========================================================================
+
+#ifndef DDG4_GEANT4HITS_H
+#define DDG4_GEANT4HITS_H
+
+// Framework include files
+#include "DD4hep/Objects.h"
+#include "DDG4/Geant4StepHandler.h"
+
+// Geant4 include files
+#include "G4VHit.hh"
+#include "G4Step.hh"
+#include "G4StepPoint.hh"
+
+/// Namespace for the AIDA detector description toolkit
+namespace dd4hep {
+
+ /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
+ namespace sim {
+
+ // Forward declarations;
+ class Geant4Hit;
+ class Geant4TrackerHit;
+ class Geant4CalorimeterHit;
+
+ /// Deprecated: Base class for hit comparisons.
+ /** @class HitCompare Geant4Hits.h DDG4/Geant4Hits.h
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ template <class HIT> class HitCompare {
+ public:
+ /// Default destructor
+ virtual ~HitCompare() {}
+ /// Comparison function
+ virtual bool operator()(const HIT* h) const = 0;
+ };
+
+ /// Deprecated: Seek the hits of an arbitrary collection for the same position.
+ /** @class HitPositionCompare Geant4Hits.h DDG4/Geant4Hits.h
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ template <class HIT> struct HitPositionCompare: public HitCompare<HIT> {
+ /// Reference to the hit position
+ const Position& pos;
+ /// Constructor
+ HitPositionCompare(const Position& p) : pos(p) {}
+ /// Default destructor
+ virtual ~HitPositionCompare() {}
+ /// Comparison function
+ virtual bool operator()(const HIT* h) const {
+ return pos == h->position;
+ }
+ };
+
+ /// Deprecated: basic geant4 hit class for deprecated sensitive detectors
+ /** @class Geant4Hit Geant4Hits.h DDG4/Geant4Hits.h
+ *
+ * Geant4 hit base class. Here only the basic
+ * quantites are stored such as the energy deposit and
+ * the time stamp.
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4Hit: public G4VHit {
+ public:
+
+ // cellID
+ unsigned long cellID = 0;
+
+ /// Deprecated!!!
+ struct MonteCarloContrib {
+ /// Geant 4 Track identifier
+ int trackID = -1;
+ /// Particle ID from the PDG table
+ int pdgID = -1;
+ /// Total energy deposit in this hit
+ double deposit = 0.0;
+ /// Timestamp when this energy was deposited
+ double time = 0.0;
+ /// Default constructor
+ MonteCarloContrib() = default;
+ /// Copy constructor
+ MonteCarloContrib(const MonteCarloContrib& c) = default;
+ /// Initializing constructor
+ MonteCarloContrib(int track_id, double dep, double time_stamp)
+ : trackID(track_id), pdgID(-1), deposit(dep), time(time_stamp) {}
+ /// Initializing constructor
+ MonteCarloContrib(int track_id, int pdg, double dep, double time_stamp)
+ : trackID(track_id), pdgID(pdg), deposit(dep), time(time_stamp) {}
+ /// Assignment operator
+ MonteCarloContrib& operator=(const MonteCarloContrib& c) = default;
+ /// Clear data
+ void clear() {
+ time = deposit = 0.0;
+ pdgID = trackID = -1;
+ }
+ };
+ typedef MonteCarloContrib Contribution;
+ typedef std::vector<MonteCarloContrib> Contributions;
+
+ public:
+ /// Standard constructor
+ Geant4Hit() = default;
+ /// Default destructor
+ virtual ~Geant4Hit() { }
+ /// Check if the Geant4 track is a Geantino
+ static bool isGeantino(G4Track* track);
+ /// Extract the MC contribution for a given hit from the step information
+ static Contribution extractContribution(G4Step* step);
+ };
+
+ /// Deprecated: Geant4 tracker hit class for deprecated sensitive detectors
+ /** @class Geant4TrackerHit Geant4Hits.h DDG4/Geant4Hits.h
+ *
+ * Geant4 tracker hit class. Tracker hits contain the momentum
+ * direction as well as the hit position.
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4TrackerHit: public Geant4Hit {
+ public:
+ /// Hit position
+ Position position;
+ /// Hit direction
+ Direction momentum;
+ /// Length of the track segment contributing to this hit
+ double length;
+ /// Monte Carlo / Geant4 information
+ Contribution truth;
+ /// Energy deposit of the hit
+ double energyDeposit;
+
+ public:
+ /// Default constructor
+ Geant4TrackerHit();
+ /// Standard initializing constructor
+ Geant4TrackerHit(int track_id, int pdg_id, double deposit, double time_stamp);
+ /// Default destructor
+ virtual ~Geant4TrackerHit() {}
+ /// Clear hit content
+ Geant4TrackerHit& clear();
+ /// Store Geant4 point and step information into tracker hit structure.
+ Geant4TrackerHit& storePoint(G4Step* step, G4StepPoint* point);
+
+ /// Assignment operator
+ Geant4TrackerHit& operator=(const Geant4TrackerHit& c);
+ /// Geant4 required object allocator
+ void *operator new(size_t);
+ /// Geat4 required object destroyer
+ void operator delete(void *ptr);
+ };
+
+ /// Deprecated: Geant4 calorimeter hit class for deprecated sensitive detectors
+ /** @class Geant4CalorimeterHit Geant4Hits.h DDG4/Geant4Hits.h
+ *
+ * Geant4 calorimeter hit class. Calorimeter hits contain the momentum
+ * direction as well as the hit position.
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4CalorimeterHit: public Geant4Hit {
+ public:
+ /// Hit position
+ Position position;
+ /// Hit contributions by individual particles
+ Contributions truth;
+ /// Total energy deposit
+ double energyDeposit;
+ public:
+ /// Standard constructor
+ Geant4CalorimeterHit(const Position& cell_pos);
+ /// Default destructor
+ virtual ~Geant4CalorimeterHit() { }
+ /// Geant4 required object allocator
+ void *operator new(size_t);
+ /// Geat4 required object destroyer
+ void operator delete(void *ptr);
+ };
+
+ } // End namespace sim
+} // End namespace dd4hep
+
+#endif // DDG4_GEANT4HITS_H
diff --git a/Simulation/DetSimSD/src/DDG4SensitiveDetector.cpp b/Simulation/DetSimSD/src/DDG4SensitiveDetector.cpp
index a1698ee3cd3cbe43b6ae856112a2edb2f9f20ed0..3d16e02be94bc3ecfe30e161e81b80bb683892f7 100644
--- a/Simulation/DetSimSD/src/DDG4SensitiveDetector.cpp
+++ b/Simulation/DetSimSD/src/DDG4SensitiveDetector.cpp
@@ -1,7 +1,6 @@
#include "DetSimSD/DDG4SensitiveDetector.h"
#include "DDG4/Geant4Converter.h"
-#include "DDG4/Geant4Hits.h"
#include "DD4hep/Segmentations.h"
#include "DD4hep/Printout.h"
diff --git a/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp b/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp
index 43dfa17828453c1545318eea45c0aec7e271976e..58c49d6c62e5d3646348995b52e815a2bc42b31f 100644
--- a/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp
+++ b/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp
@@ -5,7 +5,7 @@
#include "DD4hep/Detector.h"
#include "DriftChamberSensitiveDetector.h"
-
+#include "TrackerCombineSensitiveDetector.h"
DECLARE_COMPONENT(DriftChamberSensDetTool)
StatusCode DriftChamberSensDetTool::initialize() {
@@ -22,7 +22,6 @@ StatusCode DriftChamberSensDetTool::initialize() {
error() << "Failed to find dedx simtoo." << endmsg;
return StatusCode::FAILURE;
}
-
return sc;
}
@@ -39,12 +38,19 @@ DriftChamberSensDetTool::createSD(const std::string& name) {
G4VSensitiveDetector* sd = nullptr;
if (name == "DriftChamber") {
- DriftChamberSensitiveDetector* dcsd = new DriftChamberSensitiveDetector(name, *dd4hep_geo);
- dcsd->setDedxSimTool(m_dedx_simtool);
-
- sd = dcsd;
+ auto sens = dd4hep_geo->sensitiveDetector(name);
+ if(!sens.combineHits()){
+ DriftChamberSensitiveDetector* dcsd = new DriftChamberSensitiveDetector(name, *dd4hep_geo);
+ dcsd->setDedxSimTool(m_dedx_simtool);
+ sd = dcsd;
+ }
+ else{
+ sd = new TrackerCombineSensitiveDetector(name, *dd4hep_geo);
+ }
+ }
+ else{
+ warning() << "The SD " << name << " want to use SD for DriftChamber" << endmsg;
}
-
return sd;
}
diff --git a/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp
index c6ff3db08209c3307b8fd2c954ee9c97f5c015cc..44bf01d1016ae3578ef192f609a3163b6ea6b483 100644
--- a/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp
+++ b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp
@@ -71,5 +71,5 @@ DriftChamberSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) {
void
DriftChamberSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE) {
-
+ m_dedx_simtool->endOfEvent();
}
diff --git a/Simulation/DetSimSD/src/Geant4Hits.cpp b/Simulation/DetSimSD/src/Geant4Hits.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0358bf019d8222bc8a92ab02bc4fb2022e60b1cf
--- /dev/null
+++ b/Simulation/DetSimSD/src/Geant4Hits.cpp
@@ -0,0 +1,135 @@
+//==========================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------------
+// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+// All rights reserved.
+//
+// For the licensing terms see $DD4hepINSTALL/LICENSE.
+// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+//
+// Author : M.Frank
+//
+//==========================================================================
+
+// Framework include files
+#include "DetSimSD/Geant4Hits.h"
+
+// Geant4 include files
+#include "G4Allocator.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4ChargedGeantino.hh"
+#include "G4OpticalPhoton.hh"
+#include "G4Geantino.hh"
+
+// C/C++ include files
+#include <iostream>
+
+using namespace std;
+using namespace dd4hep::sim;
+
+G4ThreadLocal G4Allocator<Geant4TrackerHit>* TrackerHitAllocator = 0;
+G4ThreadLocal G4Allocator<Geant4CalorimeterHit>* CalorimeterHitAllocator = 0;
+
+
+/// Check if the Geant4 track is a Geantino
+bool Geant4Hit::isGeantino(G4Track* track) {
+ if (track) {
+ G4ParticleDefinition* def = track->GetDefinition();
+ if (def == G4ChargedGeantino::Definition())
+ return true;
+ if (def == G4Geantino::Definition()) {
+ return true;
+ }
+ }
+ return false;
+}
+
+Geant4Hit::Contribution Geant4Hit::extractContribution(G4Step* step) {
+ G4Track* trk = step->GetTrack();
+ double energy_deposit =
+ (trk->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) ?
+ trk->GetTotalEnergy() : step->GetTotalEnergyDeposit();
+ Contribution contrib(trk->GetTrackID(), trk->GetDefinition()->GetPDGEncoding(), energy_deposit, trk->GetGlobalTime());
+ return contrib;
+}
+
+/// Default constructor
+Geant4TrackerHit::Geant4TrackerHit()
+ : Geant4Hit(), position(), momentum(), length(0.0), truth(), energyDeposit(0.0) {
+}
+
+/// Standard initializing constructor
+Geant4TrackerHit::Geant4TrackerHit(int track_id, int pdg_id, double deposit, double time_stamp)
+ : Geant4Hit(), position(), momentum(), length(0.0), truth(track_id, pdg_id, deposit, time_stamp), energyDeposit(deposit) {
+}
+
+/// Assignment operator
+Geant4TrackerHit& Geant4TrackerHit::operator=(const Geant4TrackerHit& c) {
+ if ( this != &c ) {
+ position = c.position;
+ momentum = c.momentum;
+ length = c.length;
+ truth = c.truth;
+ energyDeposit = c.energyDeposit;
+ }
+ return *this;
+}
+
+/// Clear hit content
+Geant4TrackerHit& Geant4TrackerHit::clear() {
+ position.SetXYZ(0, 0, 0);
+ momentum.SetXYZ(0, 0, 0);
+ length = 0.0;
+ truth.clear();
+ energyDeposit = 0.0;
+ return *this;
+}
+
+/// Store Geant4 point and step information into tracker hit structure.
+Geant4TrackerHit& Geant4TrackerHit::storePoint(G4Step* step, G4StepPoint* pnt) {
+ G4Track* trk = step->GetTrack();
+ G4ThreeVector pos = pnt->GetPosition();
+ G4ThreeVector mom = pnt->GetMomentum();
+
+ truth.trackID = trk->GetTrackID();
+ truth.pdgID = trk->GetDefinition()->GetPDGEncoding();
+ truth.deposit = step->GetTotalEnergyDeposit();
+ truth.time = trk->GetGlobalTime();
+ energyDeposit = step->GetTotalEnergyDeposit();
+ position.SetXYZ(pos.x(), pos.y(), pos.z());
+ momentum.SetXYZ(mom.x(), mom.y(), mom.z());
+ length = 0;
+ return *this;
+}
+
+/// Geant4 required object allocator
+void* Geant4TrackerHit::operator new(size_t) {
+ if ( TrackerHitAllocator )
+ return TrackerHitAllocator->MallocSingle();
+ TrackerHitAllocator = new G4Allocator<Geant4TrackerHit>;
+ return TrackerHitAllocator->MallocSingle();
+}
+
+/// Geat4 required object destroyer
+void Geant4TrackerHit::operator delete(void *p) {
+ TrackerHitAllocator->FreeSingle((Geant4TrackerHit*) p);
+}
+
+/// Standard constructor
+Geant4CalorimeterHit::Geant4CalorimeterHit(const Position& pos)
+ : Geant4Hit(), position(pos), truth(), energyDeposit(0) {
+}
+
+/// Geant4 required object allocator
+void* Geant4CalorimeterHit::operator new(size_t) {
+ if ( CalorimeterHitAllocator )
+ return CalorimeterHitAllocator->MallocSingle();
+ CalorimeterHitAllocator = new G4Allocator<Geant4CalorimeterHit>;
+ return CalorimeterHitAllocator->MallocSingle();
+}
+
+/// Geat4 required object destroyer
+void Geant4CalorimeterHit::operator delete(void *p) {
+ CalorimeterHitAllocator->FreeSingle((Geant4CalorimeterHit*) p);
+}
+
diff --git a/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp b/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp
index 480f509aa57c563e90c032aac64abcaca4957a31..edb275a8d76e4d0710368611e916f4c026837e12 100644
--- a/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp
+++ b/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp
@@ -5,6 +5,7 @@
#include "DD4hep/Detector.h"
#include "GenericTrackerSensitiveDetector.h"
+#include "TrackerCombineSensitiveDetector.h"
#include "CLHEP/Units/SystemOfUnits.h"
@@ -33,8 +34,14 @@ G4VSensitiveDetector* GenericTrackerSensDetTool::createSD(const std::string& nam
debug() << "createSD for " << name << endmsg;
dd4hep::Detector* dd4hep_geo = m_geosvc->lcdd();
-
- G4VSensitiveDetector* sd = new GenericTrackerSensitiveDetector(name, *dd4hep_geo);
+ auto sens = dd4hep_geo->sensitiveDetector(name);
+ G4VSensitiveDetector* sd = nullptr;
+ if(sens.combineHits()){
+ sd = new TrackerCombineSensitiveDetector(name, *dd4hep_geo);
+ }
+ else{
+ sd = new GenericTrackerSensitiveDetector(name, *dd4hep_geo);
+ }
return sd;
}
diff --git a/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp
index d250079ef0a8aa5f15e50fb801b1de90c46d3a17..5e9ba1a6af4cf2dbe994a1bf7d686f3bfcbf19f7 100644
--- a/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp
+++ b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp
@@ -23,10 +23,11 @@ void GenericTrackerSensitiveDetector::Initialize(G4HCofThisEvent* HCE){
}
G4bool GenericTrackerSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*){
-
G4TouchableHandle touchPost = step->GetPostStepPoint()->GetTouchableHandle();
G4TouchableHandle touchPre = step->GetPreStepPoint()->GetTouchableHandle();
dd4hep::sim::Geant4StepHandler h(step);
+ if (fabs(h.trackDef()->GetPDGCharge()) < 0.01) return true;
+
dd4hep::Position prePos = h.prePos();
dd4hep::Position postPos = h.postPos();
dd4hep::Position direction = postPos - prePos;
diff --git a/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp b/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..591c3cf3273a872f8533a8d11f3301fe60a408eb
--- /dev/null
+++ b/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp
@@ -0,0 +1,58 @@
+#include "TrackerCombineSensitiveDetector.h"
+
+#include "G4Step.hh"
+#include "G4VProcess.hh"
+#include "G4SDManager.hh"
+#include "DD4hep/DD4hepUnits.h"
+
+TrackerCombineSensitiveDetector::TrackerCombineSensitiveDetector(const std::string& name,
+ dd4hep::Detector& description)
+ : DDG4SensitiveDetector(name, description),
+ m_hc(nullptr){
+ G4String CollName = m_sensitive.hitsCollection();
+ collectionName.insert(CollName);
+}
+
+void TrackerCombineSensitiveDetector::Initialize(G4HCofThisEvent* HCE){
+ userData.e_cut = m_sensitive.energyCutoff();
+
+ m_hc = new HitCollection(GetName(), collectionName[0]);
+ int HCID = G4SDManager::GetSDMpointer()->GetCollectionID(m_hc);
+ HCE->AddHitsCollection( HCID, m_hc );
+}
+
+G4bool TrackerCombineSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*){
+ dd4hep::sim::Geant4StepHandler h(step);
+ if (fabs(h.trackDef()->GetPDGCharge()) < 0.01) return true;
+
+ bool return_code = false;
+ if ( userData.current == -1 ) userData.start(getCellID(step), step, h.pre);
+ else if ( !userData.track || userData.current != h.track->GetTrackID() ) {
+ return_code = userData.extractHit(m_hc) != 0;
+ userData.start(getCellID(step), step, h.pre);
+ }
+
+ // ....update .....
+ userData.update(step);
+
+ void *prePV = h.volume(h.pre), *postPV = h.volume(h.post);
+ if ( prePV != postPV ) {
+ return_code = userData.extractHit(m_hc) != 0;
+ void* postSD = h.sd(h.post);
+ if ( 0 != postSD ) {
+ void* preSD = h.sd(h.pre);
+ if ( preSD == postSD ) {
+ // fucd: getCellID(step) for preVolume not postVolume, so should start at next step
+ //userData.start(getCellID(step), step, h.post);
+ }
+ }
+ }
+ else if ( userData.track->GetTrackStatus() == fStopAndKill ) {
+ return_code = userData.extractHit(m_hc) != 0;
+ }
+ return return_code;
+}
+
+void TrackerCombineSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE){
+ userData.clear();
+}
diff --git a/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.h b/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.h
new file mode 100644
index 0000000000000000000000000000000000000000..cf0cb844eb64a076d8c39477c94a6fd62f677402
--- /dev/null
+++ b/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.h
@@ -0,0 +1,95 @@
+// *********************************************************
+//
+// $Id: TrackerCombineSensitiveDetector.hh,v 1.0 2022/03/27
+
+#ifndef TrackerCombineSensitiveDetector_h
+#define TrackerCombineSensitiveDetector_h
+
+#include "DetSimSD/DDG4SensitiveDetector.h"
+#include "DDG4/Defs.h"
+
+namespace dd4hep {
+ namespace sim {
+ struct TrackerCombine {
+ Geant4TrackerHit pre;
+ Geant4TrackerHit post;
+ G4Track* track;
+ double e_cut;
+ int current;
+ long long int cellID;
+ TrackerCombine() : pre(), post(), track(0), e_cut(0.0), current(-1), cellID(0) {}
+ void start(long long int cell, G4Step* step, G4StepPoint* point) {
+ pre.storePoint(step,point);
+ current = pre.truth.trackID;
+ track = step->GetTrack();
+ cellID = cell;
+ post = pre;
+ //std::cout << "start: " << cellID << " " << pre.position << " current = " << current << std::endl;
+ }
+ void update(G4Step* step) {
+ post.storePoint(step,step->GetPostStepPoint());
+ pre.truth.deposit += post.truth.deposit;
+ //std::cout << "update: " << cellID << " " << post.position << " current = " << current << std::endl;
+ }
+ void clear() {
+ pre.truth.clear();
+ current = -1;
+ track = 0;
+ }
+ Geant4TrackerHit* extractHit(DDG4SensitiveDetector::HitCollection* c) {
+ //std::cout << "create Geant4TrackerHit: " << cellID << " current = " << current << " track = " << track << " de = " << pre.truth.deposit << std::endl;
+ if ( current == -1 || !track ) {
+ return 0;
+ }
+ else if ( pre.truth.deposit <= e_cut && !Geant4Hit::isGeantino(track) ) {
+ clear();
+ return 0;
+ }
+ double rho1 = pre.position.Rho();
+ double rho2 = post.position.Rho();
+ double rho = 0.5*(rho1+rho2);
+ Position pos = 0.5 * (pre.position + post.position);
+ double z = pos.z();
+ double r = sqrt(rho*rho+z*z);
+ Position path = post.position - pre.position;
+ double angle_O_pre_post = acos(-pre.position.Unit().Dot(path.Unit()));
+ double angle_O_post_pre = acos(post.position.Unit().Dot(path.Unit()));
+ double angle_O_P_pre = asin(pre.position.R()*sin(angle_O_pre_post)/r);
+ if(angle_O_pre_post>dd4hep::halfpi||angle_O_post_pre>dd4hep::halfpi){
+ bool backward = angle_O_post_pre>angle_O_pre_post;
+ double angle_O_P_pre = backward ? dd4hep::pi - asin(pre.position.R()*sin(angle_O_pre_post)/r) : asin(pre.position.R()*sin(angle_O_pre_post)/r);
+ double pre2P = r/sin(angle_O_pre_post)*sin(angle_O_pre_post+angle_O_P_pre);
+ pos = pre.position + pre2P*path.Unit();
+ }
+ Momentum mom = 0.5 * (pre.momentum + post.momentum);
+ Geant4TrackerHit* hit = new Geant4TrackerHit(pre.truth.trackID,
+ pre.truth.pdgID,
+ pre.truth.deposit,
+ pre.truth.time);
+ hit->cellID = cellID;
+ hit->position = pos;
+ hit->momentum = mom;
+ hit->length = path.R();;
+ clear();
+ c->insert(hit);
+ return hit;
+ }
+ };
+ }
+}
+
+class TrackerCombineSensitiveDetector: public DDG4SensitiveDetector {
+ public:
+ TrackerCombineSensitiveDetector(const std::string& name, dd4hep::Detector& description);
+
+ void Initialize(G4HCofThisEvent* HCE);
+ G4bool ProcessHits(G4Step* step, G4TouchableHistory* history);
+ void EndOfEvent(G4HCofThisEvent* HCE);
+
+ protected:
+ HitCollection* m_hc = nullptr;
+
+ private:
+ dd4hep::sim::TrackerCombine userData;
+};
+#endif
diff --git a/build.sh b/build.sh
index 337287fe081bdbf364c5fa2b856dc0a4534aa9df..5c7393d169fd5b35240df699e7c1314b5cc783a3 100755
--- a/build.sh
+++ b/build.sh
@@ -83,7 +83,12 @@ function run-cmake() {
}
function run-make() {
- cmake --build .
+ local njobs=-j$(nproc)
+ cmake --build . $njobs
+}
+
+function run-install() {
+ cmake --install .
}
##############################################################################
@@ -91,8 +96,8 @@ function run-make() {
##############################################################################
# The current default platform
-lcg_platform=x86_64-centos7-gcc8-opt
-lcg_version=101.0.1
+lcg_platform=x86_64-centos7-gcc11-opt
+lcg_version=103.0.2
bldtool=${CEPCSW_BLDTOOL} # make, ninja # set in env var
@@ -103,3 +108,5 @@ check-working-builddir || exit -1
run-cmake || exit -1
run-make || exit -1
+
+run-install || exit -1
diff --git a/setup.sh b/setup.sh
index a75d5f39d2d3fca4d0501c970c74648b69229014..3e09a68fff76f3d16300423d10bd8aba3d713317 100644
--- a/setup.sh
+++ b/setup.sh
@@ -9,6 +9,8 @@
# Author: Tao Lin <lintao@ihep.ac.cn>
##############################################################################
+THISSCRITDIR=$(dirname $(readlink -e "${BASH_SOURCE[0]}" 2>/dev/null) 2>/dev/null) # Darwin readlink doesnt accept -e
+
function info:() {
echo "INFO: $*" 1>&2
}
@@ -41,6 +43,25 @@ function setup-external() {
}
+function setup-install-area() {
+ local installarea=$THISSCRITDIR/InstallArea
+ if [ ! -d "$installarea" ]; then
+ return
+ fi
+
+ export PATH=$installarea/bin:$PATH
+ export LD_LIBRARY_PATH=$installarea/lib:$LD_LIBRARY_PATH
+ export PYTHONPATH=$installarea/lib:$PYTHONPATH
+ export PYTHONPATH=$installarea/python:$PYTHONPATH
+ export ROOT_INCLUDE_PATH=$installarea/include:$ROOT_INCLUDE_PATH
+
+ local extrasetupscript=$installarea/setup.sh
+ if [ -f "$extrasetupscript" ]; then
+ source $extrasetupscript
+ fi
+
+ info: "Setup CEPCSW: $installarea"
+}
##############################################################################
# Parse the command line options
@@ -49,8 +70,9 @@ function setup-external() {
# CEPCSW_LCG_VERSION=${1}; shift
if [ -z "$CEPCSW_LCG_VERSION" ]; then
- CEPCSW_LCG_VERSION=101.0.1
+ CEPCSW_LCG_VERSION=103.0.2
fi
export CEPCSW_LCG_VERSION
setup-external
+setup-install-area