diff --git a/Detector/DetDriftChamber/compact/det.xml b/Detector/DetDriftChamber/compact/det.xml
index 7ee0d9ae62e7700f6cffcb4770329e0383344907..a070d543de61ff6e44f8952b24be492cea927c9f 100644
--- a/Detector/DetDriftChamber/compact/det.xml
+++ b/Detector/DetDriftChamber/compact/det.xml
@@ -16,7 +16,7 @@
</includes>
<define>
- <constant name="world_size" value="30*m"/>
+ <constant name="world_size" value="2226*mm"/>
<constant name="world_x" value="world_size"/>
<constant name="world_y" value="world_size"/>
<constant name="world_z" value="world_size"/>
@@ -65,7 +65,7 @@
</regions>
<detectors>
- <detector id="1" name="DriftChamber" type="DriftChamber" readout="DriftChamberHitsCollection" vis="VisibleBlue" sensitive="true" region="DriftChamberRegion">
+ <detector id="7" name="DriftChamber" type="DriftChamber" readout="DriftChamberHitsCollection" vis="VisibleBlue" sensitive="true" region="DriftChamberRegion">
<envelope vis="SeeThrough">
<shape type="BooleanShape" operation="Union" material="Air">
<shape type="Tube" rmin="SDT_radius_min" rmax="909*mm" dz="SDT_half_length" />
diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
index 36783972a794fab389b648b8b36f772a3cf7ff06..6f88f5898faf2ce84fb3ae4d1f4410fcb521a0bb 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
@@ -74,14 +74,12 @@ StatusCode DCHDigiAlg::initialize()
m_tuple->addItem( "dca" , m_n_digi,m_dca ).ignore();
m_tuple->addItem( "hit_dE" , m_n_digi,m_hit_dE ).ignore();
m_tuple->addItem( "hit_dE_dx", m_n_digi,m_hit_dE_dx ).ignore();
- }
- else { // did not manage to book the N tuple....
- error() << " Cannot book N-tuple:" << long( m_tuple ) << endmsg;
- return StatusCode::FAILURE;
+ } else { // did not manage to book the N tuple....
+ info() << " Cannot book N-tuple:" << long( m_tuple ) << endmsg;
}
}
}
- std::cout<<"DCHDigiAlg::initialized"<< std::endl;
+ info()<<"DCHDigiAlg::initialized"<<endmsg;
return GaudiAlgorithm::initialize();
}
@@ -93,7 +91,7 @@ StatusCode DCHDigiAlg::execute()
edm4hep::TrackerHitCollection* Vec = w_DigiDCHCol.createAndPut();
edm4hep::MCRecoTrackerAssociationCollection* AssoVec = w_AssociationCol.createAndPut();
const edm4hep::SimTrackerHitCollection* SimHitCol = r_SimDCHCol.get();
- std::cout<<"input sim hit size="<< SimHitCol->size() <<std::endl;
+ debug()<<"input sim hit size="<< SimHitCol->size() <<endmsg;
for( int i = 0; i < SimHitCol->size(); i++ )
{
edm4hep::SimTrackerHit SimHit = SimHitCol->at(i);
@@ -110,7 +108,7 @@ StatusCode DCHDigiAlg::execute()
id_hits_map[id] = vhit ;
}
}
- if(m_WriteAna){
+ if(m_WriteAna && (nullptr!=m_tuple)){
m_n_sim = 0;
m_n_digi = 0 ;
}
@@ -164,7 +162,7 @@ StatusCode DCHDigiAlg::execute()
asso.setSim(iter->second.at(i));
asso.setWeight(iter->second.at(i).getEDep()/tot_edep);
- if(m_WriteAna){
+ if(m_WriteAna && (nullptr!=m_tuple)){
m_simhit_x[m_n_sim] = pos.x();
m_simhit_y[m_n_sim] = pos.y();
m_simhit_z[m_n_sim] = pos.z();
@@ -178,7 +176,7 @@ StatusCode DCHDigiAlg::execute()
trkHit.setPosition (edm4hep::Vector3d(pos_x, pos_y, pos_z));//position of closest sim hit
trkHit.setCovMatrix(std::array<float, 6>{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm
- if(m_WriteAna){
+ if(m_WriteAna && (nullptr!=m_tuple)){
m_chamber [m_n_digi] = chamber;
m_layer [m_n_digi] = layer ;
m_cell [m_n_digi] = cellID;
@@ -195,10 +193,10 @@ StatusCode DCHDigiAlg::execute()
}
- std::cout<<"output digi DCHhit size="<< Vec->size() <<std::endl;
+ debug()<<"output digi DCHhit size="<< Vec->size() <<endmsg;
_nEvt ++ ;
- if(m_WriteAna){
+ if(m_WriteAna && (nullptr!=m_tuple)){
StatusCode status = m_tuple->write();
if ( status.isFailure() ) {
error() << " Cannot fill N-tuple:" << long( m_tuple ) << endmsg;
diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.h b/Digitisers/DCHDigi/src/DCHDigiAlg.h
index e703998271af669f23c15a1e5482bd0d516199ad..5ddf59a5426a86358f0b2d775cef8946e46ddbc1 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.h
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.h
@@ -81,7 +81,7 @@ protected:
// Input collections
DataHandle<edm4hep::SimTrackerHitCollection> r_SimDCHCol{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this};
// Output collections
- DataHandle<edm4hep::TrackerHitCollection> w_DigiDCHCol{"DigiDCHitsCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::TrackerHitCollection> w_DigiDCHCol{"DigiDCHitCollection", Gaudi::DataHandle::Writer, this};
DataHandle<edm4hep::MCRecoTrackerAssociationCollection> w_AssociationCol{"DCHitAssociationCollection", Gaudi::DataHandle::Writer, this};
};
diff --git a/Examples/options/sim_fit_CRD.py b/Examples/options/sim_fit_CRD.py
new file mode 100644
index 0000000000000000000000000000000000000000..e877aeb4f5eef6c197cc7a939baba95e76ea46e0
--- /dev/null
+++ b/Examples/options/sim_fit_CRD.py
@@ -0,0 +1,278 @@
+#!/usr/bin/env python
+from Gaudi.Configuration import *
+
+##############################################################################
+# Event service
+##############################################################################
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc")
+
+##############################################################################
+# Random seed
+##############################################################################
+from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
+# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
+rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() # The default engine in Geant4
+rndmengine.SetSingleton = True
+rndmengine.Seeds = [10]
+
+from Configurables import MarlinEvtSeeder
+evtseeder = MarlinEvtSeeder("EventSeeder")
+
+##############################################################################
+# Detector geometry
+##############################################################################
+geometry_option = "CRD_o1_v01/CRD_o1_v01_noECAL.xml"
+
+if not os.getenv("DETCRDROOT"):
+ print("Can't find the geometry. Please setup envvar DETCRDROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+ print("Can't find the compact geometry file: %s"%geometry_path)
+ sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+from Configurables import GearSvc
+gearsvc = GearSvc("GearSvc")
+
+##############################################################################
+# Physics Generator
+##############################################################################
+from Configurables import GenAlgo
+from Configurables import GtGunTool
+from Configurables import StdHepRdr
+from Configurables import SLCIORdr
+from Configurables import HepMCRdr
+from Configurables import GenPrinter
+gun = GtGunTool("GtGunTool")
+gun.Particles = ["mu-"]
+gun.EnergyMins = [100.] # GeV
+gun.EnergyMaxs = [100.] # GeV
+gun.ThetaMins = [0] # deg
+gun.ThetaMaxs = [180] # deg
+gun.PhiMins = [0] # deg
+gun.PhiMaxs = [360] # deg
+# stdheprdr = StdHepRdr("StdHepRdr")
+# stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"
+# lciordr = SLCIORdr("SLCIORdr")
+# lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio"
+# hepmcrdr = HepMCRdr("HepMCRdr")
+# hepmcrdr.Input = "example_UsingIterators.txt"
+
+genprinter = GenPrinter("GenPrinter")
+
+genalg = GenAlgo("GenAlgo")
+genalg.GenTools = ["GtGunTool"]
+#genalg.GenTools = ["StdHepRdr"]
+# genalg.GenTools = ["StdHepRdr", "GenPrinter"]
+# genalg.GenTools = ["SLCIORdr", "GenPrinter"]
+# genalg.GenTools = ["HepMCRdr", "GenPrinter"]
+
+##############################################################################
+# Detector Simulation
+##############################################################################
+from Configurables import DetSimSvc
+detsimsvc = DetSimSvc("DetSimSvc")
+
+from Configurables import DetSimAlg
+detsimalg = DetSimAlg("DetSimAlg")
+#detsimalg.VisMacs = ["vis.mac"]
+detsimalg.RunCmds = [
+# "/tracking/verbose 1",
+]
+detsimalg.AnaElems = [
+ # example_anatool.name()
+ # "ExampleAnaElemTool"
+ "Edm4hepWriterAnaElemTool"
+]
+detsimalg.RootDetElem = "WorldDetElemTool"
+
+from Configurables import AnExampleDetElemTool
+example_dettool = AnExampleDetElemTool("AnExampleDetElemTool")
+
+from Configurables import CalorimeterSensDetTool
+from Configurables import DriftChamberSensDetTool
+
+calo_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool")
+driftchamber_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool")
+
+##############################################################################
+# DedxSimTool
+##############################################################################
+# dedxoption = "DummyDedxSimTool"
+dedxoption = "BetheBlochEquationDedxSimTool"
+
+driftchamber_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import DummyDedxSimTool
+from Configurables import BetheBlochEquationDedxSimTool
+
+if dedxoption == "DummyDedxSimTool":
+ dedx_simtool = DummyDedxSimTool("DummyDedxSimTool")
+elif dedxoption == "BetheBlochEquationDedxSimTool":
+ dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool")
+ dedx_simtool.material_Z = 2
+ dedx_simtool.material_A = 4
+ dedx_simtool.scale = 10
+ dedx_simtool.resolution = 0.0001
+
+##############################################################################
+# DCHDigiAlg
+##############################################################################
+from Configurables import DCHDigiAlg
+dCHDigiAlg = DCHDigiAlg("DCHDigiAlg")
+dCHDigiAlg.readout = "DriftChamberHitsCollection"
+dCHDigiAlg.drift_velocity = 40#um/ns
+dCHDigiAlg.mom_threshold = 0 #GeV
+dCHDigiAlg.WriteAna = True
+
+##############################################################################
+# PODIO
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "CRD-o1-v01-sim00.root"
+out.outputCommands = ["keep *"]
+
+
+##############################################################################
+# Silicon digitization
+##############################################################################
+vxdhitname = "VXDTrackerHits"
+sithitname = "SITTrackerHits"
+sitspname = "SITSpacePoints"
+tpchitname = "TPCTrackerHits"
+sethitname = "SETTrackerHits"
+setspname = "SETSpacePoints"
+ftdspname = "FTDSpacePoints"
+ftdhitname = "FTDTrackerHits"
+from Configurables import PlanarDigiAlg
+digiVXD = PlanarDigiAlg("VXDDigi")
+digiVXD.SimTrackHitCollection = "VXDCollection"
+digiVXD.TrackerHitCollection = vxdhitname
+digiVXD.ResolutionU = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004]
+digiVXD.ResolutionV = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004]
+
+digiSIT = PlanarDigiAlg("SITDigi")
+digiSIT.IsStrip = 1
+digiSIT.SimTrackHitCollection = "SITCollection"
+digiSIT.TrackerHitCollection = sithitname
+digiSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation"
+digiSIT.ResolutionU = [0.007]
+digiSIT.ResolutionV = [0.000]
+
+digiSET = PlanarDigiAlg("SETDigi")
+digiSET.IsStrip = 1
+digiSET.SimTrackHitCollection = "SETCollection"
+digiSET.TrackerHitCollection = sethitname
+digiSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation"
+digiSET.ResolutionU = [0.007]
+digiSET.ResolutionV = [0.000]
+
+digiFTD = PlanarDigiAlg("FTDDigi")
+digiFTD.SimTrackHitCollection = "FTDCollection"
+digiFTD.TrackerHitCollection = ftdhitname
+digiFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation"
+digiFTD.ResolutionU = [0.003, 0.003, 0.007, 0.007, 0.007, 0.007, 0.007, 0.007]
+digiFTD.ResolutionV = [0.003, 0.003, 0, 0, 0, 0, 0, 0 ]
+#digiFTD.OutputLevel = DEBUG
+
+from Configurables import SpacePointBuilderAlg
+spSIT = SpacePointBuilderAlg("SITBuilder")
+spSIT.TrackerHitCollection = sithitname
+spSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation"
+spSIT.SpacePointCollection = sitspname
+spSIT.SpacePointAssociationCollection = "SITSpacePointAssociation"
+#spSIT.OutputLevel = DEBUG
+
+spFTD = SpacePointBuilderAlg("FTDBuilder")
+spFTD.TrackerHitCollection = ftdhitname
+spFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation"
+spFTD.SpacePointCollection = ftdspname
+spFTD.SpacePointAssociationCollection = "FTDSpacePointAssociation"
+#spFTD.OutputLevel = DEBUG
+
+
+##############################################################################
+# Silicon reconstruction
+##############################################################################
+from Configurables import TrackSystemSvc
+tracksystemsvc = TrackSystemSvc("TrackSystemSvc")
+
+from Configurables import SiliconTrackingAlg
+tracking = SiliconTrackingAlg("SiliconTracking")
+tracking.HeaderCol = "EventHeader"
+tracking.VTXHitCollection = vxdhitname
+tracking.SITHitCollection = sitspname
+tracking.FTDPixelHitCollection = ftdhitname
+tracking.FTDSpacePointCollection = ftdspname
+tracking.SITRawHitCollection = sithitname
+tracking.FTDRawHitCollection = ftdhitname
+tracking.UseSIT = 1
+tracking.SmoothOn = 0
+#tracking.OutputLevel = DEBUG
+
+from Configurables import ForwardTrackingAlg
+forward = ForwardTrackingAlg("ForwardTracking")
+forward.FTDPixelHitCollection = ftdhitname
+forward.FTDSpacePointCollection = ftdspname
+forward.FTDRawHitCollection = ftdhitname
+forward.Chi2ProbCut = 0.0
+forward.HitsPerTrackMin = 3
+forward.BestSubsetFinder = "SubsetSimple"
+forward.Criteria = ["Crit2_DeltaPhi","Crit2_StraightTrackRatio","Crit3_3DAngle","Crit3_ChangeRZRatio","Crit3_IPCircleDist","Crit4_3DAngleChange","Crit4_DistToExtrapolation",
+ "Crit2_DeltaRho","Crit2_RZRatio","Crit3_PT"]
+forward.CriteriaMin = [0, 0.9, 0, 0.995, 0, 0.8, 0, 20, 1.002, 0.1, 0, 0.99, 0, 0.999, 0, 0.99, 0]
+forward.CriteriaMax = [30, 1.02, 10, 1.015, 20, 1.3, 1.0, 150, 1.08, 99999999, 0.8, 1.01, 0.35, 1.001, 1.5, 1.01, 0.05]
+#forward.OutputLevel = DEBUG
+
+from Configurables import TrackSubsetAlg
+subset = TrackSubsetAlg("TrackSubset")
+subset.TrackInputCollections = ["ForwardTracks", "SiTracks"]
+subset.RawTrackerHitCollections = [vxdhitname, sithitname, ftdhitname, sitspname, ftdspname]
+subset.TrackSubsetCollection = "SubsetTracks"
+#subset.OutputLevel = DEBUG
+
+##############################################################################
+# TruthTrackerAlg
+##############################################################################
+from Configurables import TruthTrackerAlg
+truthTrackerAlg = TruthTrackerAlg("TruthTrackerAlg")
+truthTrackerAlg.SiSubsetTrackCollection = "SubsetTracks"
+
+##############################################################################
+# RecGenfitAlgSDT
+##############################################################################
+from Configurables import RecGenfitAlgSDT
+recGenfitAlgSDT = RecGenfitAlgSDT("RecGenfitAlgSDT")
+recGenfitAlgSDT.debug=10
+
+##############################################################################
+# NTupleSvc
+##############################################################################
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = [
+#"MyTuples DATAFILE='DCH_digi_ana.root' OPT='NEW' TYP='ROOT'",
+ "RecGenfitAlgSDT DATAFILE='fit_SDT.root' OPT='NEW' TYP='ROOT'"]
+
+
+##############################################################################
+# ApplicationMgr
+##############################################################################
+from Configurables import ApplicationMgr
+ApplicationMgr(
+ TopAlg = [genalg, detsimalg, digiVXD, digiSIT, digiSET, digiFTD, spSIT,
+ spFTD, tracking, forward, subset, dCHDigiAlg, truthTrackerAlg,
+ recGenfitAlgSDT, out],
+ EvtSel = 'NONE',
+ EvtMax = 1,
+ ExtSvc = [rndmengine, dsvc, ntsvc, evtseeder, geosvc, gearsvc, tracksystemsvc],
+ HistogramPersistency = "ROOT",
+ OutputLevel=DEBUG
+)
diff --git a/Examples/options/tut_detsim_digi_fit_DC.py b/Examples/options/tut_detsim_digi_fit_DC.py
new file mode 100644
index 0000000000000000000000000000000000000000..52a7c026596dcee1047590b74bf08c513b1966a3
--- /dev/null
+++ b/Examples/options/tut_detsim_digi_fit_DC.py
@@ -0,0 +1,199 @@
+#!/usr/bin/env python
+
+import os
+print(os.environ["DD4HEP_LIBRARY_PATH"])
+import sys
+# sys.exit(0)
+
+from Gaudi.Configuration import *
+
+##############################################################################
+# Random Number Svc
+##############################################################################
+from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
+
+# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
+rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() # The default engine in Geant4
+rndmengine.SetSingleton = True
+rndmengine.Seeds = [42]
+
+# rndmgensvc = RndmGenSvc("RndmGenSvc")
+# rndmgensvc.Engine = rndmengine.name()
+
+
+##############################################################################
+# Event Data Svc
+##############################################################################
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc")
+
+
+##############################################################################
+# Geometry Svc
+##############################################################################
+
+# geometry_option = "CepC_v4-onlyTracker.xml"
+geometry_option = "det.xml"
+
+if not os.getenv("DETDRIFTCHAMBERROOT"):
+ print("Can't find the geometry. Please setup envvar DETDRIFTCHAMBERROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETDRIFTCHAMBERROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+ print("Can't find the compact geometry file: %s"%geometry_path)
+ sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+##############################################################################
+# Physics Generator
+##############################################################################
+from Configurables import GenAlgo
+from Configurables import GtGunTool
+from Configurables import StdHepRdr
+from Configurables import SLCIORdr
+from Configurables import HepMCRdr
+from Configurables import GenPrinter
+
+gun = GtGunTool("GtGunTool")
+# gun.Particles = ["pi+"]
+# gun.EnergyMins = [100.] # GeV
+# gun.EnergyMaxs = [100.] # GeV
+
+gun.Particles = ["e-"]
+
+# gun.PositionXs = [100.] # mm
+# gun.PositionYs = [100.] # mm
+# gun.PositionZs = [0.] # mm
+
+
+gun.EnergyMins = [10.] # GeV
+gun.EnergyMaxs = [10.] # GeV
+
+gun.ThetaMins = [90] # rad; 45deg
+gun.ThetaMaxs = [90] # rad; 45deg
+
+gun.PhiMins = [0] # rad; 0deg
+gun.PhiMaxs = [360] # rad; 360deg
+#gun.PhiMins = [0] # rad; 0deg
+#gun.PhiMaxs = [360] # rad; 360deg
+
+# stdheprdr = StdHepRdr("StdHepRdr")
+# stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"
+
+# lciordr = SLCIORdr("SLCIORdr")
+# lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio"
+
+# hepmcrdr = HepMCRdr("HepMCRdr")
+# hepmcrdr.Input = "example_UsingIterators.txt"
+
+genprinter = GenPrinter("GenPrinter")
+
+genalg = GenAlgo("GenAlgo")
+genalg.GenTools = ["GtGunTool"]
+# genalg.GenTools = ["StdHepRdr"]
+# genalg.GenTools = ["StdHepRdr", "GenPrinter"]
+# genalg.GenTools = ["SLCIORdr", "GenPrinter"]
+# genalg.GenTools = ["HepMCRdr", "GenPrinter"]
+
+##############################################################################
+# Detector Simulation
+##############################################################################
+from Configurables import DetSimSvc
+
+detsimsvc = DetSimSvc("DetSimSvc")
+
+# from Configurables import ExampleAnaElemTool
+# example_anatool = ExampleAnaElemTool("ExampleAnaElemTool")
+
+from Configurables import DetSimAlg
+
+detsimalg = DetSimAlg("DetSimAlg")
+
+if int(os.environ.get("VIS", 0)):
+ detsimalg.VisMacs = ["vis.mac"]
+
+detsimalg.RunCmds = [
+# "/tracking/verbose 1",
+]
+detsimalg.AnaElems = [
+ # example_anatool.name()
+ # "ExampleAnaElemTool"
+ "Edm4hepWriterAnaElemTool"
+]
+detsimalg.RootDetElem = "WorldDetElemTool"
+
+from Configurables import AnExampleDetElemTool
+example_dettool = AnExampleDetElemTool("AnExampleDetElemTool")
+
+from Configurables import CalorimeterSensDetTool
+from Configurables import DriftChamberSensDetTool
+
+calo_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool")
+driftchamber_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool")
+
+# dedxoption = "DummyDedxSimTool"
+dedxoption = "BetheBlochEquationDedxSimTool"
+
+driftchamber_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import DummyDedxSimTool
+from Configurables import BetheBlochEquationDedxSimTool
+
+if dedxoption == "DummyDedxSimTool":
+ dedx_simtool = DummyDedxSimTool("DummyDedxSimTool")
+elif dedxoption == "BetheBlochEquationDedxSimTool":
+ dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool")
+ dedx_simtool.material_Z = 2
+ dedx_simtool.material_A = 4
+ dedx_simtool.scale = 10
+ dedx_simtool.resolution = 0.0001
+
+##############################################################################
+from Configurables import DCHDigiAlg
+dCHDigiAlg = DCHDigiAlg("DCHDigiAlg")
+
+##############################################################################
+# POD I/O
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "rec_DCH.root"
+out.outputCommands = ["keep *"]
+
+##############################################################################
+# TruthTrackerAlg
+##############################################################################
+from Configurables import TruthTrackerAlg
+truthTrackerAlg = TruthTrackerAlg("TruthTrackerAlg")
+
+##############################################################################
+# GenfitAlg
+##############################################################################
+from Configurables import RecGenfitAlgDC
+recGenfitAlgDC = RecGenfitAlgDC("RecGenfitAlgDC")
+
+##############################################################################
+# NTupleSvc
+##############################################################################
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = ["MyTuples DATAFILE='dCHDigiAlg.root' OPT='NEW' TYP='ROOT'",
+ "RecGenfitAlgDC DATAFILE='recGenfitAlgDC.root' OPT='NEW' TYP='ROOT'"]
+
+##############################################################################
+# ApplicationMgr
+##############################################################################
+
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, truthTrackerAlg,
+ recGenfitAlgDC,out],
+ EvtSel = 'NONE',
+ EvtMax = 10,
+ ExtSvc = [rndmengine, dsvc, geosvc, ntsvc],
+ HistogramPersistency = "ROOT",
+ OutputLevel=DEBUG
+)
diff --git a/Reconstruction/CMakeLists.txt b/Reconstruction/CMakeLists.txt
index e506b90c7c3096915bb1d881d3109984958fa517..6dd63b02cda864eb8f587f8a9f8f6bc57f1b5030 100644
--- a/Reconstruction/CMakeLists.txt
+++ b/Reconstruction/CMakeLists.txt
@@ -3,3 +3,4 @@ add_subdirectory(Digi_Calo)
add_subdirectory(PFA)
add_subdirectory(SiliconTracking)
add_subdirectory(Tracking)
+add_subdirectory(RecGenfitAlg)
diff --git a/Reconstruction/RecGenfitAlg/CMakeLists.txt b/Reconstruction/RecGenfitAlg/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..630086dce7a22cef058903d2859aaed2d469959c
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/CMakeLists.txt
@@ -0,0 +1,36 @@
+# RecGenfitAlg
+
+if (GenFit_FOUND)
+gaudi_add_module(RecGenfitAlg
+ SOURCES src/RecGenfitAlgSDT.cpp
+ src/RecGenfitAlgDC.cpp
+ src/GenfitTrack.cpp
+ src/GenfitField.cpp
+ src/GenfitFitter.cpp
+ src/GenfitMaterialInterface.cpp
+ src/GenfitMsg.cpp
+ LINK GearSvc
+ Gaudi::GaudiAlgLib
+ Gaudi::GaudiKernel
+ ${GEAR_LIBRARIES}
+ ${GSL_LIBRARIES}
+ ${LCIO_LIBRARIES}
+ DetSegmentation
+ DetInterface
+ DataHelperLib
+ EDM4HEP::edm4hep
+ EDM4HEP::edm4hepDict
+ GenFit::genfit2
+)
+
+target_include_directories(RecGenfitAlg PUBLIC
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/include
+ $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+
+install(TARGETS RecGenfitAlg
+ EXPORT CEPCSWTargets
+ RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+ LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+ COMPONENT dev)
+endif()
diff --git a/Reconstruction/RecGenfitAlg/README.md b/Reconstruction/RecGenfitAlg/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..c9cddb4c8b0b9285c1a42a72140c953b9f79797a
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/README.md
@@ -0,0 +1,13 @@
+# RecGenfitAlg package
+
+This package is the track fitting algorithms to do track fitting with genfit.
+It's incudingg the interface to genfit fitter, track and measurement:
+* GenfitAlg: The algorithm for drift chamber fitting with genfit
+* GenfitFitter: fitter implementation
+* GenfitTrack: class for create and access of genfit Track and Measurement
+* GenfitField: implementation of genfit AbsField
+* GenfitGeoMaterialInterface: implementation of genfit AbsMaterialInterface
+* GenfitMsg: a class to get Gaudi message instance in Genfit classes.
+* GenfitAlgDC: a class to fit drift chamber
+* GenfitAlgSDT: a class to fit drift chamber+silicon
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..341e52f3d9c1a0e6539b62bec6225410dbbad7d4
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
@@ -0,0 +1,58 @@
+#include "GenfitField.h"
+#include "GenfitMsg.h"
+
+//External
+#include "DD4hep/DD4hepUnits.h"
+
+//genfit
+#include "FieldManager.h"
+#include "ConstField.h"
+
+///-------------------------------------------------------------
+/// GenfitField constructor
+///-------------------------------------------------------------
+GenfitField::GenfitField(dd4hep::OverlayedField dd4hepField)
+:m_dd4hepField(dd4hepField)
+{
+ genfit::FieldManager::getInstance()->init(this);
+}
+
+
+///-------------------------------------------------------------
+/// GenfitField get field value by TVector3
+///-------------------------------------------------------------
+/// unit in genfit is cm and kGauss
+TVector3 GenfitField::get(const TVector3& pos) const
+{
+ double B[3]={1e9,1e9,1e9};
+ get(pos.X(),pos.Y(),pos.Z(),B[0],B[1],B[2]);
+ return TVector3(B[0],B[1],B[2]);
+}
+
+///-------------------------------------------------------------
+/// GenfitField get field value by double
+///-------------------------------------------------------------
+/// unit in genfit for position is cm and Bxyz is kGauss
+/// unit in DD4hepUnits kilogaus and cm? FIXME
+void
+GenfitField::get(const double& posX, const double& posY, const double& posZ,
+ double& Bx, double& By, double& Bz) const
+{
+ /// get field from dd4hepField
+ const dd4hep::Direction& field=m_dd4hepField.magneticField(
+ dd4hep::Position(posX/dd4hep::cm,posY/dd4hep::cm,posZ/dd4hep::cm));
+ //m_dd4hepField.magneticField(pos,B);
+
+ Bx=field.X()/dd4hep::kilogauss;
+ By=field.Y()/dd4hep::kilogauss;
+ Bz=field.Z()/dd4hep::kilogauss;
+
+// std::cout<<"GenfitField "
+// <<Form("xyz(%f,%f,%f)cm B(%f,%f,%f)kilogauss",posX,posY,posZ,Bx,By,Bz)
+// <<std::endl;
+}
+
+double GenfitField::getBz(const TVector3& pos) const
+{
+ return get(pos).Z();
+}
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitField.h b/Reconstruction/RecGenfitAlg/src/GenfitField.h
new file mode 100644
index 0000000000000000000000000000000000000000..406d0e6b75ccc7c2f87e23fa0bbd60718d92dfd8
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitField.h
@@ -0,0 +1,38 @@
+///////////////////////////////////////////////////////////
+//// An implementation of genfit AbsBField - GenfitField
+//// Authors:
+//// Yao ZHANG (zhangyao@ihep.ac.cn)
+///////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_GENFITFIELD_H
+#define RECGENFITALG_GENFITFIELD_H
+
+#include "AbsBField.h"
+#include "DD4hep/Fields.h"
+
+class TVector3;
+
+
+/// Calss for translating Field into genfit::AbsBField
+class GenfitField : public genfit::AbsBField{
+ public:
+ GenfitField(dd4hep::OverlayedField dd4hepField);
+ virtual ~GenfitField(){;}
+
+ //Get field value from TVector3
+ TVector3 get(const TVector3& pos) const override;
+
+ //Get field value from doubles
+ void get(const double& posX, const double& posY, const double& posZ,
+ double& Bx, double& By, double& Bz) const override;
+
+ //Get Bz, kilogauss
+ double getBz(const TVector3& pos) const;
+
+ //Get dd4hep field
+ const dd4hep::OverlayedField field() const {return m_dd4hepField;}
+
+ private:
+ dd4hep::OverlayedField m_dd4hepField;
+};
+#endif
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..15d40830f9d3fa8deb71e0ca3e306ee7e5ad7c37
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
@@ -0,0 +1,612 @@
+#include "GenfitFitter.h"
+#include "GenfitTrack.h"
+#include "GenfitMsg.h"
+#include "GenfitField.h"
+#include "GenfitMaterialInterface.h"
+
+//Gaudi
+#include "GaudiKernel/StatusCode.h"
+#include "GaudiKernel/ISvcLocator.h"
+
+//External
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/Detector.h"
+
+//genfit
+#include <Track.h>
+#include <Exception.h>
+#include <FieldManager.h>
+#include <TGeoMaterialInterface.h>
+#include <TGeoManager.h>
+#include <MaterialEffects.h>
+#include <MeasuredStateOnPlane.h>
+#include <KalmanFitter.h>
+#include <KalmanFitterRefTrack.h>
+#include <StateOnPlane.h>
+#include <DAF.h>
+#include <AbsKalmanFitter.h>
+#include <KalmanFitterInfo.h>
+
+//ROOT
+#include <TVector3.h>
+#include <TGeoManager.h>
+
+//STL
+#include <iostream>
+#include <string>
+
+
+GenfitFitter::~GenfitFitter(){
+ delete m_absKalman;
+}
+
+GenfitFitter::GenfitFitter(const char* type, const char* name):
+ m_absKalman(nullptr)
+ ,m_genfitField(nullptr)
+ ,m_geoMaterial(nullptr)
+ ,m_fitterType(type)
+ ,m_name(name)
+ ,m_minIterations(4)
+ ,m_maxIterations(10)
+ ,m_deltaPval(1e-3)
+ ,m_relChi2Change(0.2)
+ ,m_blowUpFactor(1e3)
+ ,m_resetOffDiagonals(true)
+ ,m_blowUpMaxVal(1.e6)
+ ,m_multipleMeasurementHandling(genfit::unweightedClosestToPredictionWire)
+ ,m_maxFailedHits(-1)
+ ,m_deltaWeight(1e-3)
+ ,m_annealingBetaStart(100)
+ ,m_annealingBetaStop(0.01)
+ ,m_annealingNSteps(0.01)
+ ,m_noEffects(false)
+ ,m_energyLossBetheBloch(true)
+ ,m_noiseBetheBloch(true)
+ ,m_noiseCoulomb(true)
+ ,m_energyLossBrems(false)
+ ,m_noiseBrems(false)
+ ,m_ignoreBoundariesBetweenEqualMaterials(true)
+ ,m_mscModelName("GEANE")
+ ,m_debug(0)
+ ,m_hist(0)
+{
+ /// Initialize genfit fitter
+ init();
+}
+
+void GenfitFitter::setField(const GenfitField* field)
+{
+ if(nullptr==m_genfitField) m_genfitField=field;
+}
+
+/// Set geometry for material, use geometry from IOADatabase
+void GenfitFitter::setGeoMaterial(const dd4hep::Detector* dd4hepGeo)
+{
+ if(nullptr==m_geoMaterial){
+ m_geoMaterial=GenfitMaterialInterface::getInstance(dd4hepGeo);
+ }
+}
+
+/// initialize genfit fitter
+int GenfitFitter::init(bool deleteOldFitter)
+{
+ if(deleteOldFitter && m_absKalman) delete m_absKalman;
+
+ GenfitMsg::get()<<MSG::DEBUG<<"Initialize GenfitFitter with "
+ <<m_fitterType<<endmsg;
+
+ if (m_fitterType=="DAFRef") {
+ m_absKalman = new genfit::DAF(true,getDeltaPval(),
+ getConvergenceDeltaWeight());
+ }
+ else if (m_fitterType=="DAF") {
+ m_absKalman = new genfit::DAF(false,getDeltaPval(),
+ getConvergenceDeltaWeight());
+ }
+ else if (m_fitterType=="KalmanFitter") {
+ m_absKalman = new genfit::KalmanFitter(getMaxIterations());
+ }
+ else if (m_fitterType=="KalmanFitterRefTrack") {
+ m_absKalman = new genfit::KalmanFitterRefTrack(getMaxIterations());
+ }
+ else {
+ m_absKalman = nullptr;
+ GenfitMsg::get()<< MSG::DEBUG<<"Fitter type is invalid:"
+ <<m_fitterType<<endmsg;
+ return -1;
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<"Fitter type is "<<m_fitterType<<endmsg;
+ m_absKalman->setDebugLvl(m_debug);
+
+ return 0;
+}
+
+/// Fit a track from a candidate track
+int GenfitFitter::processTrackWithRep(GenfitTrack* track,int repID,bool resort)
+{
+ GenfitMsg::get()<<MSG::DEBUG<< "In ProcessTrackWithRep rep "<<repID<<endmsg;
+ if(getDebug()>2) print("");
+
+ if(getDebug()>0){
+ GenfitMsg::get()<<MSG::DEBUG<<"Print track seed "<<endmsg;
+ track->getTrack()->getStateSeed().Print();
+ }
+ /// Do the fitting
+ try{
+ m_absKalman->processTrackWithRep(track->getTrack(), track->getRep(repID),
+ resort);
+ }catch(genfit::Exception& e){
+ GenfitMsg::get()<<MSG::DEBUG<<"Genfit exception caught "<<endmsg;
+ e.what();
+ return false;
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<"End of ProcessTrackWithRep"<<endmsg;
+ return true;
+} // End of ProcessTrackWithRep
+
+/// Fit a track from a candidate track
+int GenfitFitter::processTrack(GenfitTrack* track, bool resort)
+{
+ GenfitMsg::get()<<MSG::DEBUG<<"In ProcessTrack"<<endmsg;
+ if(getDebug()>2) print("");
+
+ /// Do the fitting
+ try{
+ m_absKalman->processTrack(track->getTrack(),resort);
+ }catch(genfit::Exception& e){
+ GenfitMsg::get()<<MSG::DEBUG<<"Genfit exception caught "<<endmsg;
+ e.what();
+ return false;
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<"End of ProcessTrack"<<endmsg;
+ return true;
+} // End of ProcessTrack
+
+
+void GenfitFitter::setFitterType(const char* val)
+{
+ std::string oldFitterType=m_fitterType;
+ m_fitterType = val;
+ GenfitMsg::get()<<MSG::DEBUG<<"Fitter type is "<<m_fitterType<<endmsg;
+ init(oldFitterType==val);
+}
+
+/// Extrapolate track to the cloest point of approach(POCA) to the wire of hit,
+/// return StateOnPlane of this POCA
+/// inputs
+/// pos,mom ... position & momentum at starting point, unit= [mm]&[GeV/c]
+/// (converted to [cm]&[GeV/c] inside this function)
+/// hit ... destination
+/// outputs poca [mm] (converted from [cm] in this function) ,global
+double GenfitFitter::extrapolateToHit( TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, const GenfitTrack* track,
+ TVector3 pos, TVector3 mom,
+ TVector3 end0,//one end of the hit wire
+ TVector3 end1,//the orhter end of the hit wire
+ int debug,
+ int repID,
+ bool stopAtBoundary,
+ bool calcJacobianNoise)
+{
+
+ return track->extrapolateToHit(poca,pocaDir,pocaOnWire,doca,pos,mom,end0,
+ end1,debug,repID,stopAtBoundary,calcJacobianNoise);
+}//end of ExtrapolateToHit
+
+
+/// Extrapolate the track to the cyliner at fixed raidus
+/// position & momentum as starting point
+/// position and momentum at global coordinate in dd4hepUnit
+/// return trackLength in dd4hepUnit
+ double
+GenfitFitter::extrapolateToCylinder(TVector3& pos, TVector3& mom,
+ GenfitTrack* track, double radius, const TVector3 linePoint,
+ const TVector3 lineDirection, int hitID, int repID,
+ bool stopAtBoundary, bool calcJacobianNoise)
+{
+ double trackLength(1e9*dd4hep::cm);
+ if(!track->getFitStatus(repID)->isFitted()) return trackLength;
+ try{
+ // get track rep
+ genfit::AbsTrackRep* rep = track->getRep(repID);
+ if(nullptr == rep) {
+ GenfitMsg::get()<<MSG::DEBUG<<"In ExtrapolateToCylinder rep is null"
+ <<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+
+ // get track point with fitter info
+ genfit::TrackPoint* tp =
+ track->getTrack()->getPointWithFitterInfo(hitID,rep);
+ if(nullptr == tp) {
+ GenfitMsg::get()<<MSG::DEBUG<<
+ "In ExtrapolateToCylinder TrackPoint is null"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+
+ // get fitted state on plane of this track point
+ genfit::KalmanFittedStateOnPlane* state =
+ static_cast<genfit::KalmanFitterInfo*>(
+ tp->getFitterInfo(rep))->getBackwardUpdate();
+
+ if(nullptr == state){
+ GenfitMsg::get()<<MSG::DEBUG<<"In ExtrapolateToCylinder "
+ <<"no KalmanFittedStateOnPlane in backwardUpdate"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+ rep->setPosMom(*state, pos*(1/dd4hep::cm), mom*(1/dd4hep::GeV));
+
+ /// extrapolate
+ trackLength = rep->extrapolateToCylinder(*state,
+ radius/dd4hep::cm, linePoint*(1/dd4hep::cm), lineDirection,
+ stopAtBoundary, calcJacobianNoise);
+ // get pos&mom at extrapolated point on the cylinder
+ rep->getPosMom(*state,pos,mom);//FIXME exception exist
+ pos = pos*dd4hep::cm;
+ mom = mom*dd4hep::GeV;
+ } catch(genfit::Exception& e){
+ GenfitMsg::get() << MSG::ERROR
+ <<"Exception in GenfitFitter::extrapolateToCylinder "
+ << e.what()<<endmsg;
+ trackLength = 1e9*dd4hep::cm;
+ }
+ return trackLength*dd4hep::cm;
+}
+
+double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom,
+ const GenfitTrack* track,
+ const TVector3& point,
+ int repID,// same with pidType
+ bool stopAtBoundary,
+ bool calcJacobianNoise) const
+{
+ double trackLength(1e9*dd4hep::cm);
+ if(!track->getFitStatus(repID)->isFitted()) return trackLength;
+ try{
+ // get track rep
+ genfit::AbsTrackRep* rep = track->getRep(repID);
+ if(nullptr == rep) {
+ GenfitMsg::get()<<MSG::DEBUG<<"In ExtrapolateToPoint rep "
+ <<repID<<" not exist!"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+
+ /// extrapolate to point
+ //genfit::StateOnPlane state(*(&(track->getTrack()->getFittedState(0,rep))));
+
+ // get track point with fitter info
+ genfit::TrackPoint* tp =
+ track->getTrack()->getPointWithFitterInfo(0,rep);
+ if(nullptr == tp) {
+ GenfitMsg::get()<<MSG::DEBUG<<
+ "In ExtrapolateToPoint TrackPoint is null"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+
+ // get fitted state on plane of this track point
+ genfit::KalmanFittedStateOnPlane* state =
+ static_cast<genfit::KalmanFitterInfo*>(
+ tp->getFitterInfo(rep))->getBackwardUpdate();
+
+ if(nullptr == state) {
+ GenfitMsg::get()<<MSG::DEBUG<<
+ "In ExtrapolateToPoint KalmanFittedStateOnPlane is null"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+ trackLength = rep->extrapolateToPoint(*state,
+ point*(1/dd4hep::cm),stopAtBoundary, calcJacobianNoise);
+ rep->getPosMom(*state,pos,mom);//FIXME exception exist
+ pos = pos*dd4hep::cm;
+ mom = mom*dd4hep::GeV;
+ } catch(genfit::Exception& e){
+ GenfitMsg::get() << MSG::ERROR
+ <<"Exception in GenfitFitter::extrapolateToPoint"
+ << e.what()<<endmsg;
+ trackLength = 1e9*dd4hep::cm;
+ }
+ return trackLength*dd4hep::cm;
+}//end of ExtrapolateToPoint
+
+GenfitFitter& GenfitFitter::operator=(
+ const GenfitFitter& r)
+{
+ m_fitterType = r.m_fitterType;
+ m_minIterations = r.m_minIterations;
+ m_maxIterations = r.m_maxIterations;
+ m_deltaPval = r.m_deltaPval;
+ m_relChi2Change = r.m_relChi2Change;
+ m_blowUpFactor = r.m_blowUpFactor;
+ m_resetOffDiagonals = r.m_resetOffDiagonals;
+ m_blowUpMaxVal = r.m_blowUpMaxVal;
+ m_multipleMeasurementHandling = r.m_multipleMeasurementHandling;
+ m_maxFailedHits = r.m_maxFailedHits;
+ m_annealingNSteps = r.m_annealingNSteps;
+ m_deltaWeight = r.m_deltaWeight;
+ m_annealingBetaStart = r.m_annealingBetaStart;
+ m_annealingBetaStop = r.m_annealingBetaStop;
+ m_noEffects = r.m_noEffects;
+ m_energyLossBetheBloch= r.m_energyLossBetheBloch;
+ m_noiseBetheBloch = r.m_noiseBetheBloch;
+ m_noiseCoulomb = r.m_noiseCoulomb;
+ m_energyLossBrems = r.m_energyLossBrems;
+ m_noiseBrems = r.m_noiseBrems;
+ m_ignoreBoundariesBetweenEqualMaterials
+ = r.m_ignoreBoundariesBetweenEqualMaterials;
+ m_mscModelName = r.m_mscModelName;
+ m_debug = r.m_debug;
+ m_hist = r.m_hist;
+ return *this;
+}
+
+void GenfitFitter::print(const char* name)
+{
+
+ if(0==strlen(name)) name=m_name;
+ //TODO print all fitting parameters
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" GenfitFitter Global Parameters:"<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" Fitter type = " << m_fitterType<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" Fitting Parameters:"<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" MinIteration = " << m_minIterations<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" MaxIteration = " << m_maxIterations<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_deltaPval = " << m_deltaPval<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" Material Effect Parameters:"<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_noEffects = " << m_noEffects<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_energyLossBetheBloch= " << m_energyLossBetheBloch<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_noiseBetheBloch = " << m_noiseBetheBloch<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_noiseCoulomb = " << m_noiseCoulomb<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_energyLossBrems = " << m_energyLossBrems<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_noiseBrems = " << m_noiseBrems<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_ignoreBoundariesBetweenEqualMaterials= "
+ << m_ignoreBoundariesBetweenEqualMaterials<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_mscModelName = " << m_mscModelName<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_debug = " << m_debug<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_hist = " << m_hist<<endmsg;
+
+ if(m_fitterType=="DAF"||m_fitterType=="DAFRef"){
+ genfit::DAF* daf = getDAF();
+
+ if(nullptr != daf){
+ std::ostringstream sBetas;
+ std::vector<double> betas = daf->getBetas();
+ for (unsigned int i=0; i<betas.size(); ++i) {
+ sBetas<<betas.at(i)<<",";
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<" print "<<betas.size()<<
+ " betas: %s "<<sBetas.str()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<m_name
+ << " m_deltaWeight = " << m_deltaWeight<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF maxIterations= "
+ << daf->getMaxIterations()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF minIterations= "
+ << daf->getMinIterations()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF deltaPval= "
+ << daf->getDeltaPval()<<endmsg;
+ }
+ }
+ genfit::KalmanFitterRefTrack* ref = getKalRef();
+ if(nullptr != ref){
+ std::ostringstream sBetas;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF maxIterations= "
+ << ref->getMaxIterations()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF minIterations= "
+ << ref->getMinIterations()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF deltaPval= "
+ << ref->getDeltaPval()<<endmsg;
+ }
+}
+
+///Setters of AbsKalmanFitter
+void GenfitFitter::setMinIterations(unsigned int val)
+{
+ m_absKalman->setMinIterations(val);
+ m_minIterations = val;
+}
+
+void GenfitFitter::setMaxIterations(unsigned int val)
+{
+ if(val<=4) {
+ GenfitMsg::get()<<MSG::ERROR<<"Could NOT set MaxIteration<=4!"<<endmsg;
+ }
+ m_absKalman->setMaxIterations(val);
+ m_maxIterations = val;
+}
+
+void GenfitFitter::setMaxIterationsBetas(double bStart,double bFinal,
+ unsigned int val)
+{
+ m_absKalman->setMaxIterations(val);
+ m_maxIterations = val;
+ if(val<=4) {
+ GenfitMsg::get()<<MSG::ERROR<<"Could NOT set MaxIteration<=4!"<<endmsg;
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<"bStart "<<bStart<<" bFinal "<<bFinal
+ <<" MaxIteration "<<val<<endmsg;
+ // genfit version less than 2.2.0
+ genfit::DAF* daf = dynamic_cast<genfit::DAF*> (m_absKalman);
+ daf->setAnnealingScheme(bStart,bFinal,val);
+}
+
+void GenfitFitter::setDeltaPval(double val)
+{
+ m_absKalman->setDeltaPval(val);
+ m_deltaPval = val;
+}
+
+void GenfitFitter::setRelChi2Change(double val)
+{
+ m_absKalman->setRelChi2Change(val);
+ m_relChi2Change = val;
+}
+
+void GenfitFitter::setBlowUpFactor(double val)
+{
+ m_absKalman->setBlowUpFactor(val);
+ m_blowUpFactor = val;
+}
+
+void GenfitFitter::setResetOffDiagonals(bool val)
+{
+ m_absKalman->setResetOffDiagonals(val);
+ m_resetOffDiagonals = val;
+}
+
+void GenfitFitter::setBlowUpMaxVal(double val)
+{
+ m_absKalman->setBlowUpMaxVal(val);
+ m_blowUpMaxVal = val;
+}
+
+void GenfitFitter::setMultipleMeasurementHandling(
+ genfit::eMultipleMeasurementHandling val)
+{
+ m_absKalman->setMultipleMeasurementHandling(val);
+ m_multipleMeasurementHandling = val;
+}
+
+void GenfitFitter::setMaxFailedHits(int val)
+{
+ m_absKalman->setMaxFailedHits(val);
+ m_maxFailedHits = val;
+}
+
+///DAF setters
+void GenfitFitter::setConvergenceDeltaWeight(double val)
+{
+ genfit::DAF* daf = getDAF();
+ if(nullptr != daf) daf->setConvergenceDeltaWeight(val);
+ m_deltaWeight = val;
+}
+void GenfitFitter::setAnnealingScheme(
+ double bStart, double bFinal, unsigned int nSteps)
+{
+ genfit::DAF* daf = getDAF();
+ if(nullptr != daf) daf->setAnnealingScheme(bStart, bFinal, nSteps);
+ m_annealingBetaStart = bStart;
+ m_annealingBetaStop = bFinal;
+ m_annealingNSteps = nSteps;
+}
+
+//TODO chi2Cuts?
+///Material effects
+void GenfitFitter::setNoEffects(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setNoEffects();
+ m_noEffects = val;
+}
+
+void GenfitFitter::setEnergyLossBetheBloch(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setEnergyLossBetheBloch(val);
+ m_energyLossBetheBloch = val;
+}
+
+void GenfitFitter::setNoiseBetheBloch(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setNoiseBetheBloch(val);
+ m_noiseBetheBloch = val;
+}
+
+void GenfitFitter::setNoiseCoulomb(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setNoiseCoulomb(val);
+ m_noiseCoulomb = val;
+}
+
+void GenfitFitter::setEnergyLossBrems(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setEnergyLossBrems(val);
+ m_energyLossBrems = val;
+}
+
+void GenfitFitter::setNoiseBrems(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setNoiseBrems(val);
+ m_noiseBrems = val;
+}
+
+void GenfitFitter::setIgnoreBoundariesBetweenEqualMaterials(bool val)
+{
+ genfit::MaterialEffects::getInstance()->
+ ignoreBoundariesBetweenEqualMaterials(val);
+ m_ignoreBoundariesBetweenEqualMaterials = val;
+}
+
+void GenfitFitter::setMscModelName(std::string val)
+{
+ genfit::MaterialEffects::getInstance()->setMscModel(val);
+ m_mscModelName = val;
+}
+
+void GenfitFitter::setMaterialDebugLvl(unsigned int val)
+{
+ genfit::MaterialEffects::getInstance()->setDebugLvl(val);
+}
+
+///Set GenfitFitter parameters
+void GenfitFitter::setDebug(unsigned int val)
+{
+ GenfitMsg::get()<<MSG::DEBUG<<"set fitter debugLvl "<<val<<endmsg;
+ m_absKalman->setDebugLvl(val);
+ if(nullptr!=getDAF()) getDAF()->setDebugLvl(val);
+ if(val>10) genfit::MaterialEffects::getInstance()->setDebugLvl(val);
+ m_debug = val;
+}
+
+void GenfitFitter::setHist(unsigned int val) {m_hist = val;}
+
+genfit::DAF* GenfitFitter::getDAF()
+{
+ genfit::DAF* daf = nullptr;
+ try{
+ daf = dynamic_cast<genfit::DAF*> (m_absKalman);
+ }catch(...){
+ GenfitMsg::get()<< MSG::ERROR
+ << "dynamic_cast m_rom AbsFitter to DAF m_ailed!"<<endmsg;
+ return nullptr;
+ }
+ return daf;
+}
+
+genfit::KalmanFitterRefTrack* GenfitFitter::getKalRef()
+{
+ genfit::KalmanFitterRefTrack* ref=nullptr;
+ try{
+ ref = dynamic_cast<genfit::KalmanFitterRefTrack*> (m_absKalman);
+ }catch(...){
+ GenfitMsg::get()<< MSG::ERROR
+ << "dynamic_cast m_rom AbsFitter to KalmanFitterRefTrack m_ailed!"
+ <<endmsg;
+ }
+ return ref;
+}
+
+void GenfitFitter::initHist(std::string name)
+{
+ GenfitMsg::get()<<MSG::DEBUG<<"GenfitFitter::initHist "<<name<<endmsg;
+ //getDAF()->initHist(name);
+}
+
+void GenfitFitter::writeHist()
+{
+ GenfitMsg::get()<<MSG::DEBUG<<"GenfitFitter::writeHist "<<endmsg;
+ //getDAF()->writeHist();
+}
+
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitFitter.h b/Reconstruction/RecGenfitAlg/src/GenfitFitter.h
new file mode 100644
index 0000000000000000000000000000000000000000..9ca4664456dab2f9484f1b1e9ad2c8d33fc5474a
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitFitter.h
@@ -0,0 +1,200 @@
+//////////////////////////////////////////////////////////////////
+///
+/// This is an interface of to call genfit fitter
+///
+/// In this file, including:
+/// a genfit fitter class
+/// Fitter type is fixed after creation
+/// Magnetic field and geometry for material effect are singletons in genfit.
+/// They will be set before fitting and before each time call track rep.
+///
+/// Authors:
+/// Yao ZHANG (zhangyao@ihep.ac.cn)
+///
+//////////////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_GENFITFITTER_H
+#define RECGENFITALG_GENFITFITTER_H
+
+#include "AbsKalmanFitter.h"
+#include <string>
+
+class GenfitField;
+class GenfitMaterialInterface;
+class GenfitTrack;
+class TVector3;
+class TGeoManager;
+class ISvcLocator;
+namespace genfit{
+ class AbsKalmanFitter;
+ class KalmanFitterRefTrack;
+ class DAF;
+}
+namespace dd4hep{
+ class OverlayedField;
+ class Detector;
+}
+
+
+/// The GenfitTrack class
+class GenfitFitter{
+ public:
+ /// Type of fitters are :DAFRef,DAF,KalmanFitter,KalmanFitterRefTrack
+ GenfitFitter(const char* type="DAFRef",const char* name="GenfitFitter");
+ virtual ~GenfitFitter();
+
+ /// Magnetic field and geometry for material effect in genfit
+ /// please SET before use !!!!
+ void setField(const GenfitField* field);
+ /// please SET before use !!!!
+ void setGeoMaterial(const dd4hep::Detector* dd4hepGeo);
+
+ /// Main fitting function
+ int processTrack(GenfitTrack* track, bool resort=false);
+
+ /// fitting with rep
+ int processTrackWithRep(GenfitTrack* track,int repID=0,
+ bool resort=false);
+
+ /// Extrapolate the track to the CDC hit
+ /// Output: poca pos and dir and poca distance to the hit wire
+ /// Input: genfit track, pos and mom, two ends of a wire
+ /// pos, and mom are position & momentum at starting point
+ double extrapolateToHit(TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, const GenfitTrack* track,
+ TVector3 pos, TVector3 mom, TVector3 end0, TVector3 end1,
+ int debug=0, int repID=0, bool stopAtBoundary=false,
+ bool calcJacobianNoise=false);
+
+ /// Extrapolate the track to the cyliner at fixed raidus
+ /// Output: pos and mom at fixed radius
+ /// Input: genfitTrack, radius of cylinder at center of the origin,
+ /// repID, stopAtBoundary and calcAverageState
+ double extrapolateToCylinder(TVector3& pos, TVector3& mom,
+ GenfitTrack* track, double radius, const TVector3 linePoint,
+ const TVector3 lineDirection, int hitID =0, int repID=0,
+ bool stopAtBoundary=false, bool calcJacobianNoise=false);
+
+ /// Extrapolate the track to the point
+ /// Output: pos and mom of POCA point to point
+ /// Input: genfitTrack,point,repID,stopAtBoundary and calcAverageState
+ /// repID same with pidType
+ double extrapolateToPoint(TVector3& pos, TVector3& mom,
+ const GenfitTrack* genfitTrack, const TVector3& point,
+ int repID=0, bool stopAtBoundary = false,
+ bool calcJacobianNoise = false) const;
+
+ /// setters of fitter properties
+ void setFitterType(const char* val);
+ void setMinIterations(unsigned int val);
+ void setMaxIterations(unsigned int val);
+ void setMaxIterationsBetas(double bStart,double bFinal,unsigned int val);
+ void setDeltaPval(double val);
+ void setRelChi2Change(double val);
+ void setBlowUpFactor(double val);
+ void setResetOffDiagonals(bool val);
+ void setBlowUpMaxVal(double val);
+ void setMultipleMeasurementHandling(
+ genfit::eMultipleMeasurementHandling val);
+ void setMaxFailedHits(int val);
+ void setConvergenceDeltaWeight(double val);
+ void setAnnealingScheme(double bStart,double bFinal,unsigned int nSteps);
+ //TODO chi2cut?
+ void setNoEffects(bool val);
+ void setEnergyLossBetheBloch(bool val);
+ void setNoiseBetheBloch(bool val);
+ void setNoiseCoulomb(bool val);
+ void setEnergyLossBrems(bool val);
+ void setNoiseBrems(bool val);
+ void setIgnoreBoundariesBetweenEqualMaterials(bool val);
+ void setMscModelName(std::string val);
+ void setMaterialDebugLvl(unsigned int val);
+ void setDebug(unsigned int val);
+ void setHist(unsigned int val);
+
+ /// getters of fitter properties
+ std::string getFitterType() const {return m_fitterType;}
+ unsigned int getMinIterations() const { return m_minIterations; }
+ unsigned int getMaxIterations() const { return m_maxIterations; }
+ double getDeltaPval() const { return m_deltaPval; }
+ double getRelChi2Change() const { return m_relChi2Change; }
+ double getBlowUpFactor() const { return m_blowUpFactor; }
+ double getBlowUpMaxVal() const { return m_blowUpMaxVal; }
+ bool getResetOffDiagonals() const { return m_resetOffDiagonals; }
+ genfit::eMultipleMeasurementHandling getMultipleMeasurementHandling()
+ const { return m_multipleMeasurementHandling; }
+ int getMaxFailedHits() const { return m_maxFailedHits; }
+ double getConvergenceDeltaWeight() const { return m_deltaWeight; }
+ float getAnnealingBetaStart() const {return m_annealingBetaStart;}
+ float getAnnealingBetaStop() const {return m_annealingBetaStop;}
+ bool getNoEffects(){return m_noEffects;}
+ bool getEnergyLossBetheBloch(){return m_energyLossBetheBloch;}
+ bool getNoiseBetheBloch(){return m_noiseBetheBloch;}
+ bool getNoiseCoulomb(){return m_noiseCoulomb;}
+ bool getEnergyLossBrems(){return m_energyLossBrems;}
+ bool getNoiseBrems(){return m_noiseBrems;}
+ bool getIgnoreBoundariesBetweenEqualMaterials()
+ {return m_ignoreBoundariesBetweenEqualMaterials;}
+ std::string getMscModelName(){return m_mscModelName;}
+ int getDebug() const {return m_debug;}
+ int getHist() const {return m_hist;}
+
+ /// Printer
+ void print(const char* name="");
+ void initHist(std::string name);
+ void writeHist();
+
+ const char* getName(void) const {return m_name;}
+
+ private:
+ GenfitFitter(const GenfitFitter&){};
+ GenfitFitter& operator=(const GenfitFitter&);
+
+ /// Initialze fitter and setting fitter parameter
+ int init(bool deleteOldFitter=false);
+
+ /// Get DAF fitter
+ genfit::DAF* getDAF();
+ /// Get KalmanFitterRefTrack
+ genfit::KalmanFitterRefTrack* getKalRef();
+
+ genfit::AbsKalmanFitter* m_absKalman;/// kalman fitter object
+
+ const GenfitField* m_genfitField;//pointer to genfit field
+ GenfitMaterialInterface* m_geoMaterial;//pointer to genfit geo
+
+ ///fitting method: DAFRef,DAF,KalmanFitter,KalmanFitterRefTrack
+ std::string m_fitterType;
+
+ const char* m_name;
+
+ /// control parameters of fitter
+ unsigned int m_minIterations; /// minimum number of iterations
+ unsigned int m_maxIterations; /// maximum number of iterations
+ double m_deltaPval; /// delta pVal that converged
+ double m_relChi2Change; ///
+ double m_blowUpFactor; ///
+ bool m_resetOffDiagonals; ///
+ double m_blowUpMaxVal; ///
+ genfit::eMultipleMeasurementHandling m_multipleMeasurementHandling;///
+ int m_maxFailedHits; ///
+ double m_deltaWeight; /// delta weight that converged
+ double m_annealingBetaStart;/// start beta for annealing
+ double m_annealingBetaStop; /// stop beta for annealing
+ unsigned int m_annealingNSteps; /// n steps
+
+ /// control parAmeters of material effects
+ bool m_noEffects;
+ bool m_energyLossBetheBloch;
+ bool m_noiseBetheBloch;
+ bool m_noiseCoulomb;
+ bool m_energyLossBrems;
+ bool m_noiseBrems;
+ bool m_ignoreBoundariesBetweenEqualMaterials;
+ std::string m_mscModelName;
+
+ int m_debug; /// debug
+ int m_hist; /// hist
+};
+
+#endif
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..73a59316e0dbbe52c1b1743eaa6bb6bcc4b4b6c9
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp
@@ -0,0 +1,353 @@
+///
+/// Authors: ZHANG Yao (zhangyao@ihep.ac.cn)
+///
+
+#include "GenfitMaterialInterface.h"
+//Gaudi
+#include "GaudiKernel/Bootstrap.h"
+#include "GaudiKernel/SmartIF.h"
+#include "GaudiKernel/ISvcLocator.h"
+//CEPCSW
+#include "DetInterface/IGeomSvc.h"
+//Externals
+#include "DD4hep/Detector.h"
+//#include "DD4hep/DD4hepUnits.h"
+//genfit
+#include "Exception.h"
+//ROOT
+#include "TGeoManager.h"
+#include "TGeoNode.h"
+#include "TGeoMedium.h"
+#include "TGeoMaterial.h"
+#include "TGeoManager.h"
+#include "MaterialEffects.h"
+#include "Material.h"
+//STL
+#include <assert.h>
+#include <math.h>
+
+GenfitMaterialInterface* GenfitMaterialInterface::m_instance = nullptr;
+
+//#define DEBUG_GENFITGEO 1
+
+double MeanExcEnergy_get(int Z);
+double MeanExcEnergy_get(TGeoMaterial*);
+
+GenfitMaterialInterface::GenfitMaterialInterface(
+ const dd4hep::Detector* dd4hepGeo)
+{
+ assert(nullptr!=dd4hepGeo);
+ m_geoManager=&(dd4hepGeo->manager());
+ assert(nullptr!=m_geoManager);
+}
+
+GenfitMaterialInterface* GenfitMaterialInterface::getInstance(
+ const dd4hep::Detector* dd4hepGeo)
+{
+ if(nullptr == m_instance){
+ m_instance = new GenfitMaterialInterface(dd4hepGeo);
+ }
+ genfit::MaterialEffects::getInstance()->init(m_instance);
+ return m_instance;
+}
+
+void GenfitMaterialInterface::destruct(){
+ delete m_instance;
+}
+
+TGeoManager* GenfitMaterialInterface::getGeoManager()
+{
+ return m_geoManager;
+}
+
+ bool
+GenfitMaterialInterface::initTrack(double posX, double posY,
+ double posZ, double dirX, double dirY, double dirZ)
+{
+ //debugLvl_ = 1;
+#ifdef DEBUG_GENFITGEO
+ std::cout << "GenfitMaterialInterface::initTrack. \n";
+ std::cout << "Pos "; TVector3(posX, posY, posZ).Print();
+ std::cout << "Dir "; TVector3(dirX, dirY, dirZ).Print();
+#endif
+
+ // Move to the new point.
+ bool result = !isSameLocation(posX, posY, posZ, kTRUE);
+
+ // Set the intended direction.
+ setCurrentDirection(dirX, dirY, dirZ);
+
+ if (debugLvl_ > 0) {
+ std::cout << " GenfitMaterialInterface::initTrack at \n";
+ std::cout << " position: "; TVector3(posX, posY, posZ).Print();
+ std::cout << " direction: "; TVector3(dirX, dirY, dirZ).Print();
+ }
+
+ return result;
+}
+
+
+genfit::Material GenfitMaterialInterface::getMaterialParameters()
+{
+ TGeoMaterial* mat =
+ getGeoManager()->GetCurrentVolume()->GetMedium()->GetMaterial();
+ //Scalar density; /// Density in g / cm^3
+ //Scalar Z; /// Atomic number
+ //Scalar A; /// Mass number in g / mol
+ //Scalar radiationLength; /// Radiation Length in cm
+ //Scalar mEE; /// Mean excitaiton energy in eV
+ //Material from CEPCSW is NOT follow the dd4hep?FIXME
+
+ //std::cout<<__FILE__<<" "<<__LINE__<<" yzhang debug material "<<std::endl;
+ //mat->Print();
+ return genfit::Material(mat->GetDensity(), mat->GetZ(), mat->GetA(),
+ mat->GetRadLen(), MeanExcEnergy_get(mat));
+}
+
+ double
+GenfitMaterialInterface::findNextBoundary(const genfit::RKTrackRep* rep,
+ const genfit::M1x7& stateOrig,
+ double sMax, // signed
+ bool varField)
+{
+ //debugLvl_ = 1;
+ // cm, distance limit beneath which straight-line steps are taken.
+ const double delta(1.E-2);
+ const double epsilon(1.E-1); // cm, allowed upper bound on arch
+ // deviation from straight line
+
+ genfit::M1x3 SA;
+ genfit::M1x7 state7, oldState7;
+ oldState7 = stateOrig;
+
+ int stepSign(sMax < 0 ? -1 : 1);
+
+ double s = 0; // trajectory length to boundary
+
+ const unsigned maxIt = 300;
+ unsigned it = 0;
+
+ // Initialize the geometry to the current location (set by caller).
+ findNextBoundary(fabs(sMax) - s);
+ double safety = getSafeDistance(); // >= 0
+ double slDist = getStep();
+
+ // this should not happen, but it happens sometimes.
+ // The reason are probably overlaps in the geometry.
+ // Without this "hack" many small steps with size of minStep will be made,
+ // until eventually the maximum number of iterations are exceeded
+ // and the extrapolation fails
+ if (slDist < safety)
+ slDist = safety;
+ double step = slDist;
+
+ if (debugLvl_ > 0)
+ std::cout << " safety = " << safety << "; step, slDist = " << step << "\n";
+
+ while (1) {
+ if (++it > maxIt){
+ genfit::Exception exc("GenfitMaterialInterface::findNextBoundary ==> maximum number of iterations exceeded",__LINE__,__FILE__);
+ exc.setFatal();
+ throw exc;
+ }
+
+ // No boundary in sight?
+ if (s + safety > fabs(sMax)) {
+ if (debugLvl_ > 0)
+ std::cout << " next boundary is further away than sMax \n";
+ return stepSign*(s + safety); //sMax;
+ }
+
+ // Are we at the boundary?
+ if (slDist < delta) {
+ if (debugLvl_ > 0)
+ std::cout << " very close to the boundary -> return"
+ << " stepSign*(s + slDist) = "
+ << stepSign << "*(" << s + slDist << ")\n";
+ return stepSign*(s + slDist);
+ }
+
+ // We have to find whether there's any boundary on our path.
+
+ // Follow curved arch, then see if we may have missed a boundary.
+ // Always propagate complete way from original start to avoid
+ // inconsistent extrapolations.
+ state7 = stateOrig;
+ rep->RKPropagate(state7, nullptr, SA, stepSign*(s + step), varField);
+
+ // Straight line distance虏 between extrapolation finish and
+ // the end of the previously determined safe segment.
+ double dist2 = (pow(state7[0] - oldState7[0], 2)
+ + pow(state7[1] - oldState7[1], 2)
+ + pow(state7[2] - oldState7[2], 2));
+ // Maximal lateral deviation虏.
+ double maxDeviation2 = 0.25*(step*step - dist2);
+
+ if (step > safety
+ && maxDeviation2 > epsilon*epsilon) {
+ // Need to take a shorter step to reliably estimate material,
+ // but only if we didn't move by safety. In order to avoid
+ // issues with roundoff we check 'step > safety' instead of
+ // 'step != safety'. If we moved by safety, there couldn't have
+ // been a boundary that we missed along the path, but we could
+ // be on the boundary.
+
+ // Take a shorter step, but never shorter than safety.
+ step = std::max(step / 2, safety);
+ } else {
+ getGeoManager()->PushPoint();
+ bool volChanged = initTrack(state7[0], state7[1], state7[2],
+ stepSign*state7[3], stepSign*state7[4],
+ stepSign*state7[5]);
+
+ if (volChanged) {
+ if (debugLvl_ > 0)
+ std::cout << " volChanged\n";
+ // Move back to start.
+ getGeoManager()->PopPoint();
+
+ // Extrapolation may not take the exact step length we asked
+ // for, so it can happen that a requested step < safety takes
+ // us across the boundary. This is then the best estimate we
+ // can get of the distance to the boundary with the stepper.
+ if (step <= safety) {
+ if (debugLvl_ > 0)
+ std::cout <<
+ " step <= safety, return stepSign*(s + step) = "
+ << stepSign*(s + step) << "\n";
+ return stepSign*(s + step);
+ }
+
+ // Volume changed during the extrapolation. Take a shorter
+ // step, but never shorter than safety.
+ step = std::max(step / 2, safety);
+ } else {
+ // we're in the new place, the step was safe, advance
+ s += step;
+
+ oldState7 = state7;
+ getGeoManager()->PopDummy(); // Pop stack, but stay in place.
+
+ findNextBoundary(fabs(sMax) - s);
+ safety = getSafeDistance();
+ slDist = getStep();
+
+ // this should not happen, but it happens sometimes.
+ // The reason are probably overlaps in the geometry.
+ // Without this "hack" many small steps with size of minStep will be made,
+ // until eventually the maximum number of iterations are exceeded and the extrapolation fails
+ if (slDist < safety)
+ slDist = safety;
+ step = slDist;
+ if (debugLvl_ > 0){
+ std::cout << " safety = " << safety
+ << "; step, slDist = " << step << "\n";
+ }
+ }
+ }
+ }
+}
+
+
+/*
+ Reference for elemental mean excitation energies at:
+http://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html
+
+Code ported from GEANT 3
+*/
+
+const int MeanExcEnergy_NELEMENTS = 93; // 0 = vacuum, 1 = hydrogen, 92 = uranium
+const double MeanExcEnergy_vals[MeanExcEnergy_NELEMENTS] = {
+ 1.e15, 19.2, 41.8, 40.0, 63.7, 76.0, 78., 82.0,
+ 95.0, 115.0, 137.0, 149.0, 156.0, 166.0, 173.0, 173.0,
+ 180.0, 174.0, 188.0, 190.0, 191.0, 216.0, 233.0, 245.0,
+ 257.0, 272.0, 286.0, 297.0, 311.0, 322.0, 330.0, 334.0,
+ 350.0, 347.0, 348.0, 343.0, 352.0, 363.0, 366.0, 379.0,
+ 393.0, 417.0, 424.0, 428.0, 441.0, 449.0, 470.0, 470.0,
+ 469.0, 488.0, 488.0, 487.0, 485.0, 491.0, 482.0, 488.0,
+ 491.0, 501.0, 523.0, 535.0, 546.0, 560.0, 574.0, 580.0,
+ 591.0, 614.0, 628.0, 650.0, 658.0, 674.0, 684.0, 694.0,
+ 705.0, 718.0, 727.0, 736.0, 746.0, 757.0, 790.0, 790.0,
+ 800.0, 810.0, 823.0, 823.0, 830.0, 825.0, 794.0, 827.0,
+ 826.0, 841.0, 847.0, 878.0, 890.0, };
+// Logarithms of the previous table, used to calculate mixtures.
+const double MeanExcEnergy_logs[MeanExcEnergy_NELEMENTS] = {
+ 34.5388, 2.95491, 3.7329, 3.68888, 4.15418, 4.33073, 4.35671, 4.40672,
+ 4.55388, 4.74493, 4.91998, 5.00395, 5.04986, 5.11199, 5.15329, 5.15329,
+ 5.19296, 5.15906, 5.23644, 5.24702, 5.25227, 5.37528, 5.45104, 5.50126,
+ 5.54908, 5.6058, 5.65599, 5.69373, 5.73979, 5.77455, 5.79909, 5.81114,
+ 5.85793, 5.84932, 5.8522, 5.83773, 5.86363, 5.8944, 5.90263, 5.93754,
+ 5.97381, 6.03309, 6.04973, 6.05912, 6.08904, 6.10702, 6.15273, 6.15273,
+ 6.1506, 6.19032, 6.19032, 6.18826, 6.18415, 6.19644, 6.17794, 6.19032,
+ 6.19644, 6.21661, 6.25958, 6.28227, 6.30262, 6.32794, 6.35263, 6.36303,
+ 6.38182, 6.41999, 6.44254, 6.47697, 6.4892, 6.51323, 6.52796, 6.54247,
+ 6.5582, 6.57647, 6.58893, 6.60123, 6.61473, 6.62936, 6.67203, 6.67203,
+ 6.68461, 6.69703, 6.71296, 6.71296, 6.72143, 6.71538, 6.67708, 6.7178,
+ 6.71659, 6.73459, 6.7417, 6.77765, 6.79122, };
+
+
+double
+MeanExcEnergy_get(int Z, bool logs = false) {
+ assert(Z >= 0 && Z < MeanExcEnergy_NELEMENTS);
+ if (logs)
+ return MeanExcEnergy_logs[Z];
+ else
+ return MeanExcEnergy_vals[Z];
+}
+
+
+double
+MeanExcEnergy_get(TGeoMaterial* mat) {
+ if (mat->IsMixture()) {
+ // The mean excitation energy is calculated as the geometric mean
+ // of the mEEs of the components, weighted for abundance.
+ double logMEE = 0.;
+ double denom = 0.;
+ TGeoMixture* mix = (TGeoMixture*)mat;
+ for (int i = 0; i < mix->GetNelements(); ++i) {
+ int index = int(mix->GetZmixt()[i]);
+ double weight = mix->GetWmixt()[i] * mix->GetZmixt()[i] / mix->GetAmixt()[i];
+ logMEE += weight * MeanExcEnergy_get(index, true);
+ denom += weight;
+ }
+ logMEE /= denom;
+ return exp(logMEE);
+ }
+
+ // not a mixture
+ int index = int(mat->GetZ());
+ return MeanExcEnergy_get(index, false);
+}
+
+double GenfitMaterialInterface::getSafeDistance()
+{
+ return getGeoManager()->GetSafeDistance();
+}
+
+double GenfitMaterialInterface::getStep()
+{
+ return getGeoManager()->GetSafeDistance();
+}
+
+TGeoNode* GenfitMaterialInterface::findNextBoundary(double stepmax,
+ const char* path,bool frombdr)
+{
+ return getGeoManager()->FindNextBoundary(stepmax,path,frombdr);
+}
+
+bool GenfitMaterialInterface::isSameLocation(double posX, double posY,
+ double posZ, bool change)
+{
+ //std::cout<<" MatInt "<<__LINE__<<" posXYZ*dd4hep::cm "<<posX*dd4hep::cm
+ //<<" posY "<<posY*dd4hep::cm<<" posZ "<<posZ*dd4hep::cm<<std::endl;
+ return getGeoManager()->IsSameLocation(posX,posY,
+ posZ,change);
+}
+
+void GenfitMaterialInterface::setCurrentDirection(double nx, double ny,
+ double nz)
+{
+ return getGeoManager()->SetCurrentDirection(nx, ny, nz);
+
+}
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..83da9fc1e4cbb372e7c7d7ee35efde3bd5d22646
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h
@@ -0,0 +1,74 @@
+/////////////////////////////////////////////////////////////////////////
+// An implementation of genfit AbsMaterialInterface
+//
+// Author:
+// Yao ZHANG (zhangyao@ihep.ac.cn)
+//
+//
+/////////////////////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_GENFITMATERIALINTERFACE_H
+#define RECGENFITALG_GENFITMATERIALINTERFACE_H
+
+#include "AbsMaterialInterface.h"
+
+class RKTrackRep;
+class TGeoManager;
+class TGeoNode;
+namespace dd4hep{
+ class Detector;
+}
+
+/**
+ * @brief MaterialInterface implementation for use with ROOT's TGeoManager.
+ */
+class GenfitMaterialInterface : public genfit::AbsMaterialInterface{
+ public:
+ GenfitMaterialInterface(const dd4hep::Detector* dd4hepGeo);
+ virtual ~GenfitMaterialInterface(){};
+ static GenfitMaterialInterface* getInstance(
+ const dd4hep::Detector* dd4hepGeo);
+ void destruct();
+
+ //Set Detector pointer of dd4hep
+ void setDetector(dd4hep::Detector*);
+
+ /** @brief Initialize the navigator at given position and with given
+ direction. Returns true if the volume changed.
+ */
+ bool initTrack(double posX, double posY, double posZ,
+ double dirX, double dirY, double dirZ) override;
+
+ genfit::Material getMaterialParameters() override;
+
+
+ /** @brief Make a step (following the curvature) until step length
+ * sMax or the next boundary is reached. After making a step to a
+ * boundary, the position has to be beyond the boundary, i.e. the
+ * current material has to be that beyond the boundary. The actual
+ * step made is returned.
+ */
+ double findNextBoundary(const genfit::RKTrackRep* rep,
+ const genfit::M1x7& state7,
+ double sMax,
+ bool varField = true) override;
+
+ // ClassDef(GenfitMaterialInterface, 1);
+
+ private:
+ static GenfitMaterialInterface* m_instance;
+ TGeoManager* getGeoManager();
+ TGeoManager* m_geoManager;
+ double getSafeDistance();
+ double getStep();
+ TGeoNode* findNextBoundary(double stepmax, const char* path="",
+ bool frombdr=false);
+ bool isSameLocation(double posX, double posY, double posZ,
+ bool change=false);
+ void setCurrentDirection(double nx, double ny, double nz);
+
+};
+
+/** @} */
+
+#endif // GenfitMaterialInterface
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp b/Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6799d68e0ed8120274c7b5a0d6ac24eb8fb4b078
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp
@@ -0,0 +1,19 @@
+#include "GenfitMsg.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/IMessageSvc.h"
+
+MsgStream* GenfitMsg::m_log = nullptr;
+
+MsgStream GenfitMsg::get(){
+ if(nullptr == m_log){
+ SmartIF<IMessageSvc> msgSvc(
+ Gaudi::svcLocator()->service<IMessageSvc>("MessageSvc"));
+ m_log = new MsgStream(msgSvc, "GenfitMsg");
+ (*m_log) << MSG::DEBUG << "initialize GenfitMsg" << endmsg;
+ }
+ return (*m_log);
+}
+
+GenfitMsg::~GenfitMsg(){
+ delete m_log;
+}
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMsg.h b/Reconstruction/RecGenfitAlg/src/GenfitMsg.h
new file mode 100644
index 0000000000000000000000000000000000000000..43808b42b425ef29245445c2144c8a255facda4b
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMsg.h
@@ -0,0 +1,26 @@
+//////////////////////////////////////////////////////////////////////
+///
+/// This is a class for genfit interfaces classes to access Gaudi log
+/// system
+///
+/// Authors:
+/// Yao ZHANG(zhangyao@ihep.ac.cn)
+///
+//////////////////////////////////////////////////////////////////////
+#ifndef RECGENFITALG_GENFITMSG_H
+#define RECGENFITALG_GENFITMSG_H
+
+#include "GaudiKernel/MsgStream.h"
+
+class GenfitMsg {
+ public:
+ static MsgStream get();
+ private:
+ GenfitMsg(){};
+ virtual ~GenfitMsg();
+ static MsgStream* m_log;
+
+};
+
+#endif
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f210dc9f69c88a0cbc6407b18911c181c700accf
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
@@ -0,0 +1,838 @@
+#include "GenfitTrack.h"
+#include "GenfitMsg.h"
+#include "GenfitField.h"
+
+//CEPCSW
+#include "DataHelper/HelixClass.h"
+#include "UTIL/ILDConf.h"
+
+//Externals
+#include "DD4hep/DD4hepUnits.h"
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/Track.h"
+#include "edm4hep/TrackerHitConst.h"
+#include "edm4hep/SimTrackerHit.h"
+#include "edm4hep/ReconstructedParticle.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/Vector3d.h"
+#include "DetSegmentation/GridDriftChamber.h"
+
+//genfit
+#include "Track.h"
+#include "MeasuredStateOnPlane.h"
+#include "RKTrackRep.h"
+#include "TrackPoint.h"
+#include "StateOnPlane.h"
+#include "KalmanFitterInfo.h"
+#include "KalmanFittedStateOnPlane.h"
+#include "AbsTrackRep.h"
+#include "FitStatus.h"
+#include "SpacepointMeasurement.h"
+#include "WireMeasurementNew.h"
+
+//ROOT
+#include "TRandom.h"
+#include "TVector3.h"
+#include "TLorentzVector.h"
+#include "TMatrixDSym.h"
+
+//cpp
+#include <cfloat>
+
+const int GenfitTrack::s_PDG[2][5]
+={{-11,-13,211,321,2212},{11,13,-211,-321,-2212}};
+
+ bool
+sortDCHit(edm4hep::ConstSimTrackerHit hit1,edm4hep::ConstSimTrackerHit hit2)
+{
+ //std::cout<<"hit1"<<hit1<<std::endl;
+ //std::cout<<"hit2"<<hit2<<std::endl;
+ bool isEarly=hit1.getTime()<hit2.getTime();
+ return isEarly;
+}
+
+ GenfitTrack::GenfitTrack(const GenfitField* genfitField, const dd4hep::DDSegmentation::GridDriftChamber* seg, const char* name)
+:m_name(name),m_track(nullptr),m_reps(),m_debug(0),
+m_genfitField(genfitField),m_gridDriftChamber(seg)
+{
+
+}
+
+GenfitTrack::~GenfitTrack()
+{
+ ///Note: track reps and points will be deleted in the destructor of track
+ ///implemented in genfit::Track::Clear()
+ delete m_track;
+}
+
+/// create a Genfit track from track state, without trackRep
+/// Initialize track with seed states
+/// NO unit conversion here
+bool GenfitTrack::createGenfitTrack(int pdgType,int charge,
+ TLorentzVector posInit, TVector3 momInit, TMatrixDSym covMInit)
+{
+ TVectorD seedState(6);
+ TMatrixDSym seedCov(6);
+
+ //for(int i = 0; i < 6; ++i) {
+ // for(int j = 0; j < 6; ++j) {
+ // seedCov(i,j)=covMInit(i,j);
+ // }
+ //}
+ //yzhang FIXME
+ //seed position
+ for(int i = 0; i < 3; ++i) {
+ seedState(i)=posInit[i];
+ //yzhang TODO from covMInit to seedCov
+ double resolution = 0.1;//*dd4hep::mm/dd4hep::cm;
+ seedCov(i,i)=resolution*resolution;
+ if(i==2) seedCov(i,i)=0.5*0.5;
+ }
+ //seed momentum
+ for(int i = 3; i < 6; ++i){
+ //seedState(i)=momInit[i-3]*(dd4hep::GeV);
+ seedState(i)=momInit[i-3];
+ //yzhang TODO from covMInit to seedCov
+ seedCov(i,i)=0.01;//pow(resolution / sqrt(3),2);
+ }
+
+ if(nullptr==m_track) m_track=new genfit::Track();
+ m_track->setStateSeed(seedState);
+ m_track->setCovSeed(seedCov);
+
+ /// new a track representation and add to the track
+ int chargeId=0;
+ charge>0 ? chargeId=0 : chargeId=1;
+
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" CreateGenfitTrack seed pos("
+ <<seedState[0]<<" "<<seedState[1]<<" "<<seedState[2]<<")cm ("
+ <<seedState[3]<<" "<<seedState[4]<<" "<<seedState[5]<<")GeV charge "
+ <<charge<<" pdg "<<s_PDG[chargeId][pdgType]<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<"seedCov "<<endmsg;
+
+ addTrackRep(s_PDG[chargeId][pdgType]);
+
+ if(m_debug>0) seedCov.Print();
+ return true;
+}
+
+///Create a Genfit track with MCParticle, unit conversion here
+bool GenfitTrack::createGenfitTrackFromMCParticle(int pidType,
+ const edm4hep::MCParticle& mcParticle, double eventStartTime)
+{
+ ///get track parameters from McParticle
+ edm4hep::Vector3d mcPocaPos = mcParticle.getVertex();//mm
+ edm4hep::Vector3f mcPocaMom = mcParticle.getMomentum();//GeV
+ GenfitMsg::get()<<MSG::DEBUG<<"seedPos poca "<< mcPocaPos.x
+ <<" "<<mcPocaPos.y<<" "<<mcPocaPos.z<<" mm "<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<"seedMom poca "<< mcPocaMom.x
+ <<" "<<mcPocaMom.y<<" "<<mcPocaMom.z<<" GeV "<<endmsg;
+
+ ///Pivot to first layer to avoid correction of beam pipe
+ edm4hep::Vector3d firstLayerPos(1e9,1e9,1e9);
+ edm4hep::Vector3f firstLayerMom(1e9,1e9,1e9);
+ pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom);
+
+ //TODO convert unit
+ ///Get seed position and momentum
+ TLorentzVector seedPos(firstLayerPos.x,firstLayerPos.y,firstLayerPos.z,
+ eventStartTime);
+ TVector3 seedMom(firstLayerMom.x,firstLayerMom.y,firstLayerMom.z);
+ GenfitMsg::get()<<MSG::DEBUG<<"seedPos "<< firstLayerPos.x
+ <<" "<<firstLayerPos.y<<" "<<firstLayerPos.z<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<"seedMom "<< firstLayerMom.x
+ <<" "<<firstLayerMom.y<<" "<<firstLayerMom.z<<endmsg;
+
+ ///Get error matrix of seed track
+ TMatrixDSym covM(5);//FIXME, TODO
+
+ ///Create a genfit track with seed
+ GenfitMsg::get()<<MSG::DEBUG<<"createGenfitTrack " ;
+ if(!GenfitTrack::createGenfitTrack(pidType,mcParticle.getCharge(),
+ seedPos,seedMom,covM)){
+ GenfitMsg::get()<<MSG::ERROR<<"GenfitTrack"
+ <<" Error in createGenfitTrack" <<endmsg;
+ return false;
+ }else{
+ GenfitMsg::get()<<MSG::DEBUG<<"GenfitTrack "
+ <<"createGenfitTrackFromMCParticle track created" <<endmsg;
+ }
+ return true;
+}//end of createGenfitTrackFromMCParticle
+
+///Create a Genfit track with MCParticle, unit conversion here
+bool GenfitTrack::createGenfitTrackFromEDM4HepTrack(int pidType,
+ const edm4hep::Track& track, double eventStartTime)
+{
+ //std::cout<<__FILE__<<" "<<__LINE__<<" bz kilogauss "<<m_genfitField->getBz({0.,0.,0.})/dd4hep::kilogauss<<std::endl;
+ //std::cout<<__FILE__<<" "<<__LINE__<<" bz tesla "<<m_genfitField->getBz({0.,0.,0.})/dd4hep::tesla<<std::endl;
+ //std::cout<<__FILE__<<" "<<__LINE__<<" bz "<<m_genfitField->getBz({0.,0.,0.})<< dd4hep::kilogauss <<" "<<dd4hep::tesla<<std::endl;
+ //TODO
+ //pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom);
+ //Get track parameters
+ edm4hep::TrackState trackStat=track.getTrackStates(0);//FIXME?
+ HelixClass helixClass;
+ helixClass.Initialize_Canonical(trackStat.phi,trackStat.D0,
+ trackStat.Z0,trackStat.omega,trackStat.tanLambda,
+ m_genfitField->getBz({0.,0.,0.})*dd4hep::kilogauss/dd4hep::tesla);
+ TLorentzVector posInit(helixClass.getReferencePoint()[0],
+ helixClass.getReferencePoint()[1],
+ helixClass.getReferencePoint()[2],eventStartTime);
+ posInit.SetX(posInit.X()*dd4hep::mm);
+ posInit.SetY(posInit.Y()*dd4hep::mm);
+ posInit.SetZ(posInit.Z()*dd4hep::mm);
+ TVector3 momInit(helixClass.getMomentum()[0],
+ helixClass.getMomentum()[1],helixClass.getMomentum()[2]);
+ momInit.SetX(momInit.x()*dd4hep::GeV);
+ momInit.SetY(momInit.y()*dd4hep::GeV);
+ momInit.SetZ(momInit.z()*dd4hep::GeV);
+ TMatrixDSym covMInit;
+ if(!createGenfitTrack(pidType,
+ int(trackStat.omega/fabs(trackStat.omega)),//charge,FIXME?
+ posInit,momInit,covMInit)){
+ return false;
+ }
+ return true;
+}
+
+/// Add a 3d SpacepointMeasurement on TrackerHit
+bool GenfitTrack::addSpacePointTrakerHit(edm4hep::ConstTrackerHit& hit,
+ int hitID)
+{
+ edm4hep::Vector3d pos=hit.getPosition();
+ TVectorD p(3);
+ p[0]=pos.x*dd4hep::mm;
+ p[1]=pos.y*dd4hep::mm;
+ p[2]=pos.z*dd4hep::mm;
+
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" addSpacePointTrakerHit"<<hitID
+ <<"pos "<<p[0]<<" "<<p[1]<<" "<<p[2]<<" cm"<<endmsg;
+ /// New a SpacepointMeasurement
+ double cov[6];
+ for(int i=0;i<6;i++) {
+ cov[i]=hit.getCovMatrix(i);
+ GenfitMsg::get()<<MSG::DEBUG<<"cov "<<cov[i]<<endmsg;
+ }
+ TMatrixDSym hitCov(3);//FIXME?
+ //in SimpleDigi/src/PlanarDigiAlg.cpp
+ //cov[0] = u_direction[0];
+ //cov[1] = u_direction[1];
+ //cov[2] = resU;
+ //cov[3] = v_direction[0];
+ //cov[4] = v_direction[1];
+ //cov[5] = resV;
+ hitCov(0,0)=cov[2]*dd4hep::mm*dd4hep::mm;
+ hitCov(1,1)=cov[2]*dd4hep::mm*dd4hep::mm;
+ hitCov(2,2)=cov[2]*dd4hep::mm*dd4hep::mm;
+
+ genfit::SpacepointMeasurement* sMeas =
+ new genfit::SpacepointMeasurement(p,hitCov,(int) hit.getCellID(),hitID,
+ nullptr);
+ genfit::TrackPoint* trackPoint = new genfit::TrackPoint(sMeas,m_track);
+ m_track->insertPoint(trackPoint);
+
+ GenfitMsg::get()<<MSG::DEBUG<<"end of addSpacePointTrakerHit"<<endmsg;
+ return true;
+}
+
+/// Add a 3d SpacepointMeasurement with MC truth position smeared by sigma
+bool GenfitTrack::addSpacePointMeasurement(const TVectorD& pos,
+ double sigma, int detID, int hitID, bool smear)
+{
+ double sigma_t=sigma*dd4hep::mm;
+ /// Convert from CEPCSW unit to genfit unit, cm
+ TVectorD pos_t(3);
+ pos_t(0)=pos(0)*dd4hep::mm;
+ pos_t(1)=pos(1)*dd4hep::mm;
+ pos_t(2)=pos(2)*dd4hep::mm;
+
+ /// smear hit position with same weight
+ TVectorD pos_smeared(3);
+ for (int i=0;i<3;i++){
+ pos_smeared[i]=pos_t(i);
+ if(smear) pos_smeared[i]+=gRandom->Gaus(0,sigma_t/TMath::Sqrt(3.));
+ }
+
+ /// New a SpacepointMeasurement
+ TMatrixDSym hitCov(3);
+ hitCov(0,0)=sigma_t*sigma_t;
+ hitCov(1,1)=sigma_t*sigma_t;
+ hitCov(2,2)=sigma_t*sigma_t;
+
+ GenfitMsg::get()<< MSG::DEBUG<<m_name<<" addSpacePointMeasurement detID "
+ <<detID<<" hitId "<<hitID<<" " <<pos_t[0]<<" "<<pos_t[1]<<" "<<pos_t[2]
+ <<" cm smeared "<<pos_smeared[0]<<" "<<pos_smeared[1]<<" "
+ <<pos_smeared[2]<<" sigma_t "<<sigma_t<<" cm"<<endmsg;
+
+ genfit::SpacepointMeasurement* sMeas =
+ new genfit::SpacepointMeasurement(pos_smeared,hitCov,detID,hitID,nullptr);
+ genfit::TrackPoint* trackPoint = new genfit::TrackPoint(sMeas,m_track);
+ m_track->insertPoint(trackPoint);
+
+ return true;
+}
+
+
+/// Add a WireMeasurement, no Unit conversion here
+void GenfitTrack::addWireMeasurement(double driftDistance,
+ double driftDistanceError, const TVector3& endPoint1,
+ const TVector3& endPoint2, int lrAmbig, int detID, int hitID)
+{
+ try{
+ /// New a WireMeasurement
+ genfit::WireMeasurementNew* wireMeas = new genfit::WireMeasurementNew(
+ driftDistance, driftDistanceError, endPoint1, endPoint2, detID,
+ hitID, nullptr);
+ wireMeas->setMaxDistance(0.5);//0.5 cm FIXME
+ wireMeas->setLeftRightResolution(lrAmbig);
+
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" Add wire measurement(cm) "<<hitID
+ <<" ep1("<<endPoint1[0]<<" "<<endPoint1[1]<<" "<<endPoint1[2]
+ <<") ep2("<<endPoint2[0]<<" "<<endPoint2[1]<<" "<<endPoint2[2]
+ <<") drift "<<driftDistance<<" driftErr "<<driftDistanceError
+ <<" lr "<<lrAmbig<<" detId "<<detID << " hitId "<< hitID
+ <<endmsg;
+
+ ///New a TrackPoint,create connection between meas. and trackPoint
+ genfit::TrackPoint* trackPoint=new genfit::TrackPoint(wireMeas,m_track);
+ wireMeas->setTrackPoint(trackPoint);
+
+ m_track->insertPoint(trackPoint);
+
+ }catch(genfit::Exception& e){
+ GenfitMsg::get()<<MSG::DEBUG<<m_name
+ <<"Add wire measurement exception"<<endmsg;
+ e.what();
+ }
+}//end of addWireMeasurementOnTrack
+
+//Add wire measurement on wire, unit conversion here
+bool GenfitTrack::addWireMeasurementOnTrack(edm4hep::Track& track,double sigma)
+{
+ for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){
+ edm4hep::ConstTrackerHit hit=track.getTrackerHits(iHit);
+
+ double driftVelocity=40;//FIXME, TODO, um/ns
+ double driftDistance=hit.getTime()*driftVelocity*dd4hep::um;
+ TVector3 endPointStart(0,0,0);
+ TVector3 endPointEnd(0,0,0);
+ m_gridDriftChamber->cellposition(hit.getCellID(),endPointStart,
+ endPointEnd);
+ int lrAmbig=0;
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" time "<<hit.getTime()
+ <<" driftVelocity " <<driftVelocity<<std::endl;
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" wire pos " <<endPointStart.X()
+ <<" "<<endPointStart.Y()<<" " <<endPointStart.Z()<<" "
+ <<endPointEnd.X()<<" " <<endPointEnd.Y()<<" "
+ <<endPointEnd.Z()<<" " <<std::endl;
+ endPointStart.SetX(endPointStart.x()*dd4hep::cm);
+ endPointStart.SetY(endPointStart.y()*dd4hep::cm);
+ endPointStart.SetZ(endPointStart.z()*dd4hep::cm);
+ endPointEnd.SetX(endPointEnd.x()*dd4hep::cm);
+ endPointEnd.SetY(endPointEnd.y()*dd4hep::cm);
+ endPointEnd.SetZ(endPointEnd.z()*dd4hep::cm);
+ addWireMeasurement(driftDistance,sigma*dd4hep::cm,endPointStart,
+ endPointEnd,lrAmbig,hit.getCellID(),iHit);
+ }
+ return true;
+}//end of addWireMeasurementOnTrack of Track
+
+/// Get MOP
+bool GenfitTrack::getMOP(int hitID,
+ genfit::MeasuredStateOnPlane& mop, genfit::AbsTrackRep* trackRep) const
+{
+ if(nullptr == trackRep) trackRep = getRep();
+ try{
+ mop = m_track->getFittedState(hitID,trackRep);
+ }catch(genfit::Exception& e){
+ e.what();
+ return false;
+ }
+ return true;
+}
+
+/// New and add a track representation to track
+genfit::RKTrackRep* GenfitTrack::addTrackRep(int pdg)
+{
+ /// create a new track representation
+ genfit::RKTrackRep* rep = new genfit::RKTrackRep(pdg);
+ m_reps.push_back(rep);
+ m_track->addTrackRep(rep);
+ //try{
+ // genfit::MeasuredStateOnPlane stateInit(rep);
+ // rep->setPosMomCov(stateInit, pos, mom, covM);
+ //}catch(genfit::Exception e){
+ // GenfitMsg::get() << MSG::DEBUG<<m_name<<" Exception in set track status"
+ // <<endmsg ;
+ // std::cout<<e.what()<<std::endl;
+ // return false;
+ //}
+ return rep;
+}
+
+/// Get the position from genfit::Track::getStateSeed
+const TLorentzVector GenfitTrack::getSeedStatePos()const
+{
+ TVectorD seedStat(6);
+ seedStat = m_track->getStateSeed();
+ TVector3 p(seedStat[0],seedStat[1],seedStat[2]);
+ p = p*dd4hep::cm;
+ TLorentzVector pos(p[0],p[1],p[2],9999);//FIXME
+ return pos;
+}
+
+/// Get the momentum from genfit::Track::getStateSeed
+const TVector3 GenfitTrack::getSeedStateMom() const
+{
+ TVectorD seedStat(6); seedStat = m_track->getStateSeed();
+ TVector3 mom(seedStat[3],seedStat[4],seedStat[5]);
+ return mom*dd4hep::GeV;
+}
+
+/// Get the seed states of momentum and position
+void GenfitTrack::getSeedStateMom(TLorentzVector& pos, TVector3& mom) const
+{
+ TVectorD seedStat(6); seedStat = m_track->getStateSeed();
+ mom = TVector3(seedStat[3],seedStat[4],seedStat[5])*dd4hep::GeV;
+ seedStat = m_track->getStateSeed();
+ TVector3 p = TVector3(seedStat[0],seedStat[1],seedStat[2])*dd4hep::cm;
+ pos.SetXYZT(p[0],p[1],p[2],9999);//FIXME time
+}
+
+unsigned int GenfitTrack::getNumPoints() const
+{
+ return m_track->getNumPoints();
+}
+
+/// Test the fit result FIXME
+bool GenfitTrack::fitSuccess(int repID) const
+{
+
+ genfit::FitStatus* fitStatus = m_track->getFitStatus(getRep(repID));
+
+ /// Test fitting converged
+ if (!fitStatus->isFitted()||!fitStatus->isFitConverged()
+ ||fitStatus->isFitConvergedFully()) {
+ GenfitMsg::get() << MSG::DEBUG<<m_name<< "Fitting is failed... isFitted"
+ <<fitStatus->isFitted()<<" , isFitConverged "
+ <<fitStatus->isFitConverged()<<", isFitConvergedFully "
+ <<fitStatus->isFitConvergedFully()<<endmsg;
+ return false;
+ }
+
+ double chi2 = fitStatus->getChi2();
+ double ndf = fitStatus->getNdf();
+ GenfitMsg::get() << MSG::INFO << "Fit result: chi2 "<<chi2 <<" ndf "<<ndf
+ << " chi2/ndf = " << chi2/ndf<<endmsg;
+
+ /// Test fitting chi2
+ if (chi2<= 0) {
+ GenfitMsg::get() << MSG::DEBUG<<m_name<< "Fit chi2<0 (chi2,ndf) = (" <<
+ chi2 << "," << ndf << ")"<<endmsg;
+ return false;
+ }
+ return true;
+}
+
+void GenfitTrack::setDebug(int debug)
+{
+ m_debug = debug;
+ for(unsigned int i=0;i<m_reps.size();i++){ m_reps[i]->setDebugLvl(debug); }
+}
+
+void GenfitTrack::printSeed() const
+{
+ TLorentzVector pos = getSeedStatePos();
+ TVector3 mom = getSeedStateMom();
+ print(pos,mom);
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" NumPoints "<<getNumPoints()<<endmsg;
+}
+
+void GenfitTrack::printFitted(int repID) const
+{
+ TLorentzVector fittedPos;
+ TVector3 fittedMom;
+ TMatrixDSym cov;
+
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<< "printFitted nHit="
+ <<m_track->getNumPoints()<<endmsg;
+ for(unsigned int iHit=0; iHit<m_track->getNumPoints(); iHit++){
+ if (getPosMomCovMOP((int) iHit, fittedPos, fittedMom, cov, repID)){
+ //print(fittedPos,fittedMom,to_string(iHit).c_str());//TODO
+ }else{
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<"Hit "<<iHit
+ <<" have no fitter info"<<endmsg;
+ }
+ }
+}
+
+/// Print track information
+void GenfitTrack::print( TLorentzVector pos, TVector3 mom,
+ const char* comment) const
+{
+ TVector3 pglo = pos.Vect();
+ TVector3 mglo = mom;
+
+ // TODO
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" "<<comment<<endmsg;
+
+ if(m_debug>1){
+ for(unsigned int i=0;i<m_reps.size();i++){ m_reps[i]->Print(); }
+ }
+ //for(unsigned int i=0; i<m_track->getNumPoints(); i++){
+ // m_track->getPoint(i)->print();
+ //}
+}
+
+/// Get position, momentum, cov on plane of hitID-th hit
+bool GenfitTrack::getPosMomCovMOP(int hitID, TLorentzVector& pos,
+ TVector3& mom, TMatrixDSym& cov, int repID) const
+{
+ TVector3 p;
+ genfit::MeasuredStateOnPlane mop;
+ if(!getMOP(hitID,mop,getRep(repID))) return false;
+ mop.getPosMomCov(p,mom,cov);
+ pos.SetVect(p*dd4hep::cm);
+ pos.SetT(9999);//FIXME
+ mom = mom*(dd4hep::GeV);
+ //FIXME
+ //TrackingUtils::CovConvertUnit(cov, dd4hep::cm, dd4hep::GeV);
+ return true;
+}
+
+int GenfitTrack::getNumPointsWithFittedInfo(int repID) const
+{
+ int nHitWithFittedInfo = 0;
+ int nHit = m_track->getNumPointsWithMeasurement();
+ for(int i=0; i<nHit; i++){
+ if(nullptr != m_track->getPointWithFitterInfo(i,getRep(repID))){
+ nHitWithFittedInfo++;
+ }
+ }
+ return nHitWithFittedInfo;
+}
+
+int GenfitTrack::getFittedState(TLorentzVector& pos, TVector3& mom,
+ TMatrixDSym& cov, int repID, bool biased) const
+{
+ //check number of hit with fitter info
+ if(getNumPointsWithFittedInfo(repID)<=2) return 1;
+
+ //get track rep
+ genfit::AbsTrackRep* rep = getRep(repID);
+ if(nullptr == rep) return 2;
+
+ //get first or last measured state on plane
+ genfit::MeasuredStateOnPlane mop;
+ try{
+ mop = m_track->getFittedState(biased);
+ }catch(genfit::Exception& e){
+ std::cout<<" getNumPointsWithFittedInfo "
+ <<getNumPointsWithFittedInfo(repID)
+ <<" no TrackPoint with fitted info "<<std::endl;
+ GenfitMsg::get()<<MSG::DEBUG<<m_name
+ <<"Exception in getFittedState"<<endmsg;
+ std::cout<<e.what()<<std::endl;
+ return 3;
+ }
+
+ //get state
+ TVector3 p;
+ mop.getPosMomCov(p,mom,cov);
+ pos.SetVect(p*dd4hep::cm);
+ pos.SetT(9999);//FIXME
+ mom = mom*(dd4hep::GeV);
+
+ return 0;//success
+}
+
+// Get point with fitter info
+int GenfitTrack::getDetIDWithFitterInfo(int hitID, int idRaw) const
+{
+ return m_track->getPointWithFitterInfo(hitID)->
+ getRawMeasurement(idRaw)->getDetId();
+}
+
+int GenfitTrack::getPDG(int id) const
+{
+ return m_reps[id]->getPDG();
+}
+
+int GenfitTrack::getPDGCharge(int id) const
+{
+ return m_reps[id]->getPDGCharge();
+}
+
+const genfit::FitStatus*
+GenfitTrack::getFitStatus(int repID) const
+{
+ return m_track->getFitStatus(getRep(repID));
+}
+
+/// Extrapolate track to the cloest point of approach(POCA) to the wire of hit,
+/// return StateOnPlane of this POCA
+/// inputs
+/// pos,mom ... position & momentum at starting point, unit= [mm]&[GeV/c]
+/// (converted to [cm]&[GeV/c] inside this function)
+/// hit ... destination
+/// outputs poca [mm] (converted from [cm] in this function) ,global
+double GenfitTrack::extrapolateToHit( TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, TVector3 pos, TVector3 mom,
+ TVector3 end0,//one end of the hit wire
+ TVector3 end1,//the orhter end of the hit wire
+ int debug,
+ int repID,
+ bool stopAtBoundary,
+ bool calcJacobianNoise)const
+{
+
+ //genfit::MeasuredStateOnPlane state = getMOP(iPoint); // better?
+ //genfit::MeasuredStateOnPlane state = getMOP(0); // better?
+ //To do the extrapolation without IHitSelection,above 2 lines are not used.
+ pos = pos*dd4hep::cm;//FIXME
+ mom = mom*dd4hep::GeV;
+
+ //std::cout<<__LINE__<<" extrapolate to Hit pos and mom"<<std::endl;
+ //pos.Print();
+ //mom.Print();
+
+ genfit::AbsTrackRep* rep = new genfit::RKTrackRep(getRep(repID)->getPDG());
+ rep->setDebugLvl(debug);
+ genfit::MeasuredStateOnPlane state(rep);
+ rep->setPosMom(state, pos, mom);
+
+
+ //genfit::MeasuredStateOnPlane state(m_track->getRep(repID));
+ //m_track->getRep(repID)->setPosMom(state, pos, mom);
+
+ //m_track->PrintSeed();
+ double extrapoLen(0);
+ //std::cout<<" wire1 "<<std::endl;
+ //end0.Print();
+ //std::cout<<" wire0 "<<std::endl;
+ //end1.Print();
+ //std::cout<<" state "<<std::endl;
+ //state.Print();
+
+
+ // forth
+ try {
+ //genfit::RKTrackRep* rkRep =
+ // dynamic_cast<genfit::RKTrackRep*> (m_track->getRep(repID));
+ //extrapoLen = rkRep->extrapolateToLine(state, end0, end1, poca,
+ // pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise);
+ extrapoLen = rep->extrapolateToLine(state, end0, end1, poca,
+ pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise);
+ }catch(genfit::Exception& e) {
+ GenfitMsg::get() << MSG::ERROR
+ <<"Exception in GenfigFitter::ExtrapolateToHit"<<e.what()<<endmsg;
+ return extrapoLen;
+ }
+
+ //poca = poca*(dd4hep::cm);
+ //pocaOnWire = pocaOnWire*(dd4hep::cm);
+ pocaOnWire = pocaOnWire;
+ doca = (pocaOnWire-poca).Mag();
+ //TODO: debug pocaOnWire
+ GenfitMsg::get() << MSG::DEBUG << " poca "<<poca.x()<<","<<poca.y()
+ <<" "<<poca.z()<<" doca "<<doca<<endmsg;
+ GenfitMsg::get() << MSG::DEBUG << " pocaOnWire "<<pocaOnWire.x()
+ <<","<<pocaOnWire.y()<<" "<<pocaOnWire.z()<<" doca "<<doca<<endmsg;
+
+ return extrapoLen*(dd4hep::cm);
+}//end of ExtrapolateToHit
+
+
+///Add space point measurement from edm4hep::Track to genfit track
+int GenfitTrack::addSimTrackerHits(const edm4hep::Track& track,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ float sigma,bool smear){
+ //A TrakerHit collection
+ std::vector<edm4hep::ConstSimTrackerHit> sortedDCTrackHitCol;
+
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" VXD "
+ <<lcio::ILDDetID::VXD<<" SIT "
+ <<lcio::ILDDetID::SIT<<" SET "
+ <<lcio::ILDDetID::SET<<" FTD "
+ <<lcio::ILDDetID::FTD<<" "<<endmsg;
+ ///Get TrackerHit on Track
+ int hitID=0;
+ for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){
+ edm4hep::ConstTrackerHit hit=track.getTrackerHits(iHit);
+
+ UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
+ encoder.setValue(hit.getCellID());
+ int detID=encoder[lcio::ILDCellID0::subdet];
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" "<<iHit<<" hit "<<hit
+ <<" detID "<<detID<<endmsg;
+
+ if(detID==lcio::ILDDetID::VXD || detID==lcio::ILDDetID::SIT ||
+ detID==lcio::ILDDetID::SET || detID==lcio::ILDDetID::FTD){
+ if(addSpacePointTrakerHit(hit,hitID)){
+ GenfitMsg::get()<<MSG::DEBUG<<"add slicon space point"<<endmsg;
+ hitID++;
+ }else{
+ GenfitMsg::get()<<MSG::ERROR<<"silicon addSpacePointTrakerHit"
+ <<detID<<" "<<hit.getCellID()<<" faieled"<<endmsg;
+ }
+ }else if(detID==7){
+ //if(addSpacePointMeasurement(p,sigma,hit.getCellID(),hitID)){
+ // GenfitMsg::get()<<MSG::DEBUG<<"add DC space point"<<endmsg;
+ // hitID++;
+ //}else{
+ // GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurement"
+ // <<detID<<" faieled" <<endmsg;
+ //}
+ float minTime=FLT_MAX;
+ edm4hep::ConstSimTrackerHit minTimeSimHit;
+ //Select the SimTrakerHit with least time
+ for(int iSimHit=0;iSimHit<(int) assoHits->size();iSimHit++){
+ if(assoHits->at(iSimHit).getRec()==hit &&
+ assoHits->at(iSimHit).getSim().getTime()<minTime){
+ minTimeSimHit=assoHits->at(iSimHit).getSim();
+ minTime=assoHits->at(iSimHit).getSim().getTime();
+ }
+ }
+ //std::cout<<"minTimeSimHit "<<minTimeSimHit<<std::endl;
+ sortedDCTrackHitCol.push_back(minTimeSimHit);
+ }else{
+ GenfitMsg::get()<<MSG::DEBUG<<" Skip add this hit!"<<endmsg;
+ }
+ }//end loop over hit on track
+
+ GenfitMsg::get()<<MSG::DEBUG<<" addSimTrakerHits trackerHits_size="
+ <<track.trackerHits_size()<<endmsg;
+
+ ///Add DC hits to track
+ //Sort sim DC hits by time
+ //std::sort(sortedDCTrackHitCol.begin(),sortedDCTrackHitCol.end(),sortDCHit);
+ for(auto dCTrackerHit: sortedDCTrackHitCol){
+ edm4hep::Vector3d pos=dCTrackerHit.getPosition();
+ TVectorD p(3);
+ p[0]=pos.x;
+ p[1]=pos.y;
+ p[2]=pos.z;
+ unsigned long long detID = dCTrackerHit.getCellID();
+ if(addSpacePointMeasurement(p,sigma,detID,hitID,smear)){
+ GenfitMsg::get()<<MSG::DEBUG<<"add DC space point"<<endmsg;
+ hitID++;
+ }else{
+ GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurement"
+ <<detID<<" faieled" <<endmsg;
+ }
+ }
+
+ GenfitMsg::get()<<MSG::DEBUG<<"GenfitTrack nHitAdded "<<hitID<<endmsg;
+ return hitID;
+}
+
+bool GenfitTrack::storeTrack(edm4hep::ReconstructedParticle& recParticle,
+ int pidType, int ndfCut, double chi2Cut)
+{
+
+ /// Get fit status
+ const genfit::FitStatus* fitState = m_track->getFitStatus();
+ double ndf = fitState->getNdf();
+ if(ndf>ndfCut) return false;
+ double chi2 = fitState->getChi2();
+ if(chi2>chi2Cut) return false;
+ double charge = fitState->getCharge();
+ int isFitted = fitState->isFitted();
+ int isConverged = fitState->isFitConverged();
+ int isConvergedFully = fitState->isFitConvergedFully();
+ TMatrixDSym fittedCov;
+ TLorentzVector fittedPos;
+ TVector3 fittedMom;
+ int fittedState=getFittedState(fittedPos,fittedMom,fittedCov);
+ GenfitMsg::get()<<MSG::DEBUG<<"fit result: get status OK? pidType "
+ <<pidType<<" fittedState==0 " <<(0==fittedState)<<" isFitted "<<isFitted
+ <<" isConverged "<<isConverged<<" ndf "<<ndf<<endmsg;
+ if((0!=fittedState) || (!isFitted) || !isConvergedFully || (ndf<ndfCut)){
+ GenfitMsg::get()<<MSG::DEBUG<<" fitting failed"<<endmsg;
+ }else{
+ GenfitMsg::get()<<MSG::DEBUG<<"fit result: Pos("<<
+ fittedPos.X()<<" "<<
+ fittedPos.Y()<<" "<<
+ fittedPos.Z()<<") mom("<<
+ fittedMom.X()<<" "<<
+ fittedMom.Y()<<" "<<
+ fittedMom.Z()<<") p_tot "<<
+ fittedMom.Mag()<<" pt "<<
+ fittedMom.Perp()<<
+ " chi2 "<<chi2<<
+ " ndf "<<ndf
+ <<endmsg;
+ }
+ float pos[3]={float(fittedPos.X()),float(fittedPos.Y()),float(fittedPos.Z())};
+ float mom[3]={float(fittedMom.X()),float(fittedMom.Y()),float(fittedMom.Z())};
+ HelixClass helix;
+ helix.Initialize_VP(pos,mom,charge,m_genfitField->getBz(fittedPos.Vect()));
+
+
+ /////track status at POCA to origin
+ //TVector3 origin(pivot);
+ //TVector3 pocaToOrigin_pos(1e9*dd4hep::cm,1e9*dd4hep::cm,1e9*dd4hep::cm);
+ //TVector3 pocaToOrigin_mom(1e9*dd4hep::GeV,1e9*dd4hep::GeV,1e9*dd4hep::GeV);
+
+ //if(ExtrapolateToPoint(pocaToOrigin_pos,pocaToOrigin_mom,
+ // m_track,origin) > 1e6*dd4hep::cm){
+ // log<<"extrapolate to origin failed";
+ // return false;
+ //}
+
+
+ //new TrackState
+ edm4hep::TrackState* trackState = new edm4hep::TrackState();
+ trackState->location=0;//edm4hep::AtIP;
+ trackState->D0=helix.getD0();
+ trackState->phi=helix.getPhi0();
+ trackState->omega=helix.getOmega();
+ trackState->Z0=helix.getZ0();
+ trackState->tanLambda=helix.getTanLambda();
+ trackState->referencePoint=helix.getReferencePoint();
+
+ // std::array<float,15> covMatrix;
+ // int k=0;
+ // for(int i=0;i<5;i++){
+ // for(int j=0;j<5;j++){
+ // if(i<=j) covMatrix[k]=;//FIXME
+ // }
+ // }
+ // trackState.covMatrix=
+
+ //new Track
+ edm4hep::Track* track = new edm4hep::Track();
+ //track->setType();
+ track->setChi2(fitState->getChi2());
+ track->setNdf(fitState->getNdf());
+ //track->setDEdx();
+ //track->setRadiusOfInnermostHit();//FIXME
+ //track->addToTrackerHits();
+
+ //new ReconstructedParticle
+ //recParticle->setType();
+ //dcRecParticle->setEnergy();
+
+ edm4hep::Vector3f momVec3(helix.getMomentum()[0],
+ helix.getMomentum()[1],helix.getMomentum()[2]);
+ recParticle.setMomentum(momVec3);
+ //recParticle->setReferencePoint(referencePoint);
+ recParticle.setCharge(helix.getCharge());
+ // recParticle->setMass();
+ // recParticle->setCovMatrix();
+ // rcRecParticle->setStartVertex();
+ //recParticle->addToTracks(track);
+
+ return true;
+}
+
+void GenfitTrack::pivotToFirstLayer(edm4hep::Vector3d& pos,
+ edm4hep::Vector3f& mom, edm4hep::Vector3d& firstPos,
+ edm4hep::Vector3f& firstMom)
+{
+ //FIXME, TODO
+ firstPos=pos;
+ firstMom=mom;
+}
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
new file mode 100644
index 0000000000000000000000000000000000000000..708c33e38ba282c92f5d0f8ef579fde688d150f4
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
@@ -0,0 +1,204 @@
+//////////////////////////////////////////////////////////////////
+///
+/// This is an interface of to call genfit
+/// A genfit track can be created, fitted and extrapolated
+/// Track with only track representation(s)
+///
+/// In this file, including:
+/// a genfit track class
+///
+/// Units are following DD4hepUnits
+///
+/// Authors:
+/// Zhang Yao (zhangyao@ihep.ac.cn)
+/// Y.Fujii (yfujii@ihep.ac.cn)
+/// Yohei Nakatsugawa (yohei@ihep.ac.cn)
+///
+//////////////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_GENFITTRACK_H
+#define RECGENFITALG_GENFITTRACK_H
+
+#include "GenfitFitter.h"
+
+//ROOT
+#include "TVectorD.h"
+#include "TMatrixDSym.h"
+
+//STL
+#include <vector>
+
+class TLorentzVector;
+
+namespace genfit{
+ class Track;
+ class FitStatus;
+ class AbsTrackRep;
+ class RKTrackRep;
+ class KalmanFittedStateOnPlane;
+}
+namespace edm4hep{
+ class MCParticle;
+ class SimTrackerHitCollection;
+ class ReconstructedParticle;
+ class MCRecoTrackerAssociationCollection;
+ class Track;
+ class ConstTrackerHit;
+ class Vector3d;
+ class Vector3f;
+}
+namespace dd4hep {
+ namespace DDSegmentation{
+ class GridDriftChamber;
+ }
+}
+
+class GenfitTrack {
+ friend int GenfitFitter::processTrack(
+ GenfitTrack* track, bool resort=false);
+
+ friend int GenfitFitter::processTrackWithRep(
+ GenfitTrack* track, int repID=0, bool resort=false);
+
+ friend double GenfitFitter::extrapolateToHit(TVector3& poca,
+ TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, const GenfitTrack* track,
+ TVector3 pos, TVector3 mom, TVector3 end0, TVector3 end1, int debug,
+ int repID=0, bool stopAtBoundary=false, bool calcJacobianNoise=false);
+
+ friend double GenfitFitter::extrapolateToCylinder(TVector3& pos,
+ TVector3& mom,
+ GenfitTrack* track, double radius, const TVector3 linePoint,
+ const TVector3 lineDirection, int hitID =0, int repID=0,
+ bool stopAtBoundary=false, bool calcJacobianNoise=false);
+
+ friend double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom,
+ const GenfitTrack* genfitTrack, const TVector3& point, int repID=0,
+ bool stopAtBoundary = false, bool calcJacobianNoise = false) const;
+
+ public:
+ GenfitTrack(const GenfitField* field,
+ const dd4hep::DDSegmentation::GridDriftChamber* seg,
+ const char* name="GenfitTrack");
+ virtual ~GenfitTrack();
+
+ /// Add a Genfit track
+ virtual bool createGenfitTrack(int pdgType,int charge,TLorentzVector pos, TVector3 mom,
+ TMatrixDSym covM);
+ //virtual bool createGenfitTrack(TLorentzVector posInit,TVector3 momInit,
+ //TMatrixDSym covMInit);
+
+ ///Create genfit track from MCParticle
+ bool createGenfitTrackFromMCParticle(int pidTyep,const edm4hep::MCParticle&
+ mcParticle, double eventStartTime=0.);
+ bool createGenfitTrackFromEDM4HepTrack(int pidType,const edm4hep::Track& track,
+ double eventStartTime);
+
+ // /// Prepare a hit list, return number of hits on track
+ // int PrepareHits();//TODO
+
+ /// Add a space point measurement, return number of hits on track
+ bool addSpacePointTrakerHit(edm4hep::ConstTrackerHit& hit, int hitID);
+
+ /// Add a space point measurement, return number of hits on track
+ virtual bool addSpacePointMeasurement(const TVectorD&, double,
+ int detID=-1, int hitID=-1, bool smear=false);
+
+ /// Add a WireMeasurement with MC truth position smeared by sigma
+ virtual void addWireMeasurement(double driftDistance,
+ double driftDistanceError, const TVector3& endPoint1,
+ const TVector3& endPoint2, int lrAmbig, int detID, int hitID);
+
+ /// Add a WireMeasurement with DC digi
+ virtual bool addWireMeasurementOnTrack(edm4hep::Track& track, double sigma);
+
+ ///Add space point from truth to track
+ int addSimTrackerHits(const edm4hep::Track& track,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ float sigma,bool smear=false);// float nSigmaSelection
+
+ ///Store track to ReconstructedParticle
+ bool storeTrack(edm4hep::ReconstructedParticle& dcRecParticle,int pidType,
+ int ndfCut=1e9, double chi2Cut=1.e9);
+
+ ///A tool to convert track to the first layer of DC
+ void pivotToFirstLayer(edm4hep::Vector3d& pos,edm4hep::Vector3f& mom,
+ edm4hep::Vector3d& firstPos, edm4hep::Vector3f& firstMom);
+
+ /// Copy a track to event
+ //void CopyATrack()const;
+
+ ///Extrapolate to Hit
+ double extrapolateToHit( TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, TVector3 pos, TVector3 mom,
+ TVector3 end0,//one end of the hit wire
+ TVector3 end1,//the orhter end of the hit wire
+ int debug,
+ int repID,
+ bool stopAtBoundary,
+ bool calcJacobianNoise) const;
+
+ bool getPosMomCovMOP(int hitID, TLorentzVector& pos, TVector3& mom,
+ TMatrixDSym& cov, int repID=0) const;
+
+ /// get the seed position and momentum from track
+ const TLorentzVector getSeedStatePos() const;
+ const TVector3 getSeedStateMom() const;
+ void getSeedStateMom(TLorentzVector& pos, TVector3& mom) const;
+ unsigned int getNumPoints() const;
+
+ /// get the fitted track status
+ const genfit::FitStatus* getFitStatus(int repID=0) const;
+ int getFittedState(TLorentzVector& pos, TVector3& mom, TMatrixDSym& cov,
+ int repID=0, bool biased=true) const;
+ int getNumPointsWithFittedInfo(int repID=0) const;
+ bool getFirstPointWithFittedInfo(int repID=0) const;
+ bool fitSuccess(int repID) const;
+
+ /// get the wire infomation
+ int getDetIDWithFitterInfo(int hitID, int idRaw=0) const;
+
+ int getPDG(int id=0) const;
+ int getPDGCharge(int id=0) const;
+
+ /// print genfit track
+ void printSeed() const;//print seed
+ void printFitted(int repID=0) const;//print seed
+ void print(TLorentzVector pos, TVector3 mom, const char* comment="") const;
+
+ /// set and get debug level
+ void setDebug(int debug);
+ int getDebug(void) const { return m_debug;}
+
+ /// get name of this object
+ const char* getName() const {return m_name;}
+ genfit::Track* getTrack() const{return m_track;}
+
+ /// Add a track representation
+ genfit::RKTrackRep* addTrackRep(int pdg);
+
+ protected:
+ //genfit::Track* getTrack() {return m_track;}
+
+ private:
+
+ const char* m_name;
+
+ ///Note! private functions are using genfit unit, cm and MeV
+
+ genfit::AbsTrackRep* getRep(int id=0) const {return m_reps[id];}
+ bool getMOP(int hitID, genfit::MeasuredStateOnPlane& mop,
+ genfit::AbsTrackRep* trackRep=nullptr) const;
+
+ genfit::Track* m_track;/// track
+ std::vector<genfit::AbsTrackRep*> m_reps;/// track repesentations
+ int m_debug;/// debug level
+
+ const GenfitField* m_genfitField;//pointer to genfit field
+ const dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber;
+
+ static const int s_PDG[2][5];
+
+};
+
+#endif
diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cef756adc1d6cc624bf8519d69b5e367423413d9
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
@@ -0,0 +1,409 @@
+#include "RecGenfitAlgDC.h"
+#include "GenfitTrack.h"
+#include "GenfitFitter.h"
+#include "GenfitField.h"
+
+//genfit
+#include "EventDisplay.h"
+
+//cepcsw
+#include "DetInterface/IGeomSvc.h"
+#include "DataHelper/HelixClass.h"
+#include "DetSegmentation/GridDriftChamber.h"
+
+//externals
+#include "edm4hep/EventHeaderCollection.h"
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/ReconstructedParticle.h"
+#include "edm4hep/ReconstructedParticleCollection.h"
+#include "edm4hep/Track.h"
+#include "edm4hep/TrackCollection.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/Detector.h"
+#include "UTIL/BitField64.h"
+#include "DDSegmentation/Segmentation.h"
+#include "TRandom.h"
+#include "TLorentzVector.h"
+
+//stl
+#include "time.h"
+
+
+DECLARE_COMPONENT( RecGenfitAlgDC )
+
+/////////////////////////////////////////////////////////////////////
+RecGenfitAlgDC::RecGenfitAlgDC(const std::string& name, ISvcLocator* pSvcLocator):
+ GaudiAlgorithm(name, pSvcLocator),m_nPDG(5),m_dd4hep(nullptr),
+ m_gridDriftChamber(nullptr),m_decoder(nullptr)
+{
+ //declareProperty("EventHeaderCol", _headerCol);
+ declareProperty("MCParticle", m_mcParticleCol,
+ "Handle of the input MCParticle collection");
+ declareProperty("DriftChamberHitsCollection", m_simDCHitCol,
+ "Handle of the input SimHit collection");
+ declareProperty("DigiDCHitCollection", m_digiDCHitsCol,
+ "Handle of digi DCHit collection");
+ declareProperty("DCTrackCollection", m_dcTrackCol,
+ "Handle of DC track collection");
+ declareProperty("DCHitAssociationCollection", m_dcHitAssociationCol,
+ "Handle of simTrackerHit and TrackerHit association collection");
+ declareProperty("DCRecParticleCollection", m_dcRecParticleCol,
+ "Handle of drift chamber reconstructed particle collection");
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode RecGenfitAlgDC::initialize()
+{
+ MsgStream log(msgSvc(), name());
+ info()<<"RecGenfitAlgDC initialize()"<<endmsg;
+
+ ///Get GeomSvc
+ m_geomSvc=Gaudi::svcLocator()->service("GeomSvc");
+ if (nullptr==m_geomSvc) {
+ std::cout<<"Failed to find GeomSvc"<<std::endl;
+ return StatusCode::FAILURE;
+ }
+ ///Get Detector
+ m_dd4hep = m_geomSvc->lcdd();
+ ///Get Field
+ m_dd4hepField=m_geomSvc->lcdd()->field();
+
+ /// New a genfit fitter
+ m_genfitFitter=new GenfitFitter(m_fitterType.toString().c_str());
+ m_genfitField=new GenfitField(m_dd4hepField);
+ m_genfitFitter->setField(m_genfitField);
+ m_genfitFitter->setGeoMaterial(m_geomSvc->lcdd());
+ m_genfitFitter->setEnergyLossBrems(m_correctBremsstrahlung);
+ m_genfitFitter->setNoiseBrems(m_correctBremsstrahlung);
+ if(m_debug>10) m_genfitFitter->setDebug(m_debug-10);
+ if(m_noMaterialEffects) m_genfitFitter->setNoEffects(true);
+ if(-1==m_debugPid) m_genfitFitter->setNoEffects(true);
+ if(-1==m_debugPid) m_debugPid=0;//charged geantino with electron pid
+ if(m_fitterType=="DAF"||m_fitterType=="DafRef"){
+ m_genfitFitter->setMaxIterationsBetas(m_bStart,m_bFinal,m_maxIteration);
+ } else {
+ m_genfitFitter->setMaxIterations(m_maxIteration);
+ }
+ //print genfit parameters
+ if(m_debug) m_genfitFitter->print();
+ if(""!=m_genfitHistRootName) m_genfitFitter->initHist(m_genfitHistRootName);
+
+ //initialize member vairables
+ for(int i=0;i<5;i++) m_fitSuccess[i]=0;
+ m_nDCTrack=0;
+ ///Get Readout
+ dd4hep::Readout readout=m_dd4hep->readout(m_readout_name);
+ ///Get Segmentation
+ m_gridDriftChamber=dynamic_cast<dd4hep::DDSegmentation::GridDriftChamber*>
+ (readout.segmentation().segmentation());
+ if(nullptr==m_gridDriftChamber){
+ error() << "Failed to get the GridDriftChamber" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ ///Get Decoder
+ m_decoder = m_geomSvc->getDecoder(m_readout_name);
+ if (nullptr==m_decoder) {
+ error() << "Failed to get the decoder" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+
+ ///book tuple
+ NTuplePtr nt(ntupleSvc(), "RecGenfitAlgDC/recGenfitAlgDC");
+ if(nt){
+ m_tuple=nt;
+ }else{
+ m_tuple=ntupleSvc()->book("RecGenfitAlgDC/recGenfitAlgDC",
+ CLID_ColumnWiseTuple,"RecGenfitAlgDC");
+ if(m_tuple){
+ StatusCode sc;
+ sc=m_tuple->addItem("run",m_run);
+ sc=m_tuple->addItem("evt",m_evt);
+ sc=m_tuple->addItem("tkId",m_tkId);
+ sc=m_tuple->addItem("mcIndex",m_mcIndex,0,100);//max. 100 particles
+ sc=m_tuple->addItem("truthPocaMc",m_mcIndex,m_truthPocaMc,3);
+ sc=m_tuple->addItem("pocaPosMc",m_mcIndex,m_pocaPosMc,3);
+ sc=m_tuple->addItem("pocaMomMc",m_mcIndex,m_pocaMomMc,3);
+ sc=m_tuple->addItem("pocaMomMcP",m_mcIndex,m_pocaMomMcP);
+ sc=m_tuple->addItem("pocaMomMcPt",m_mcIndex,m_pocaMomMcPt);
+ sc=m_tuple->addItem("pocaPosMdc",3,m_pocaPosMdc);
+ sc=m_tuple->addItem("pocaMomMdc",3,m_pocaMomMdc);
+ sc=m_tuple->addItem("index",m_pidIndex, 0, 5);
+ sc=m_tuple->addItem("firstPosKalP",5,3,m_firstPosKal);
+ sc=m_tuple->addItem("firstMomKalP",5,m_firstMomKalP);
+ sc=m_tuple->addItem("firstMomKalPt",5,m_firstMomKalPt);
+ sc=m_tuple->addItem("pocaPosKal",5,3,m_pocaPosKal);
+ sc=m_tuple->addItem("pocaMomKal",5,3,m_pocaMomKal);
+ sc=m_tuple->addItem("pocaMomKalP",5,m_pocaMomKalP);
+ sc=m_tuple->addItem("pocaMomKalPt",5,m_pocaMomKalPt);
+ sc=m_tuple->addItem("chargeKal",5,m_chargeKal);
+ sc=m_tuple->addItem("nDofKal",5,m_nDofKal);
+ sc=m_tuple->addItem("chi2Kal",5,m_chi2Kal);
+ sc=m_tuple->addItem("isFitted",5,m_isFitted);
+ sc=m_tuple->addItem("isFitConverged",5,m_isFitConverged);
+ sc=m_tuple->addItem("isFitConvergedFully",5,
+ m_isFitConvergedFully);
+ sc=m_tuple->addItem("nHitFailedKal",5,m_nHitFailedKal);
+ sc=m_tuple->addItem("nHitFitted",5,m_nHitFitted);
+ sc=m_tuple->addItem("nDigi",m_nDigi);
+ sc=m_tuple->addItem("nHitMc",m_nHitMc);
+ sc=m_tuple->addItem("nHitKalInput",m_nHitKalInput);
+ sc=m_tuple->addItem("nHitWithFitInfo",5,m_nHitWithFitInfo);
+ sc=m_tuple->addItem("nSimDCHit",m_nSimDCHit,0,50000);
+ sc=m_tuple->addItem("mdcHitDriftT",m_nSimDCHit,m_mdcHitDriftT);
+ sc=m_tuple->addItem("mdcHitDriftDl",m_nSimDCHit,m_mdcHitDriftDl);
+ sc=m_tuple->addItem("mdcHitDriftDr",m_nSimDCHit,m_mdcHitDriftDr);
+ sc=m_tuple->addItem("mdcHitLr",m_nSimDCHit,m_mdcHitLr);
+ sc=m_tuple->addItem("mdcHitLayer",m_nSimDCHit,m_mdcHitLayer);
+ sc=m_tuple->addItem("mdcHitWire",m_nSimDCHit,m_mdcHitWire);
+ sc=m_tuple->addItem("mdcHitExpDoca",m_nSimDCHit,m_mdcHitExpDoca);
+ sc=m_tuple->addItem("mdcHitExpMcDoca",m_nSimDCHit,m_mdcHitExpMcDoca);
+ sc=m_tuple->addItem("mdcHitErr",m_nSimDCHit,m_mdcHitErr);
+ sc=m_tuple->addItem("time",5,m_time);
+ sc=m_tuple->addItem("mdcHitMcTkId",m_nSimDCHit,m_mdcHitMcTkId);
+ sc=m_tuple->addItem("mdcHitMcLr",m_nSimDCHit,m_mdcHitMcLr);
+ sc=m_tuple->addItem("mdcHitMcDrift",m_nSimDCHit,m_mdcHitMcDrift);
+ sc=m_tuple->addItem("mdcHitMcX",m_nSimDCHit,m_mdcHitMcX);
+ sc=m_tuple->addItem("mdcHitMcY",m_nSimDCHit,m_mdcHitMcY);
+ sc=m_tuple->addItem("mdcHitMcZ",m_nSimDCHit,m_mdcHitMcZ);
+ sc=m_tuple->addItem("mcPocaX",m_nSimDCHit,m_mdcHitExpMcPocaX);
+ sc=m_tuple->addItem("mcPocaY",m_nSimDCHit,m_mdcHitExpMcPocaY);
+ sc=m_tuple->addItem("mcPocaZ",m_nSimDCHit,m_mdcHitExpMcPocaZ);
+ sc=m_tuple->addItem("mcPocaWireX",m_nSimDCHit,m_mdcHitExpMcPocaWireX);
+ sc=m_tuple->addItem("mcPocaWireY",m_nSimDCHit,m_mdcHitExpMcPocaWireY);
+ sc=m_tuple->addItem("mcPocaWireZ",m_nSimDCHit,m_mdcHitExpMcPocaWireZ);
+ debug()<< "Book tuple RecGenfitAlgDC/genfit" << endmsg;
+ }else{
+ error()<< "Cannot book tuple RecGenfitAlgDC/genfit" << endmsg;
+ }
+ }//end of book tuple
+
+ //init genfit event display
+ //if(m_showDisplay) m_genfitDisplay = genfit::EventDisplay::getInstance();
+
+ return StatusCode::SUCCESS;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode RecGenfitAlgDC::execute()
+{
+ StatusCode sc;
+ m_timer=clock();
+ info()<<"RecGenfitAlgDC in execute()"<<endmsg;
+
+ /////retrieve EventHeader
+ //auto header = _headerCol.get()->at(0);
+ //int evtNo = header.getEventNumber();
+ //int runNo = header.getRunNumber();
+
+ //std::cout<<"run "<<header.getEventNumber()
+ // <<" "<<header.getRunNumber()<<std::endl;
+
+ ///retrieve Track and TrackHits
+ const edm4hep::TrackCollection* dcTrackCol=nullptr;
+ dcTrackCol=m_dcTrackCol.get();
+ if(nullptr==dcTrackCol) {
+ debug()<<"TrackCollection not found"<<endmsg;
+ return StatusCode::SUCCESS;
+ }
+ const edm4hep::TrackerHitCollection* didiDCHitsCol=nullptr;
+ didiDCHitsCol=m_digiDCHitsCol.get();
+ if(nullptr==didiDCHitsCol) {
+ debug()<<"DigiDCHitCollection not found"<<endmsg;
+ return StatusCode::SUCCESS;
+ }
+ ///retrieve DC Hit Association
+ auto assoDCHitsCol=m_dcHitAssociationCol.get();
+
+ edm4hep::ReconstructedParticleCollection* dcRecParticleCol=
+ m_dcRecParticleCol.createAndPut();
+
+ ///----------------------------------------------------
+ ///Loop over Track and do fitting for each track
+ ///----------------------------------------------------
+ debug()<<"DCTrackCol size="<<dcTrackCol->size()<<endmsg;
+ for(auto dcTrack: *dcTrackCol){
+ ///Loop over 5 particle hypothesis(0-4): e,mu,pi,K,p
+ ///-1 for chargedgeantino
+ for(unsigned int pidType=0;pidType<m_nPDG;pidType++){
+ if((m_debugPid>=0) && (m_debugPid!=pidType)) continue;
+ ///-----------------------------------
+ ///Create a GenFit track
+ ///-----------------------------------
+ GenfitTrack* genfitTrack=new GenfitTrack(m_genfitField,
+ m_gridDriftChamber);
+ genfitTrack->setDebug(m_debug);
+ double eventStartTime=0;
+ if(!genfitTrack->createGenfitTrackFromEDM4HepTrack(pidType,dcTrack,
+ eventStartTime)){
+ debug()<<"createGenfitTrackFromEDM4HepTrack failed!"<<endmsg;
+ return StatusCode::SUCCESS;
+ }
+ if(m_useTruthHit){
+ if(0==genfitTrack->addSimTrackerHits(dcTrack,assoDCHitsCol,
+ m_sigmaHit.value(),m_smearHit)){
+ debug()<<"addSimTrackerHits failed!"<<endmsg;
+ return StatusCode::FAILURE;
+ }
+ }else{
+ if(0==genfitTrack->addWireMeasurementOnTrack(dcTrack,
+ m_sigmaHit.value())){
+ debug()<<"addWireMeasurementOnTrack failed!"<<endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ if(m_debug) genfitTrack->printSeed();
+
+ ///-----------------------------------
+ ///call genfit fitting procedure
+ ///-----------------------------------
+ m_genfitFitter->processTrack(genfitTrack,m_resortHits);
+ m_genfitFitter->setDebug(m_debug);
+
+ ///-----------------------------------
+ ///Store track
+ ///-----------------------------------
+ auto dcRecParticle=dcRecParticleCol->create();
+ genfitTrack->storeTrack(dcRecParticle,pidType,m_ndfCut,
+ m_chi2Cut);
+ if(m_debug) genfitTrack->printSeed();
+
+ if(m_tuple) debugTrack(pidType,genfitTrack);
+ if(m_showDisplay) {
+ //m_genfitDisplay->addEvent(genfitTrack->getTrack());
+ //m_genfitDisplay->open();
+ }else{
+ delete genfitTrack;
+ }
+ ++m_fitSuccess[pidType];
+ }//end loop over particle type
+ }//end loop over a track
+
+ if(m_tuple) debugEvent();
+
+
+ //if(m_genfitDisplay) while(1){
+ // std::cout<<"Press any key to finish..."<<std::endl;
+ // //system ("pause");
+ //}
+
+ if(m_tuple) sc=m_tuple->write();
+
+ return StatusCode::SUCCESS;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode RecGenfitAlgDC::finalize()
+{
+ MsgStream log(msgSvc(), name());
+ info()<< "RecGenfitAlgDC in finalize()" << endmsg;
+
+ m_genfitFitter->writeHist();
+ delete m_genfitFitter;
+ info()<<"RecGenfitAlgDC nRecTrack="<<m_nDCTrack<<" success e "
+ <<m_fitSuccess[0]<<" mu "<<m_fitSuccess[1]<<" pi "<<m_fitSuccess[2]
+ <<" K "<<m_fitSuccess[3]<<" p "<<m_fitSuccess[4]<<std::endl;
+ if(m_nDCTrack>0){
+ std::cout<<"RecGenfitAlgDC Success rate = "<<std::endl;
+ for (int i=0;i<5;i++){
+ std::cout<<Form("%d %2.2f",i,((float) m_fitSuccess[i])/m_nDCTrack)
+ <<std::endl;
+ }
+ }
+ return StatusCode::SUCCESS;
+}
+
+void RecGenfitAlgDC::debugTrack(int pidType,const GenfitTrack* genfitTrack)
+{
+ /// Get fit status
+ const genfit::FitStatus* fitState = genfitTrack->getFitStatus();
+ int charge= fitState->getCharge();
+
+ m_chargeKal[pidType]= charge;
+ m_nHitWithFitInfo[pidType]=genfitTrack->getNumPointsWithFittedInfo();
+ m_chi2Kal[pidType]=fitState->getChi2();
+ m_nDofKal[pidType]=fitState->getNdf();
+ m_isFitted[pidType]=(int)fitState->isFitted();
+ m_isFitConverged[pidType]=(int) fitState->isFitConverged();
+ m_isFitConvergedFully[pidType]=(int) fitState->isFitConvergedFully();
+
+ ///get fitted state of track
+ TMatrixDSym fittedCov;
+ TLorentzVector fittedPos;
+ TVector3 fittedMom;
+ int fittedState=genfitTrack->getFittedState(fittedPos,fittedMom,fittedCov);
+ HelixClass helix;//mm and GeV
+ float pos[3]={float(fittedPos.X()/dd4hep::mm),float(fittedPos.Y()/dd4hep::mm),
+ float(fittedPos.Z()/dd4hep::mm)};
+ float mom[3]={float(fittedMom.X()),float(fittedMom.Y()),float(fittedMom.Z())};
+ helix.Initialize_VP(pos,mom,charge,m_genfitField->getBz(fittedPos.Vect()));
+ m_pocaMomKalP[pidType]=fittedMom.Mag();
+
+ if(m_debug>0){
+ /// Get fit status
+ debug()<<"evt "<<m_evt<<" fit result: get status OK? pidType "
+ <<pidType<<" fittedState "<<fittedState<<" isFitted "
+ <<m_isFitted[pidType]<<" isConverged "<<m_isFitConverged[pidType]
+ <<" isFitConvergedFully "<<m_isFitConvergedFully[pidType]
+ <<" ndf "<<m_nDofKal[pidType]
+ <<" chi2 "<<m_chi2Kal[pidType]<<endmsg;
+ if((0!=fittedState)||(!m_isFitted[pidType])||(m_nDofKal[pidType]<m_ndfCut)){
+ debug()<<"fitting failed"<<endmsg;
+ }else{
+ debug()<<"evt "<<m_evt<<" fit result: Pos("<<
+ fittedPos.X()<<" "<<
+ fittedPos.Y()<<" "<<
+ fittedPos.Z()<<") mom("<<
+ fittedMom.X()<<" "<<
+ fittedMom.Y()<<" "<<
+ fittedMom.Z()<<") p_tot "<<
+ fittedMom.Mag()<<" pt "<<
+ fittedMom.Perp()<<endmsg;
+ }
+ }
+
+}
+
+void RecGenfitAlgDC::debugEvent()
+{
+ const edm4hep::MCParticleCollection* mcParticleCol = nullptr;
+ const edm4hep::SimTrackerHitCollection* simDCHitCol=nullptr;
+
+ m_pidIndex=5;
+
+ mcParticleCol=m_mcParticleCol.get();
+ simDCHitCol=m_simDCHitCol.get();
+ m_nSimDCHit=simDCHitCol->size();
+ int iMcParticle=0;
+ int iHit=0;
+ for(auto mcParticle : *mcParticleCol){
+ for(auto simDCHit: *simDCHitCol){
+ edm4hep::Vector3d pos=simDCHit.position();
+ TVectorD p(3);
+ p[0]=pos.x;//no unit conversion here
+ p[1]=pos.y;
+ p[2]=pos.z;
+ m_mdcHitMcX[iHit]=pos.x;
+ m_mdcHitMcY[iHit]=pos.y;
+ m_mdcHitMcZ[iHit]=pos.z;
+ iHit++;
+ }
+ edm4hep::Vector3f mcPocaMom = mcParticle.getMomentum();//GeV
+ float px=mcPocaMom.x;
+ float py=mcPocaMom.y;
+ float pz=mcPocaMom.z;
+ debug()<<" "<<px<<" "<<py<<" "<<pz<<endmsg;
+ m_pocaMomMcP[iMcParticle]=sqrt(px*px+py*py+pz*pz);
+ iMcParticle++;
+ }
+ m_mcIndex=iHit;
+
+}
diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h
new file mode 100644
index 0000000000000000000000000000000000000000..14922bc694d5c6909232a73795d9afbfd2338630
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h
@@ -0,0 +1,190 @@
+//////////////////////////////////////////////////////////////////////
+///
+/// This is an algorithm for track fitting for CEPC track with genfit.
+///
+/// In this file, including:
+/// An algorithm for drift chamber track fitting with genfit with 5 hypothesis
+///
+/// Units are following DD4hepUnits
+///
+/// Authors:
+/// Yao ZHANG(zhangyao@ihep.ac.cn)
+///
+/////////////////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_RECGENFITALGDC_H
+#define RECGENFITALG_RECGENFITALGDC_H
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/NTuple.h"
+#include "k4FWCore/DataHandle.h"
+#include "DD4hep/Fields.h"
+#include <string>
+
+class GenfitFitter;
+class GenfitField;
+class GenfitTrack;
+class IGeomSvc;
+class time;
+namespace genfit{
+ class EventDisplay;
+}
+namespace dd4hep {
+ class Detector;
+ //class rec::CellIDPositionConverter;
+ namespace DDSegmentation{
+ class GridDriftChamber;
+ class BitFieldCoder;
+ }
+}
+namespace edm4hep{
+ class EventHeaderCollection;
+ class MCParticleCollection;
+ class SimTrackerHitCollection;
+ class TrackCollection;
+ class TrackerHitCollection;
+ class MCRecoTrackerAssociationCollection;
+ class ReconstructedParticle;
+ class ReconstructedParticleCollection;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+class RecGenfitAlgDC:public GaudiAlgorithm {
+ public:
+ RecGenfitAlgDC (const std::string& name, ISvcLocator* pSvcLocator);
+ StatusCode initialize() override;
+ StatusCode execute() override;
+ StatusCode finalize() override;
+
+ private:
+ GenfitFitter* m_genfitFitter;//The pointer to a GenfitFitter
+ const GenfitField* m_genfitField;//The pointer to a GenfitField
+
+ void debugTrack(int pidType,const GenfitTrack* genfitTrack);
+ void debugEvent();
+
+ DataHandle<edm4hep::EventHeaderCollection> _headerCol{
+ "EventHeaderCol", Gaudi::DataHandle::Reader, this};
+ //Drift chamber rec hit and trac
+ DataHandle<edm4hep::TrackerHitCollection> m_digiDCHitsCol{
+ "DigiDCHitCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackCollection> m_dcTrackCol{
+ "DCTrackCollection", Gaudi::DataHandle::Reader, this};
+ //Mc truth
+ DataHandle<edm4hep::MCParticleCollection> m_mcParticleCol{
+ "MCParticle", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_simDCHitCol{
+ "DriftChamberHitsCollection" , Gaudi::DataHandle::Reader, this};
+
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection>
+ m_dcHitAssociationCol{"DCHitAssociationCollection",
+ Gaudi::DataHandle::Reader, this};
+ //output hits and particles
+ DataHandle<edm4hep::TrackerHitCollection> m_dcFitRecHitCol{
+ "DCFitRecHitsCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::ReconstructedParticleCollection> m_dcRecParticleCol{
+ "DCRecParticleCollection", Gaudi::DataHandle::Writer, this};
+
+ const unsigned int m_nPDG;//5:e,mu,pi,K,proton
+ SmartIF<IGeomSvc> m_geomSvc;
+ dd4hep::OverlayedField m_dd4hepField;
+ dd4hep::Detector* m_dd4hep;
+ dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+ Gaudi::Property<std::string> m_readout_name{this,
+ "readout", "DriftChamberHitsCollection"};
+ Gaudi::Property<int> m_debug{this,"debug",false};
+ Gaudi::Property<bool> m_smearHit{this,"smearHit",true};
+ Gaudi::Property<float> m_sigmaHit{this,"sigmaHit",0.11};//mm
+ Gaudi::Property<float> m_nSigmaHit{this,"nSigmaHit",5};
+ Gaudi::Property<double> m_initCovResPos{this,"initCovResPos",1};
+ Gaudi::Property<double> m_initCovResMom{this,"initCovResMom",0.1};
+ //Fitter type default is DAFRef.
+ //Candidates are DAF,DAFRef,KalmanFitter and KalmanFitterRefTrack.
+ Gaudi::Property<std::string> m_fitterType{this,"fitterTyep","DAFRef"};
+ Gaudi::Property<bool> m_correctBremsstrahlung{this,
+ "correctBremsstrahlung",false};
+ Gaudi::Property<bool> m_noMaterialEffects{this,
+ "noMaterialEffects",false};
+ Gaudi::Property<int> m_maxIteration{this,"maxIteration",20};
+ Gaudi::Property<int> m_resortHits{this,"resortHits",true};
+ Gaudi::Property<double> m_bStart{this,"bStart",100};
+ Gaudi::Property<double> m_bFinal{this,"bFinal",0.01};
+ Gaudi::Property<double> m_dcCornerCuts{this,"dcCornerCuts",-999};
+ Gaudi::Property<double> m_ndfCut{this,"ndfCut",1};
+ Gaudi::Property<double> m_chi2Cut{this,"chi2Cut",1000};
+ //-1,chargedGeantino;0,1,2,3,4:e,mu,pi,K,proton
+ Gaudi::Property<int> m_debugPid{this,"debugPid",-99};
+ Gaudi::Property<bool> m_useTruthTrack{this,"useTruthTrack",true};
+ Gaudi::Property<bool> m_useTruthHit{this,"useTruthHit",true};
+ Gaudi::Property<std::string> m_genfitHistRootName{this,
+ "genfitHistRootName",""};
+ Gaudi::Property<bool> m_showDisplay{this,"showDisplay",false};
+ int m_fitSuccess[5];
+ int m_nDCTrack;
+ //bool m_useRecLRAmbig;
+
+ genfit::EventDisplay* m_genfitDisplay;
+ clock_t m_timer;
+
+ /// tuples
+ NTuple::Tuple* m_tuple;
+ NTuple::Item<int> m_run;
+ NTuple::Item<int> m_evt;
+ NTuple::Item<int> m_tkId;
+ NTuple::Item<int> m_mcIndex;//number of navigated mcParicle
+ NTuple::Matrix<double> m_truthPocaMc;//2 dim matched particle and 3 pos.
+ NTuple::Matrix<double> m_pocaPosMc;//2 dim matched particle and 3 pos.
+ NTuple::Matrix<double> m_pocaMomMc;//2 dim matched particle and 3 mom.
+ NTuple::Array<double> m_pocaMomMcP;//2 dim matched particle and p
+ NTuple::Array<double> m_pocaMomMcPt;//2 dim matched particle and pt
+ NTuple::Array<double> m_pocaPosMdc;//pos 0:x,1:y,2:z
+ NTuple::Array<double> m_pocaMomMdc;//mom. 0:px,1:py,2:pz
+ NTuple::Item<int> m_pidIndex;
+ NTuple::Matrix<double> m_firstPosKal;//5 hyposis and pos. at first
+ NTuple::Array<double> m_firstMomKalP;//5 hyposis and mom. at first
+ NTuple::Array<double> m_firstMomKalPt;//5 hyposis and mom. at first
+ NTuple::Matrix<double> m_pocaPosKal;//5 hyposis and 3 mom.
+ NTuple::Matrix<double> m_pocaMomKal;//5 hyposis and 3 mom.
+ NTuple::Array<double> m_pocaMomKalP;//5 hyposis and p
+ NTuple::Array<double> m_pocaMomKalPt;//5 hyposis and pt
+ NTuple::Array<int> m_chargeKal;
+ NTuple::Array<double> m_chi2Kal;
+ NTuple::Array<double> m_nDofKal;
+ NTuple::Array<int> m_isFitConverged;
+ NTuple::Array<int> m_isFitConvergedFully;
+ NTuple::Array<int> m_isFitted;
+ NTuple::Item<int> m_nDigi;
+ NTuple::Item<int> m_nHitMc;
+ NTuple::Item<int> m_nSimDCHit;
+ NTuple::Array<int> m_nHitWithFitInfo;
+ NTuple::Item<int> m_nHitKalInput;
+ NTuple::Array<double> m_mdcHitDriftT;
+ NTuple::Array<double> m_mdcHitDriftDl;
+ NTuple::Array<double> m_mdcHitDriftDr;
+ NTuple::Array<int> m_mdcHitLr;
+ NTuple::Array<int> m_mdcHitLayer;
+ NTuple::Array<int> m_mdcHitWire;
+ NTuple::Array<double> m_mdcHitExpDoca;
+ NTuple::Array<double> m_mdcHitExpMcDoca;
+ NTuple::Array<double> m_mdcHitErr;
+ NTuple::Array<int> m_nHitFailedKal;
+ NTuple::Array<int> m_nHitFitted;
+ NTuple::Array<double> m_time;
+ //truth
+ NTuple::Array<int> m_mdcHitMcLr;
+ NTuple::Array<int> m_mdcHitMcTkId;
+ NTuple::Array<double> m_mdcHitMcDrift;
+ NTuple::Array<double> m_mdcHitMcX;
+ NTuple::Array<double> m_mdcHitMcY;
+ NTuple::Array<double> m_mdcHitMcZ;
+ NTuple::Array<double> m_mdcHitExpMcPocaX;
+ NTuple::Array<double> m_mdcHitExpMcPocaY;
+ NTuple::Array<double> m_mdcHitExpMcPocaZ;
+ NTuple::Array<double> m_mdcHitExpMcPocaWireX;
+ NTuple::Array<double> m_mdcHitExpMcPocaWireY;
+ NTuple::Array<double> m_mdcHitExpMcPocaWireZ;
+
+};
+#endif
diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..18f19c8eac116e44c259725e39951ab97b481bd4
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp
@@ -0,0 +1,404 @@
+#include "RecGenfitAlgSDT.h"
+#include "GenfitTrack.h"
+#include "GenfitFitter.h"
+#include "GenfitField.h"
+
+//genfit
+#include "EventDisplay.h"
+
+//cepcsw
+#include "DetInterface/IGeomSvc.h"
+#include "DataHelper/HelixClass.h"
+#include "DetSegmentation/GridDriftChamber.h"
+
+//externals
+#include "edm4hep/EventHeaderCollection.h"
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/ReconstructedParticle.h"
+#include "edm4hep/ReconstructedParticleCollection.h"
+#include "edm4hep/Track.h"
+#include "edm4hep/TrackCollection.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/Detector.h"
+#include "UTIL/BitField64.h"
+#include "DDSegmentation/Segmentation.h"
+#include "TRandom.h"
+#include "TLorentzVector.h"
+
+//stl
+#include "time.h"
+
+
+DECLARE_COMPONENT( RecGenfitAlgSDT )
+
+/////////////////////////////////////////////////////////////////////
+RecGenfitAlgSDT::RecGenfitAlgSDT(const std::string& name,
+ ISvcLocator* pSvcLocator):
+ GaudiAlgorithm(name, pSvcLocator),m_nPDG(5),m_dd4hep(nullptr),
+ m_gridDriftChamber(nullptr),m_decoder(nullptr)
+{
+ declareProperty("EventHeaderCollection", m_headerCol);
+ declareProperty("MCParticleCollection", m_mcParticleCol,
+ "Handle of the input MCParticle collection");
+ declareProperty("DriftChamberHitCollection", m_DCSimHitCol,
+ "Handle of the input DC SimTrakerHit collection");
+ declareProperty("DriftChamberDigiCollection", m_DCDigiCol,
+ "Handle of DC digi(TrakerHit) collection");
+ declareProperty("DCTrackCollection", m_DCTrackCol,
+ "Handle of DC track collection");
+ declareProperty("DCHitAssociationCollection", m_DCHitAssociationCol,
+ "Handle of simTrackerHit and TrackerHit association collection");
+ declareProperty("SDTTrackCollection", m_SDTTrackCol,
+ "Handle of input silicon track collection");
+ declareProperty("SDTRecParticleCollection", m_SDTRecParticleCol,
+ "Handle of silicon+drift chamber rec. particle collection");
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode RecGenfitAlgSDT::initialize()
+{
+ MsgStream log(msgSvc(), name());
+ info()<<" RecGenfitAlgSDT initialize()"<<endmsg;
+
+ ///Get GeomSvc
+ m_geomSvc=Gaudi::svcLocator()->service("GeomSvc");
+ if (nullptr==m_geomSvc) {
+ std::cout<<"Failed to find GeomSvc"<<std::endl;
+ return StatusCode::FAILURE;
+ }
+ ///Get Detector
+ m_dd4hep = m_geomSvc->lcdd();
+ ///Get Field
+ m_dd4hepField=m_geomSvc->lcdd()->field();
+
+ /// New a genfit fitter
+ m_genfitFitter=new GenfitFitter(m_fitterType.toString().c_str());
+ m_genfitField=new GenfitField(m_dd4hepField);
+ m_genfitFitter->setField(m_genfitField);
+ m_genfitFitter->setGeoMaterial(m_geomSvc->lcdd());
+ m_genfitFitter->setEnergyLossBrems(m_correctBremsstrahlung);
+ m_genfitFitter->setNoiseBrems(m_correctBremsstrahlung);
+ if(m_debug>10) m_genfitFitter->setDebug(m_debug-10);
+ if(m_noMaterialEffects) m_genfitFitter->setNoEffects(true);
+ if(-1==m_debugPid) m_genfitFitter->setNoEffects(true);
+ if(-1==m_debugPid) m_debugPid=0;//charged geantino with electron pid
+ if(m_fitterType=="DAF"||m_fitterType=="DafRef"){
+ m_genfitFitter->setMaxIterationsBetas(m_bStart,m_bFinal,m_maxIteration);
+ } else {
+ m_genfitFitter->setMaxIterations(m_maxIteration);
+ }
+ //print genfit parameters
+ if(m_debug) m_genfitFitter->print();
+ if(""!=m_genfitHistRootName) m_genfitFitter->initHist(m_genfitHistRootName);
+
+ //initialize member vairables
+ for(int i=0;i<5;i++) m_fitSuccess[i]=0;
+ m_nRecTrack=0;
+ ///Get Readout
+ dd4hep::Readout readout=m_dd4hep->readout(m_readout_name);
+ ///Get Segmentation
+ m_gridDriftChamber=dynamic_cast<dd4hep::DDSegmentation::GridDriftChamber*>
+ (readout.segmentation().segmentation());
+ if(nullptr==m_gridDriftChamber){
+ error() << "Failed to get the GridDriftChamber" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ ///Get Decoder
+ m_decoder = m_geomSvc->getDecoder(m_readout_name);
+ if (nullptr==m_decoder) {
+ error() << "Failed to get the decoder" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+
+ ///book tuple
+ NTuplePtr nt(ntupleSvc(), "RecGenfitAlgSDT/recGenfitAlgSDT");
+ if(nt){
+ m_tuple=nt;
+ }else{
+ m_tuple=ntupleSvc()->book("RecGenfitAlgSDT/recGenfitAlgSDT",
+ CLID_ColumnWiseTuple,"RecGenfitAlgSDT");
+ if(m_tuple){
+ StatusCode sc;
+ sc=m_tuple->addItem("run",m_run);
+ sc=m_tuple->addItem("evt",m_evt);
+ sc=m_tuple->addItem("tkId",m_tkId);
+ sc=m_tuple->addItem("mcIndex",m_mcIndex,0,100);//max. 100 particles
+ sc=m_tuple->addItem("truthPocaMc",m_mcIndex,m_truthPocaMc,3);
+ sc=m_tuple->addItem("pocaPosMc",m_mcIndex,m_pocaPosMc,3);
+ sc=m_tuple->addItem("pocaMomMc",m_mcIndex,m_pocaMomMc,3);
+ sc=m_tuple->addItem("pocaMomMcP",m_mcIndex,m_pocaMomMcP);
+ sc=m_tuple->addItem("pocaMomMcPt",m_mcIndex,m_pocaMomMcPt);
+ sc=m_tuple->addItem("pocaPosMdc",3,m_pocaPosMdc);
+ sc=m_tuple->addItem("pocaMomMdc",3,m_pocaMomMdc);
+ sc=m_tuple->addItem("index",m_pidIndex, 0, 5);
+ sc=m_tuple->addItem("firstPosKalP",5,3,m_firstPosKal);
+ sc=m_tuple->addItem("firstMomKalP",5,m_firstMomKalP);
+ sc=m_tuple->addItem("firstMomKalPt",5,m_firstMomKalPt);
+ sc=m_tuple->addItem("pocaPosKal",5,3,m_pocaPosKal);
+ sc=m_tuple->addItem("pocaMomKal",5,3,m_pocaMomKal);
+ sc=m_tuple->addItem("pocaMomKalP",5,m_pocaMomKalP);
+ sc=m_tuple->addItem("pocaMomKalPt",5,m_pocaMomKalPt);
+ sc=m_tuple->addItem("chargeKal",5,m_chargeKal);
+ sc=m_tuple->addItem("nDofKal",5,m_nDofKal);
+ sc=m_tuple->addItem("chi2Kal",5,m_chi2Kal);
+ sc=m_tuple->addItem("isFitted",5,m_isFitted);
+ sc=m_tuple->addItem("isFitConverged",5,m_isFitConverged);
+ sc=m_tuple->addItem("isFitConvergedFully",5,
+ m_isFitConvergedFully);
+ sc=m_tuple->addItem("nHitFailedKal",5,m_nHitFailedKal);
+ sc=m_tuple->addItem("nHitFitted",5,m_nHitFitted);
+ sc=m_tuple->addItem("nDigi",m_nDigi);
+ sc=m_tuple->addItem("nHitMc",m_nHitMc);
+ sc=m_tuple->addItem("nHitKalInput",m_nHitKalInput);
+ sc=m_tuple->addItem("nHitWithFitInfo",5,m_nHitWithFitInfo);
+ sc=m_tuple->addItem("nSimDCHit",m_nSimDCHit,0,50000);
+ sc=m_tuple->addItem("mdcHitDriftT",m_nSimDCHit,m_mdcHitDriftT);
+ sc=m_tuple->addItem("mdcHitDriftDl",m_nSimDCHit,m_mdcHitDriftDl);
+ sc=m_tuple->addItem("mdcHitDriftDr",m_nSimDCHit,m_mdcHitDriftDr);
+ sc=m_tuple->addItem("mdcHitLr",m_nSimDCHit,m_mdcHitLr);
+ sc=m_tuple->addItem("mdcHitLayer",m_nSimDCHit,m_mdcHitLayer);
+ sc=m_tuple->addItem("mdcHitWire",m_nSimDCHit,m_mdcHitWire);
+ sc=m_tuple->addItem("mdcHitExpDoca",m_nSimDCHit,m_mdcHitExpDoca);
+ sc=m_tuple->addItem("mdcHitExpMcDoca",m_nSimDCHit,m_mdcHitExpMcDoca);
+ sc=m_tuple->addItem("mdcHitErr",m_nSimDCHit,m_mdcHitErr);
+ sc=m_tuple->addItem("time",5,m_time);
+ sc=m_tuple->addItem("mdcHitMcTkId",m_nSimDCHit,m_mdcHitMcTkId);
+ sc=m_tuple->addItem("mdcHitMcLr",m_nSimDCHit,m_mdcHitMcLr);
+ sc=m_tuple->addItem("mdcHitMcDrift",m_nSimDCHit,m_mdcHitMcDrift);
+ sc=m_tuple->addItem("mdcHitMcX",m_nSimDCHit,m_mdcHitMcX);
+ sc=m_tuple->addItem("mdcHitMcY",m_nSimDCHit,m_mdcHitMcY);
+ sc=m_tuple->addItem("mdcHitMcZ",m_nSimDCHit,m_mdcHitMcZ);
+ sc=m_tuple->addItem("mcPocaX",m_nSimDCHit,m_mdcHitExpMcPocaX);
+ sc=m_tuple->addItem("mcPocaY",m_nSimDCHit,m_mdcHitExpMcPocaY);
+ sc=m_tuple->addItem("mcPocaZ",m_nSimDCHit,m_mdcHitExpMcPocaZ);
+ sc=m_tuple->addItem("mcPocaWireX",m_nSimDCHit,m_mdcHitExpMcPocaWireX);
+ sc=m_tuple->addItem("mcPocaWireY",m_nSimDCHit,m_mdcHitExpMcPocaWireY);
+ sc=m_tuple->addItem("mcPocaWireZ",m_nSimDCHit,m_mdcHitExpMcPocaWireZ);
+ debug()<< "Book tuple RecGenfitAlgSDT/genfit" << endmsg;
+ }else{
+ error()<< "Cannot book tuple RecGenfitAlgSDT/genfit" << endmsg;
+ }
+ }//end of book tuple
+
+ //init genfit event display
+ //if(m_showDisplay) m_genfitDisplay = genfit::EventDisplay::getInstance();
+
+ return StatusCode::SUCCESS;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode RecGenfitAlgSDT::execute()
+{
+ StatusCode sc=StatusCode::SUCCESS;
+ m_timer=clock();
+ info()<<"RecGenfitAlgSDT in execute()"<<endmsg;
+
+ /////retrieve EventHeader
+ //auto header = _headerCol.get()->at(0);
+ //int evtNo = header.getEventNumber();
+ //int runNo = header.getRunNumber();
+
+ //std::cout<<"run "<<header.getEventNumber()
+ // <<" "<<header.getRunNumber()<<std::endl;
+
+ ///retrieve silicon Track and TrackHits
+ const edm4hep::TrackCollection* sdtTrackCol=nullptr;
+ sdtTrackCol=m_SDTTrackCol.get();
+ if(nullptr==sdtTrackCol) {
+ debug()<<"TrackCollection not found"<<endmsg;
+ return StatusCode::SUCCESS;
+ }
+ //const edm4hep::TrackerHitCollection* dcDigiCol=nullptr;
+ //dcDigiCol=m_DCDigiCol.get();
+ //if(nullptr==dcDigiCol) {
+ // debug()<<"DigiDCHitCollection not found"<<endmsg;
+ // return StatusCode::SUCCESS;
+ //}
+ ///retrieve DC Hit Association
+ auto dcHitAssociationCol=m_DCHitAssociationCol.get();
+
+
+ edm4hep::ReconstructedParticleCollection* sdtRecParticleCol=
+ m_SDTRecParticleCol.createAndPut();
+
+ ///----------------------------------------------------
+ ///Loop over Track and do fitting for each track
+ ///----------------------------------------------------
+ debug()<<"SDTTrackCol size="<<sdtTrackCol->size()<<endmsg;
+ for(auto sdtTrack: *sdtTrackCol){
+ ///Loop over 5 particle hypothesis(0-4): e,mu,pi,K,p
+ ///-1 for chargedgeantino
+ for(unsigned int pidType=0;pidType<m_nPDG;pidType++){
+ if((m_debugPid>=0) && (m_debugPid!=pidType)) continue;
+ ///-----------------------------------
+ ///Create a GenFit track
+ ///-----------------------------------
+ GenfitTrack* genfitTrack=new GenfitTrack(m_genfitField,
+ m_gridDriftChamber);
+ genfitTrack->setDebug(m_debug);
+ double eventStartTime=0;
+ if(!genfitTrack->createGenfitTrackFromEDM4HepTrack(pidType,sdtTrack,
+ eventStartTime)){
+ debug()<<"createGenfitTrack failed!"<<endmsg;
+ return StatusCode::SUCCESS;
+ }
+ if(0==genfitTrack->addSimTrackerHits(sdtTrack,dcHitAssociationCol,
+ m_sigmaHit.value())){
+ debug()<<"addSimTrackerHits failed!"<<endmsg;
+ return StatusCode::SUCCESS;
+ }
+ if(m_debug) genfitTrack->printSeed();
+
+ ///-----------------------------------
+ ///call genfit fitting procedure
+ ///-----------------------------------
+ m_genfitFitter->processTrack(genfitTrack,m_resortHits);
+ m_genfitFitter->setDebug(m_debug);
+
+ ///-----------------------------------
+ ///Store track
+ ///-----------------------------------
+ auto dcRecParticle=sdtRecParticleCol->create();
+ genfitTrack->storeTrack(dcRecParticle,pidType,m_ndfCut,
+ m_chi2Cut);
+ if(m_debug) genfitTrack->printSeed();
+
+ if(m_tuple) debugTrack(pidType,genfitTrack);
+ if(m_showDisplay) {
+ //m_genfitDisplay->addEvent(genfitTrack->getTrack());
+ //m_genfitDisplay->open();
+ }else{
+ delete genfitTrack;
+ }
+ ++m_fitSuccess[pidType];
+ }//end loop over particle type
+ }//end loop over a track
+ m_nRecTrack++;
+
+ if(m_tuple) debugEvent();
+
+
+ //if(m_genfitDisplay) while(1){
+ // std::cout<<"Press any key to finish..."<<std::endl;
+ // //system ("pause");
+ //}
+
+ if(m_tuple) sc=m_tuple->write();
+
+ return StatusCode::SUCCESS;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode RecGenfitAlgSDT::finalize()
+{
+ MsgStream log(msgSvc(), name());
+ info()<< " RecGenfitAlgSDT in finalize()" << endmsg;
+
+ m_genfitFitter->writeHist();
+ delete m_genfitFitter;
+ info()<<"RecGenfitAlgSDT nRecTrack="<<m_nRecTrack<<" success e "
+ <<m_fitSuccess[0]<<" mu "<<m_fitSuccess[1]<<" pi "<<m_fitSuccess[2]
+ <<" K "<<m_fitSuccess[3]<<" p "<<m_fitSuccess[4]<<std::endl;
+ if(m_nRecTrack>0){
+ std::cout<<"RecGenfitAlgSDT Success rate = "<<std::endl;
+ for (int i=0;i<5;i++){
+ std::cout<<Form("%d %2.2f",i,((float) m_fitSuccess[i])/m_nRecTrack)
+ <<std::endl;
+ }
+ }
+ return StatusCode::SUCCESS;
+}
+
+void RecGenfitAlgSDT::debugTrack(int pidType,const GenfitTrack* genfitTrack)
+{
+ /// Get fit status
+ const genfit::FitStatus* fitState = genfitTrack->getFitStatus();
+ int charge= fitState->getCharge();
+
+ m_chargeKal[pidType]= charge;
+ m_nHitWithFitInfo[pidType]=genfitTrack->getNumPointsWithFittedInfo();
+ m_chi2Kal[pidType]=fitState->getChi2();
+ m_nDofKal[pidType]=fitState->getNdf();
+ m_isFitted[pidType]=(int)fitState->isFitted();
+ m_isFitConverged[pidType]=(int) fitState->isFitConverged();
+ m_isFitConvergedFully[pidType]=(int) fitState->isFitConvergedFully();
+
+ ///get fitted state of track
+ TMatrixDSym fittedCov;
+ TLorentzVector fittedPos;
+ TVector3 fittedMom;
+ int fittedState=genfitTrack->getFittedState(fittedPos,fittedMom,fittedCov);
+ HelixClass helix;//mm and GeV
+ float pos[3]={float(fittedPos.X()/dd4hep::mm),float(fittedPos.Y()/dd4hep::mm),
+ float(fittedPos.Z()/dd4hep::mm)};
+ float mom[3]={float(fittedMom.X()),float(fittedMom.Y()),float(fittedMom.Z())};
+ helix.Initialize_VP(pos,mom,charge,m_genfitField->getBz(fittedPos.Vect()));
+ m_pocaMomKalP[pidType]=fittedMom.Mag();
+
+ if(m_debug>0){
+ /// Get fit status
+ debug()<<"evt "<<m_evt<<" fit result: get status OK? pidType "
+ <<pidType<<" fittedState "<<fittedState<<" isFitted "
+ <<m_isFitted[pidType]<<" isConverged "<<m_isFitConverged[pidType]
+ <<" isFitConvergedFully "<<m_isFitConvergedFully[pidType]
+ <<" ndf "<<m_nDofKal[pidType]
+ <<" chi2 "<<m_chi2Kal[pidType]<<endmsg;
+ if((0!=fittedState)||(!m_isFitted[pidType])||(m_nDofKal[pidType]<m_ndfCut)){
+ debug()<<"fitting failed"<<endmsg;
+ }else{
+ debug()<<"evt "<<m_evt<<" fit result: Pos("<<
+ fittedPos.X()<<" "<<
+ fittedPos.Y()<<" "<<
+ fittedPos.Z()<<") mom("<<
+ fittedMom.X()<<" "<<
+ fittedMom.Y()<<" "<<
+ fittedMom.Z()<<") p_tot "<<
+ fittedMom.Mag()<<" pt "<<
+ fittedMom.Perp()<<endmsg;
+ }
+ }
+}
+
+void RecGenfitAlgSDT::debugEvent()
+{
+ const edm4hep::MCParticleCollection* mcParticleCol = nullptr;
+ const edm4hep::SimTrackerHitCollection* simDCHitCol=nullptr;
+
+ m_pidIndex=5;
+
+ mcParticleCol=m_mcParticleCol.get();
+ simDCHitCol=m_simDCHitCol.get();
+ m_nSimDCHit=simDCHitCol->size();
+ int iMcParticle=0;
+ int iHit=0;
+ for(auto mcParticle : *mcParticleCol){
+ for(auto simDCHit: *simDCHitCol){
+ edm4hep::Vector3d pos=simDCHit.position();
+ TVectorD p(3);
+ p[0]=pos.x;//no unit conversion here
+ p[1]=pos.y;
+ p[2]=pos.z;
+ m_mdcHitMcX[iHit]=pos.x;
+ m_mdcHitMcY[iHit]=pos.y;
+ m_mdcHitMcZ[iHit]=pos.z;
+ iHit++;
+ }
+ edm4hep::Vector3f mcPocaMom = mcParticle.getMomentum();//GeV
+ float px=mcPocaMom.x;
+ float py=mcPocaMom.y;
+ float pz=mcPocaMom.z;
+ debug()<<" "<<px<<" "<<py<<" "<<pz<<endmsg;
+ m_pocaMomMcP[iMcParticle]=sqrt(px*px+py*py+pz*pz);
+ iMcParticle++;
+ }
+ m_mcIndex=iHit;
+}
diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h
new file mode 100644
index 0000000000000000000000000000000000000000..1a1031e1b3db3b17d5c509ebbed76304f09a1600
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h
@@ -0,0 +1,192 @@
+//////////////////////////////////////////////////////////////////////
+///
+/// This is an algorithm for track fitting for CEPC track with genfit.
+///
+/// In this file, including:
+/// An algorithm for combined silicon and drift chamber track fitting
+/// with genfit for 5 particle hypothesis
+///
+/// Units are following DD4hepUnits
+///
+/// Authors:
+/// Yao ZHANG(zhangyao@ihep.ac.cn)
+///
+/////////////////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_RECGENFITALGSDT_H
+#define RECGENFITALG_RECGENFITALGSDT_H
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/NTuple.h"
+#include "k4FWCore/DataHandle.h"
+#include "DD4hep/Fields.h"
+#include <string>
+
+class GenfitFitter;
+class GenfitField;
+class GenfitTrack;
+class IGeomSvc;
+class time;
+namespace genfit{
+ class EventDisplay;
+}
+namespace dd4hep {
+ class Detector;
+ //class rec::CellIDPositionConverter;
+ namespace DDSegmentation{
+ class GridDriftChamber;
+ class BitFieldCoder;
+ }
+}
+namespace edm4hep{
+ class EventHeaderCollection;
+ class MCParticleCollection;
+ class SimTrackerHitCollection;
+ class TrackCollection;
+ class TrackerHitCollection;
+ class MCRecoTrackerAssociationCollection;
+ class ReconstructedParticle;
+ class ReconstructedParticleCollection;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+class RecGenfitAlgSDT:public GaudiAlgorithm {
+ public:
+ RecGenfitAlgSDT (const std::string& name, ISvcLocator* pSvcLocator);
+ StatusCode initialize() override; StatusCode execute() override; StatusCode finalize() override;
+
+ private:
+ GenfitFitter* m_genfitFitter;//The pointer to a GenfitFitter
+ const GenfitField* m_genfitField;//The pointer to a GenfitField
+
+ void debugTrack(int pidType,const GenfitTrack* genfitTrack);
+ void debugEvent();
+
+ DataHandle<edm4hep::EventHeaderCollection> m_headerCol{
+ "EventHeaderCol", Gaudi::DataHandle::Reader, this};
+ //Drift chamber rec hit and trac
+ DataHandle<edm4hep::TrackerHitCollection> m_DCDigiCol{
+ "DigiDCHitCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_DCSimHitCol{
+ "DCSimHitsCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackCollection> m_DCTrackCol{
+ "DCTrackCollection", Gaudi::DataHandle::Reader, this};
+ //Mc truth
+ DataHandle<edm4hep::MCParticleCollection> m_mcParticleCol{
+ "MCParticle", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_simDCHitCol{
+ "DriftChamberHitsCollection" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection>
+ m_DCHitAssociationCol{"DCHitAssociationCollection",
+ Gaudi::DataHandle::Reader, this};
+
+ //Track from silicon detectors
+ DataHandle<edm4hep::TrackCollection> m_SDTTrackCol{"SDTTrackCollection",
+ Gaudi::DataHandle::Reader, this};
+
+ //Output hits and particles
+ DataHandle<edm4hep::ReconstructedParticleCollection> m_SDTRecParticleCol{
+ "SDTRecParticleCollection", Gaudi::DataHandle::Writer, this};
+
+ const unsigned int m_nPDG;//5:e,mu,pi,K,proton
+ SmartIF<IGeomSvc> m_geomSvc;
+ dd4hep::OverlayedField m_dd4hepField;
+ dd4hep::Detector* m_dd4hep;
+ dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+ Gaudi::Property<std::string> m_readout_name{this,
+ "readout", "DriftChamberHitsCollection"};
+ Gaudi::Property<int> m_debug{this,"debug",0};
+ Gaudi::Property<float> m_sigmaHit{this,"sigmaHit",0.11};//mm
+ Gaudi::Property<float> m_nSigmaHit{this,"nSigmaHit",5};
+ Gaudi::Property<double> m_initCovResPos{this,"initCovResPos",1};
+ Gaudi::Property<double> m_initCovResMom{this,"initCovResMom",0.1};
+ //Fitter type default is DAFRef.
+ //Candidates are DAF,DAFRef,KalmanFitter and KalmanFitterRefTrack.
+ Gaudi::Property<std::string> m_fitterType{this,"fitterTyep","DAFRef"};
+ Gaudi::Property<bool> m_correctBremsstrahlung{this,
+ "correctBremsstrahlung",false};
+ Gaudi::Property<bool> m_noMaterialEffects{this,
+ "noMaterialEffects",false};
+ Gaudi::Property<int> m_maxIteration{this,"maxIteration",20};
+ Gaudi::Property<int> m_resortHits{this,"resortHits",true};
+ Gaudi::Property<double> m_bStart{this,"bStart",100};
+ Gaudi::Property<double> m_bFinal{this,"bFinal",0.01};
+ Gaudi::Property<double> m_DCCornerCuts{this,"dcCornerCuts",-999};
+ Gaudi::Property<double> m_ndfCut{this,"ndfCut",1};
+ Gaudi::Property<double> m_chi2Cut{this,"chi2Cut",1000};
+ //-1,chargedGeantino;0,1,2,3,4:e,mu,pi,K,proton
+ Gaudi::Property<int> m_debugPid{this,"debugPid",-99};
+ Gaudi::Property<bool> m_useTruthTrack{this,"useTruthTrack",true};
+ Gaudi::Property<bool> m_useTruthHit{this,"useTruthHit",true};
+ Gaudi::Property<std::string> m_genfitHistRootName{this,
+ "genfitHistRootName",""};
+ Gaudi::Property<bool> m_showDisplay{this,"showDisplay",false};
+ int m_fitSuccess[5];
+ int m_nRecTrack;
+ //bool m_useRecLRAmbig;
+
+ genfit::EventDisplay* m_genfitDisplay;
+ clock_t m_timer;
+
+ /// tuples
+ NTuple::Tuple* m_tuple;
+ NTuple::Item<int> m_run;
+ NTuple::Item<int> m_evt;
+ NTuple::Item<int> m_tkId;
+ NTuple::Item<int> m_mcIndex;//number of navigated mcParicle
+ NTuple::Matrix<double> m_truthPocaMc;//2 dim matched particle and 3 pos.
+ NTuple::Matrix<double> m_pocaPosMc;//2 dim matched particle and 3 pos.
+ NTuple::Matrix<double> m_pocaMomMc;//2 dim matched particle and 3 mom.
+ NTuple::Array<double> m_pocaMomMcP;//2 dim matched particle and p
+ NTuple::Array<double> m_pocaMomMcPt;//2 dim matched particle and pt
+ NTuple::Array<double> m_pocaPosMdc;//pos 0:x,1:y,2:z
+ NTuple::Array<double> m_pocaMomMdc;//mom. 0:px,1:py,2:pz
+ NTuple::Item<int> m_pidIndex;
+ NTuple::Matrix<double> m_firstPosKal;//5 hyposis and pos. at first
+ NTuple::Array<double> m_firstMomKalP;//5 hyposis and mom. at first
+ NTuple::Array<double> m_firstMomKalPt;//5 hyposis and mom. at first
+ NTuple::Matrix<double> m_pocaPosKal;//5 hyposis and 3 mom.
+ NTuple::Matrix<double> m_pocaMomKal;//5 hyposis and 3 mom.
+ NTuple::Array<double> m_pocaMomKalP;//5 hyposis and p
+ NTuple::Array<double> m_pocaMomKalPt;//5 hyposis and pt
+ NTuple::Array<int> m_chargeKal;
+ NTuple::Array<double> m_chi2Kal;
+ NTuple::Array<double> m_nDofKal;
+ NTuple::Array<int> m_isFitConverged;
+ NTuple::Array<int> m_isFitConvergedFully;
+ NTuple::Array<int> m_isFitted;
+ NTuple::Item<int> m_nDigi;
+ NTuple::Item<int> m_nHitMc;
+ NTuple::Item<int> m_nSimDCHit;
+ NTuple::Array<int> m_nHitWithFitInfo;
+ NTuple::Item<int> m_nHitKalInput;
+ NTuple::Array<double> m_mdcHitDriftT;
+ NTuple::Array<double> m_mdcHitDriftDl;
+ NTuple::Array<double> m_mdcHitDriftDr;
+ NTuple::Array<int> m_mdcHitLr;
+ NTuple::Array<int> m_mdcHitLayer;
+ NTuple::Array<int> m_mdcHitWire;
+ NTuple::Array<double> m_mdcHitExpDoca;
+ NTuple::Array<double> m_mdcHitExpMcDoca;
+ NTuple::Array<double> m_mdcHitErr;
+ NTuple::Array<int> m_nHitFailedKal;
+ NTuple::Array<int> m_nHitFitted;
+ NTuple::Array<double> m_time;
+ //truth
+ NTuple::Array<int> m_mdcHitMcLr;
+ NTuple::Array<int> m_mdcHitMcTkId;
+ NTuple::Array<double> m_mdcHitMcDrift;
+ NTuple::Array<double> m_mdcHitMcX;
+ NTuple::Array<double> m_mdcHitMcY;
+ NTuple::Array<double> m_mdcHitMcZ;
+ NTuple::Array<double> m_mdcHitExpMcPocaX;
+ NTuple::Array<double> m_mdcHitExpMcPocaY;
+ NTuple::Array<double> m_mdcHitExpMcPocaZ;
+ NTuple::Array<double> m_mdcHitExpMcPocaWireX;
+ NTuple::Array<double> m_mdcHitExpMcPocaWireY;
+ NTuple::Array<double> m_mdcHitExpMcPocaWireZ;
+
+};
+#endif
diff --git a/Reconstruction/Tracking/CMakeLists.txt b/Reconstruction/Tracking/CMakeLists.txt
index 3dd11b3332d73041ae8ca50f21bd48012423c9cd..5123c2922caa0f44e9346f3e11b8bc0710fb04ea 100644
--- a/Reconstruction/Tracking/CMakeLists.txt
+++ b/Reconstruction/Tracking/CMakeLists.txt
@@ -1,4 +1,3 @@
-
# Modules
gaudi_add_module(Tracking
SOURCES src/Clupatra/ClupatraAlg.cpp
diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
index 612fae8e2a0de342361d2ad0527a5c90b9055629..d9bd55f85c3afc6528d45f3b530e5c666b821ec3 100644
--- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
+++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
@@ -1,5 +1,15 @@
#include "TruthTrackerAlg.h"
#include "DataHelper/HelixClass.h"
+#include "DDSegmentation/Segmentation.h"
+#include "DetSegmentation/GridDriftChamber.h"
+#include "DetInterface/IGeomSvc.h"
+#include "DD4hep/Detector.h"
+#include "DD4hep/IDDescriptor.h"
+#include "DD4hep/Plugins.h"
+#include "DD4hep/DD4hepUnits.h"
+
+//external
+#include "CLHEP/Random/RandGauss.h"
#include "edm4hep/EventHeaderCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimTrackerHitCollection.h"
@@ -8,35 +18,30 @@
#include "edm4hep/MCRecoTrackerAssociationCollection.h"
#include "edm4hep/MCRecoParticleAssociationCollection.h"
#include "edm4hep/ReconstructedParticleCollection.h"
-#include "DDSegmentation/Segmentation.h"
-#include "DetSegmentation/GridDriftChamber.h"
-#include "DetInterface/IGeomSvc.h"
-//#include "edm4hep/SimCalorimeterHitCollection.h"
-//#include "edm4hep/CaloHitContributionCollection.h"
-
-#include "DD4hep/Detector.h"
-#include "DD4hep/IDDescriptor.h"
-#include "DD4hep/Plugins.h"
-#include "DD4hep/DD4hepUnits.h"
DECLARE_COMPONENT(TruthTrackerAlg)
TruthTrackerAlg::TruthTrackerAlg(const std::string& name, ISvcLocator* svcLoc)
-: GaudiAlgorithm(name, svcLoc), m_dd4hep(nullptr), m_gridDriftChamber(nullptr),m_decoder(nullptr)
+: GaudiAlgorithm(name, svcLoc),m_dd4hep(nullptr),m_gridDriftChamber(nullptr),
+ m_decoder(nullptr)
{
declareProperty("MCParticle", m_mcParticleCol,
"Handle of the input MCParticle collection");
- declareProperty("DigiDCHitCollection", m_digiDCHitsCol,
- "Handle of digi DCHit collection");
- declareProperty("DCHitAssociationCollection", m_dcHitAssociationCol,
- "Handle of association collection");
- declareProperty("DCTrackCollection", m_dcTrackCol,
- "Handle of association collection");
- declareProperty("DCRecParticleCollection", m_dcRecParticleCol,
+ declareProperty("DigiDCHitCollection", m_DCDigiCol,
+ "Handle of DC digi(TrackerHit) collection");
+ declareProperty("DCHitAssociationCollection", m_DCHitAssociationCol,
+ "Handle of association collection");
+ declareProperty("DCTrackCollection", m_DCTrackCol,
+ "Handle of DC track collection");
+ declareProperty("SiSubsetTrackCollection", m_siSubsetTrackCol,
+ "Handle of silicon subset track collection");
+ declareProperty("SDTTrackCollection", m_SDTTrackCol,
+ "Handle of SDT track collection");
+ declareProperty("DCRecParticleCollection", m_DCRecParticleCol,
"Handle of drift chamber reconstructed particle collection");
- declareProperty("DCRecParticleAssociationCollection", m_dcRecParticleAssociationCol,
+ declareProperty("DCRecParticleAssociationCollection",
+ m_DCRecParticleAssociationCol,
"Handle of drift chamber reconstructed particle collection");
- declareProperty("debug", m_debug=false);
}
StatusCode TruthTrackerAlg::initialize()
@@ -85,16 +90,22 @@ StatusCode TruthTrackerAlg::execute()
}
///Retrieve DC digi
const edm4hep::TrackerHitCollection* digiDCHitsCol=nullptr;
- digiDCHitsCol=m_digiDCHitsCol.get();//FIXME DEBUG
+ digiDCHitsCol=m_DCDigiCol.get();//FIXME DEBUG
if(nullptr==digiDCHitsCol){
debug()<<"TrackerHitCollection not found"<<endmsg;
- //return StatusCode::SUCCESS;
+ //return StatusCode::SUCCESS;//FIXME return when no hits in DC + silicon
}
+ if((int) digiDCHitsCol->size()>m_maxDCDigiCut) return StatusCode::SUCCESS;
///Output Track collection
- edm4hep::TrackCollection* dcTrackCol=m_dcTrackCol.createAndPut();
+ edm4hep::TrackCollection* dcTrackCol=m_DCTrackCol.createAndPut();
+ edm4hep::TrackCollection* sdtTrackCol=m_SDTTrackCol.createAndPut();
+ edm4hep::ReconstructedParticleCollection* dcRecParticleCol(nullptr);
+ if(m_writeRecParticle){
+ dcRecParticleCol=m_DCRecParticleCol.createAndPut();
+ }
////TODO
- /////Output MCRecoTrackerAssociationCollection collection
+ //Output MCRecoTrackerAssociationCollection collection
//const edm4hep::MCRecoTrackerAssociationCollection*
// mcRecoTrackerAssociationCol=nullptr;
//if(nullptr==mcRecoTrackerAssociationCol){
@@ -104,29 +115,91 @@ StatusCode TruthTrackerAlg::execute()
//}
//mcRecoTrackerAssociationCol=m_mcRecoParticleAssociation.get();
- ///Convert MCParticle to Track and ReconstructedParticle
- if(m_debug){
- debug()<<"MCParticleCol size="<<mcParticleCol->size()<<endmsg;
+ ///Retrieve silicon Track
+ const edm4hep::TrackCollection* siTrackCol=nullptr;
+ if(m_siSubsetTrackCol.exist()) siTrackCol=m_siSubsetTrackCol.get();
+ if(nullptr!=siTrackCol) {
+ ///New SDT track
+ for(auto siTrack:*siTrackCol){
+ edm4hep::Track sdtTrack=sdtTrackCol->create();
+ edm4hep::TrackState sdtTrackState;
+ edm4hep::TrackState siTrackStat=siTrack.getTrackStates(0);//FIXME?
+ sdtTrackState.location=siTrackStat.location;
+ sdtTrackState.D0=siTrackStat.D0;
+ sdtTrackState.phi=siTrackStat.phi;
+ sdtTrackState.omega=siTrackStat.omega;
+ sdtTrackState.Z0=siTrackStat.Z0;
+ sdtTrackState.tanLambda=siTrackStat.tanLambda;
+ sdtTrackState.referencePoint=siTrackStat.referencePoint;
+ for(int k=0;k<15;k++){
+ sdtTrackState.covMatrix[k]=siTrackStat.covMatrix[k];
+ }
+ sdtTrack.addToTrackStates(sdtTrackState);
+ sdtTrack.setType(siTrack.getType());
+ sdtTrack.setChi2(siTrack.getChi2());
+ sdtTrack.setNdf(siTrack.getNdf());
+ sdtTrack.setDEdx(siTrack.getDEdx());
+ sdtTrack.setDEdxError(siTrack.getDEdxError());
+ sdtTrack.setRadiusOfInnermostHit(siTrack.getRadiusOfInnermostHit());
+ debug()<<"siTrack trackerHits_size="<<siTrack.trackerHits_size()<<endmsg;
+ for(unsigned int iSiTackerHit=0;iSiTackerHit<siTrack.trackerHits_size();
+ iSiTackerHit++){
+ sdtTrack.addToTrackerHits(siTrack.getTrackerHits(iSiTackerHit));
+ }
+ //TODO tracks
+ for(auto digiDC:*digiDCHitsCol){
+ //if(Sim->MCParti!=current) continue;//TODO
+ sdtTrack.addToTrackerHits(digiDC);
+ }
+ }
}
+
+ ///Convert MCParticle to Track and ReconstructedParticle
+ debug()<<"MCParticleCol size="<<mcParticleCol->size()<<endmsg;
for(auto mcParticle : *mcParticleCol){
- //if(fabs(mcParticle.getCharge()<1e-6) continue;//Skip neutral particles
- edm4hep::Vector3d posV= mcParticle.getVertex();
- edm4hep::Vector3f momV= mcParticle.getMomentum();//GeV
- float pos[3]={posV.x, posV.y, posV.z};
- float mom[3]={momV.x, momV.y, momV.z};
- //FIXME
- //pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom);
+ /// skip mcParticleVertex do not have enough associated hits TODO
+
+ ///Vertex
+ const edm4hep::Vector3d mcParticleVertex=mcParticle.getVertex();//mm
+ edm4hep::Vector3f mcParticleVertexSmeared;//mm
+ mcParticleVertexSmeared.x=
+ CLHEP::RandGauss::shoot(mcParticleVertex.x,m_resVertexX);
+ mcParticleVertexSmeared.y=
+ CLHEP::RandGauss::shoot(mcParticleVertex.y,m_resVertexY);
+ mcParticleVertexSmeared.z=
+ CLHEP::RandGauss::shoot(mcParticleVertex.z,m_resVertexZ);
+ ///Momentum
+ edm4hep::Vector3f mcParticleMom=mcParticle.getMomentum();//GeV
+ double mcParticlePt=sqrt(mcParticleMom.x*mcParticleMom.x+
+ mcParticleMom.y*mcParticleMom.y);
+ //double mcParticlePtSmeared=
+ // CLHEP::RandGauss::shoot(mcParticlePt,m_resPT);
+ double mcParticleMomPhi=atan2(mcParticleMom.y,mcParticleMom.x);
+ double mcParticleMomPhiSmeared=
+ CLHEP::RandGauss::shoot(mcParticleMomPhi,m_resMomPhi);
+ edm4hep::Vector3f mcParticleMomSmeared;
+ mcParticleMomSmeared.x=mcParticlePt*cos(mcParticleMomPhiSmeared);
+ mcParticleMomSmeared.y=mcParticlePt*sin(mcParticleMomPhiSmeared);
+ mcParticleMomSmeared.z=
+ CLHEP::RandGauss::shoot(mcParticleMom.z,m_resPz);
+
+ ///Converted to Helix
double B[3]={1e9,1e9,1e9};
m_dd4hepField.magneticField({0.,0.,0.},B);
- if(m_debug)std::cout<<" PDG "<<mcParticle.getPDG()<<" charge "
- << mcParticle.getCharge()
- <<" pos "<<pos[0]<<" "<<pos[1]<<" "<<pos[2]
- <<" mom "<<mom[0]<<" "<<mom[1]<<" "<<mom[2]
- <<" Bxyz "<<B[0]/dd4hep::tesla<<" "<<B[1]/dd4hep::tesla
- <<" "<<B[2]/dd4hep::tesla<<" tesla"<<std::endl;
-
HelixClass helix;
- helix.Initialize_VP(pos,mom,(int) mcParticle.getCharge(),B[2]/dd4hep::tesla);
+ //float pos[3]={mcParticleVertexSmeared.x,
+ // mcParticleVertexSmeared.y,mcParticleVertexSmeared.z};
+ //float mom[3]={mcParticleMomSmeared.x,mcParticleMomSmeared.y,
+ // mcParticleMomSmeared.z};
+ ////FIXME DEBUG
+ float pos[3]={(float)mcParticleVertex.x,
+ (float)mcParticleVertex.y,(float)mcParticleVertex.z};
+ float mom[3]={(float)mcParticleMom.x,(float)mcParticleMom.y,
+ (float)mcParticleMom.z};
+ helix.Initialize_VP(pos,mom,mcParticle.getCharge(),B[2]/dd4hep::tesla);
+
+ ///new Track
+ edm4hep::Track dcTrack=dcTrackCol->create();
edm4hep::TrackState trackState;
trackState.D0=helix.getD0();
trackState.phi=helix.getPhi0();
@@ -137,40 +210,61 @@ StatusCode TruthTrackerAlg::execute()
std::array<float,15> covMatrix;
for(int i=0;i<15;i++){covMatrix[i]=999.;}//FIXME
trackState.covMatrix=covMatrix;
-
- if(m_debug){
- std::cout<<" helix d0 "<<trackState.D0<<
- " phi "<<trackState.phi<<
- " omega "<<trackState.omega<<
- " z0 "<<trackState.Z0<<
- " tanLambda "<<trackState.tanLambda<<
- " referencePoint "<<trackState.referencePoint<< std::endl;
- }
- //new Track
- edm4hep::Track track = dcTrackCol->create();
- track.addToTrackStates(trackState);
- if(nullptr!=digiDCHitsCol) {
- //track.setType();//TODO
- //track.setChi2(trackState->getChi2());
- track.setNdf(digiDCHitsCol->size()-5);
- //track.setDEdx();//TODO
- //track.setRadiusOfInnermostHit();//TODO
- for(int i=0; i<digiDCHitsCol->size(); i++ ){
- edm4hep::TrackerHit digiDC=digiDCHitsCol->at(i);
- //if(Sim->MCParti!=current) continue;//TODO
- track.addToTrackerHits(digiDC);
- }
+ dcTrack.addToTrackStates(trackState);
+ //dcTrack.setType();//TODO
+ //dcTrack.setChi2(gauss(digiDCHitsCol->size-5(),1));//FIXME
+ dcTrack.setNdf(digiDCHitsCol->size()-5);
+ //dcTrack.setDEdx();//TODO
+ //set hits
+ double radiusOfInnermostHit=1e9;
+ debug()<<"digiDCHitsCol size"<<digiDCHitsCol->size()<<endmsg;
+ for(auto digiDC : *digiDCHitsCol){
+ //if(Sim->MCParti!=current) continue;//TODO
+ edm4hep::Vector3d digiPos=digiDC.getPosition();
+ double r=sqrt(digiPos.x*digiPos.x+digiPos.y*digiPos.y);
+ if(r<radiusOfInnermostHit) radiusOfInnermostHit=r;
+ dcTrack.addToTrackerHits(digiDC);
}
- //TODO
- //new ReconstructedParticle
- //recParticle->setType();
- //dcRecParticle->setEnergy();
- //edm4hep::Vector3f momVec3(helix.getMomentum()[0],
- // helix.getMomentum()[1],helix.getMomentum()[2]);
- //recParticle->setMomentum(momVec3);
- }
+ dcTrack.setRadiusOfInnermostHit(radiusOfInnermostHit);//TODO
+
+ edm4hep::ReconstructedParticle dcRecParticle;
+ if(m_writeRecParticle){
+ dcRecParticle=dcRecParticleCol->create();
+ ///new ReconstructedParticle
+ //dcRecParticle.setType();//TODO
+ double mass=mcParticle.getMass();
+ double p=sqrt(mcParticleMomSmeared.x*mcParticleMomSmeared.x+
+ mcParticleMomSmeared.y*mcParticleMomSmeared.y+
+ mcParticleMomSmeared.z*mcParticleMomSmeared.z);
+ dcRecParticle.setEnergy(sqrt(mass*mass+p*p));
+ dcRecParticle.setMomentum(mcParticleMomSmeared);
+ dcRecParticle.setReferencePoint(mcParticleVertexSmeared);
+ dcRecParticle.setCharge(mcParticle.getCharge());
+ dcRecParticle.setMass(mass);
+ //dcRecParticle.setGoodnessOfPID();//TODO
+ //std::array<float>,10> covMat=?;//TODO
+ //dcRecParticle.setCovMatrix(covMat);//TODO
+ //dcRecParticle.setStartVertex();//TODO
+ edm4hep::ParticleID particleID(0,mcParticle.getPDG(),0,1);//FIXME
+ dcRecParticle.setParticleIDUsed(particleID);
+ dcRecParticle.addToTracks(dcTrack);
+ }//end of write RecParticle
+
+ debug()<<"mcParticle "<<mcParticle
+ <<" momPhi "<<mcParticleMomPhi
+ <<" mcParticleVertex("<<mcParticleVertex<<")mm "
+ <<" mcParticleVertexSmeared("<<mcParticleVertexSmeared<<")mm "
+ <<" mcParticleMom("<<mcParticleMom<<")GeV "
+ <<" mcParticleMomSmeared("<<mcParticleMomSmeared<<")GeV "
+ <<" Bxyz "<<B[0]/dd4hep::tesla<<" "<<B[1]/dd4hep::tesla
+ <<" "<<B[2]/dd4hep::tesla<<" tesla"<<endmsg;
+ debug()<<"trackState:location,D0,phi,omega,Z0,tanLambda"
+ <<",referencePoint,cov"<<std::endl<<trackState<<std::endl;
+ debug()<<"dcTrack"<<dcTrack<<endmsg;
+ if(m_writeRecParticle) debug()<<"dcRecParticle"<<dcRecParticle<<endmsg;
+ }//end loop over MCParticleCol
- if(m_debug) std::cout<<"Output DCTrack size="<<dcTrackCol->size()<<std::endl;
+ debug()<<"Output DCTrack size="<<dcTrackCol->size()<<endmsg;
return StatusCode::SUCCESS;
}
diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
index 23e581e1b73a4e1f962ddc3e53cfb8e525258de6..30efd0a8a71f7d4ec4211fcb6e074e781ecc9fa0 100644
--- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
+++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
@@ -27,9 +27,9 @@ class TruthTrackerAlg: public GaudiAlgorithm
public:
TruthTrackerAlg(const std::string& name, ISvcLocator* svcLoc);
- virtual StatusCode initialize();
- virtual StatusCode execute();
- virtual StatusCode finalize();
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute() override;
+ virtual StatusCode finalize() override;
private:
SmartIF<IGeomSvc> m_geomSvc;
@@ -41,25 +41,38 @@ class TruthTrackerAlg: public GaudiAlgorithm
//reader
DataHandle<edm4hep::MCParticleCollection> m_mcParticleCol{
"MCParticle", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_digiDCHitsCol{
- "DigiDCHitsCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_DCDigiCol{
+ "DigiDCHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::MCRecoTrackerAssociationCollection>
- m_dcHitAssociationCol{ "DCHitAssociationCollection",
+ m_DCHitAssociationCol{ "DCHitAssociationCollection",
+ Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackCollection>
+ m_siSubsetTrackCol{ "SiSubsetTrackCollection",
Gaudi::DataHandle::Reader, this};
//writer
- DataHandle<edm4hep::TrackCollection> m_dcTrackCol{"DCTrackCollection",
- Gaudi::DataHandle::Writer, this};
- DataHandle<edm4hep::ReconstructedParticleCollection> m_dcRecParticleCol{
+ DataHandle<edm4hep::TrackCollection> m_DCTrackCol{
+ "DCTrackCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::TrackCollection> m_SDTTrackCol{
+ "SDTTrackCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::ReconstructedParticleCollection> m_DCRecParticleCol{
"DCRecParticleCollection", Gaudi::DataHandle::Writer, this};
DataHandle<edm4hep::MCRecoParticleAssociationCollection>
- m_dcRecParticleAssociationCol{"DCRecMCRecoParticleAssociationCollection",
+ m_DCRecParticleAssociationCol{
+ "DCRecMCRecoParticleAssociationCollection",
Gaudi::DataHandle::Writer, this};
//readout for getting segmentation
Gaudi::Property<std::string> m_readout_name{this, "readout",
"DriftChamberHitsCollection"};
-
- Gaudi::Property<int> m_debug{ this, "debug", false};
+ Gaudi::Property<bool> m_writeRecParticle{this,"writeRecParticle",false};
+ Gaudi::Property<float> m_resPT{this,"resPT",0};//ratio
+ Gaudi::Property<float> m_resPz{this,"resPz",0};//ratio
+ Gaudi::Property<float> m_resMomPhi{this,"resMomPhi",0};//radian
+ Gaudi::Property<float> m_resMomTheta{this,"resMomTheta",0};//radian
+ Gaudi::Property<float> m_resVertexX{this,"resVertexX",0.003};//3um
+ Gaudi::Property<float> m_resVertexY{this,"resVertexY",0.003};//3um
+ Gaudi::Property<float> m_resVertexZ{this,"resVertexZ",0.003};//3um
+ Gaudi::Property<int> m_maxDCDigiCut{this,"maxDigiCut",1e6};
};
#endif