diff --git a/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01_noECAL.xml b/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01_noECAL.xml
new file mode 100644
index 0000000000000000000000000000000000000000..256dcf001804c7c6d46b0ff47d5f2d403615115c
--- /dev/null
+++ b/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01_noECAL.xml
@@ -0,0 +1,51 @@
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
+ xmlns:xs="http://www.w3.org/2001/XMLSchema"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
+ <info name="CRD_i1_v01"
+ title="CepC reference detctor with coil inside Hcal"
+ author="C.D.Fu, "
+ url="http://cepc.ihep.ac.cn"
+ status="developing"
+ version="v01">
+ <comment>CepC reference detector simulation models used for detector study </comment>
+ </info>
+
+ <includes>
+ <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/>
+ <gdmlFile ref="../CRD_common_v01/materials.xml"/>
+ </includes>
+
+ <define>
+ <constant name="world_size" value="2.6*m"/>
+ <constant name="world_x" value="world_size"/>
+ <constant name="world_y" value="world_size"/>
+ <constant name="world_z" value="world_size"/>
+
+ <include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml"/>
+ </define>
+
+ <include ref="./CRD_Dimensions_v01_01.xml"/>
+
+ <include ref="../CRD_common_v01/Beampipe_v01_01.xml"/>
+ <include ref="../CRD_common_v01/VXD_v01_01.xml"/>
+ <include ref="../CRD_common_v01/FTD_SimpleStaggered_v01_01.xml"/>
+ <include ref="../CRD_common_v01/SIT_SimplePlanar_v01_01.xml"/>
+ <include ref="../CRD_common_v01/DC_Simple_v01_01.xml"/>
+ <include ref="../CRD_common_v01/SET_SimplePlanar_v01_01.xml"/>
+ <!-- <include ref="../CRD_common_v01/Ecal_Crystal_Barrel_v01_01.xml"/> -->
+ <!--include ref="../CRD_common_v01/Ecal_Crystal_Endcap_v01_01.xml"/-->
+ <!--include ref="../CRD_common_v01/Coil_v01_01.xml"/-->
+ <!--include ref="../CRD_common_v01/Hcal_v01_01.xml"/-->
+ <!--include ref="../CRD_common_v01/Yoke_v01_01.xml"/-->
+ <!--include ref="../CRD_common_v01/Lcal_v01_01.xml"/-->
+
+ <fields>
+ <field name="GlobalSolenoid" type="solenoid"
+ inner_field="Field_nominal_value"
+ outer_field="Field_outer_nominal_value"
+ zmax="SolenoidCoil_half_length"
+ outer_radius="SolenoidCoil_center_radius">
+ </field>
+ </fields>
+
+</lccdd>
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_truthTracker_fit_DC.py b/Examples/options/tut_detsim_digi_truthTracker_fit_DC.py
new file mode 100644
index 0000000000000000000000000000000000000000..52a7c026596dcee1047590b74bf08c513b1966a3
--- /dev/null
+++ b/Examples/options/tut_detsim_digi_truthTracker_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
index 920ae7b968f2573c5ec3b9d20630c77409765b01..758bdef08a6a8468731996e1c3696c30f913f5fe 100644
--- a/Reconstruction/RecGenfitAlg/CMakeLists.txt
+++ b/Reconstruction/RecGenfitAlg/CMakeLists.txt
@@ -1,39 +1,33 @@
-gaudi_subdir(RecGenfitAlg v0r0)
-
-find_package(CLHEP REQUIRED;CONFIG)
-find_package(GSL REQUIRED )
-find_package(LCIO REQUIRED )
-find_package(podio REQUIRED )
-find_package(EDM4HEP REQUIRED)
-find_package(ROOT REQUIRED Geom)
-find_package(DD4hep COMPONENTS DDCore DDRec DDParsers REQUIRED)
-
-gaudi_depends_on_subdirs(
- Detector/DetInterface
- Detector/DetSegmentation
- Utilities/DataHelper
+# RecGenfitAlg
+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
)
-set (GenFit_INCLUDE_DIRS $ENV{GENFIT_ROOT}/include)
-set (GenFit_LIB_DIRS $ENV{GENFIT_ROOT}/lib64)
-set (Eigen_INCLUDE_DIRS $ENV{Eigen_ROOT}/include/eigen3)
-
-set(RecGenfitAlg_srcs
- src/RecGenfitAlgDC.cpp
- src/GenfitTrack.cpp
- src/GenfitField.cpp
- src/GenfitFitter.cpp
- src/GenfitMaterialInterface.cpp
- src/GenfitMsg.cpp
- )
+target_include_directories(RecGenfitAlg PUBLIC
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/include
+ $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
-# Modules
-gaudi_add_module(RecGenfitAlg ${RecGenfitAlg_srcs}
- INCLUDE_DIRS k4FWCore GaudiKernel GaudiAlgLib CLHEP ROOT gear
- ${GSL_INCLUDE_DIRS} ${LCIO_INCLUDE_DIRS} ${GenFit_INCLUDE_DIRS}
- ${Eigen_INCLUDE_DIRS}
- LINK_LIBRARIES k4FWCore GaudiKernel GaudiAlgLib CLHEP ROOT DataHelperLib
- DetSegmentation $ENV{GEAR}/lib/libgearsurf.so ${GSL_LIBRARIES} ${LCIO_LIBRARIES}
- ${GenFit_LIB_DIRS}/libgenfit2.so DD4hep ${DD4hep_COMPONENT_LIBRARIES}
- -Wl,--no-as-needed EDM4HEP::edm4hep EDM4HEP::edm4hepDict -Wl,--as-needed
- )
+install(TARGETS RecGenfitAlg
+ EXPORT CEPCSWTargets
+ RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+ LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+ COMPONENT dev)
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
index de6df54222b37684face8a615ea82dad76879f84..341e52f3d9c1a0e6539b62bec6225410dbbad7d4 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
+++ b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
@@ -43,7 +43,6 @@ GenfitField::get(const double& posX, const double& posY, const double& posZ,
dd4hep::Position(posX/dd4hep::cm,posY/dd4hep::cm,posZ/dd4hep::cm));
//m_dd4hepField.magneticField(pos,B);
- double factor = 1./dd4hep::kilogauss;
Bx=field.X()/dd4hep::kilogauss;
By=field.Y()/dd4hep::kilogauss;
Bz=field.Z()/dd4hep::kilogauss;
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
index 2a7745b60d3bea05429188bc5c6d7a9e8992ee78..15d40830f9d3fa8deb71e0ca3e306ee7e5ad7c37 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
+++ b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
@@ -116,9 +116,7 @@ int GenfitFitter::init(bool deleteOldFitter)
return -1;
}
GenfitMsg::get()<<MSG::DEBUG<<"Fitter type is "<<m_fitterType<<endmsg;
-std::cout<<__FILE__<<" "<<__LINE__<<m_absKalman<<std::endl;
m_absKalman->setDebugLvl(m_debug);
-std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
return 0;
}
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
index a73f6108e66874090ca1313be8ddd748282fc9e2..f210dc9f69c88a0cbc6407b18911c181c700accf 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
@@ -4,13 +4,17 @@
//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/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
#include "edm4hep/Vector3d.h"
#include "DetSegmentation/GridDriftChamber.h"
@@ -33,17 +37,24 @@
#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::SimTrackerHit hit1,edm4hep::SimTrackerHit hit2)
+sortDCHit(edm4hep::ConstSimTrackerHit hit1,edm4hep::ConstSimTrackerHit hit2)
{
- return (hit1.getTime()<hit2.getTime());
+ //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_genfitField(genfitField),m_gridDriftChamber(seg),m_name(name), m_track(nullptr) ,m_reps() ,m_debug(0)
+:m_name(name),m_track(nullptr),m_reps(),m_debug(0),
+m_genfitField(genfitField),m_gridDriftChamber(seg)
{
}
@@ -92,12 +103,12 @@ bool GenfitTrack::createGenfitTrack(int pdgType,int charge,
/// new a track representation and add to the track
int chargeId=0;
- charge<0 ? chargeId=0 : chargeId=1;
+ 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 pdg "
- <<s_PDG[chargeId][pdgType]<<endmsg;
+ <<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]);
@@ -154,9 +165,9 @@ bool GenfitTrack::createGenfitTrackFromMCParticle(int pidType,
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;
+ //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
@@ -185,34 +196,74 @@ bool GenfitTrack::createGenfitTrackFromEDM4HepTrack(int pidType,
return true;
}
-/// Add a 3d SpacepointMeasurement with MC truth position smeared by sigma
-bool GenfitTrack::addSpacePointMeasurementOnTrack(const TVectorD& pos,
- bool smear, double resolution, int detID, int hitID)
+/// 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;
+ 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(pos_t);
- if (smear) {
- for (int ii=0;ii<3;ii++){
- pos_smeared[ii] += gRandom->Gaus(0, resolution/TMath::Sqrt(3.));
- }
+ 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) = resolution*resolution;
- hitCov(1,1) = resolution*resolution;
- hitCov(2,2) = resolution*resolution;
+ 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<<" addSpacePointMeasurementOnTrack "
- <<hitID<<" " <<pos_t[0]<<" "<<pos_t[1]<<" "<<pos_t[2]<<"cm smeared "
- <<pos_smeared[0]<<" "<<pos_smeared[1]<<" "<<pos_smeared[2]
- <<" res "<<resolution<<" cm"<<endmsg;
+ 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);
@@ -249,7 +300,7 @@ void GenfitTrack::addWireMeasurement(double driftDistance,
m_track->insertPoint(trackPoint);
- }catch(genfit::Exception e){
+ }catch(genfit::Exception& e){
GenfitMsg::get()<<MSG::DEBUG<<m_name
<<"Add wire measurement exception"<<endmsg;
e.what();
@@ -259,7 +310,7 @@ void GenfitTrack::addWireMeasurement(double driftDistance,
//Add wire measurement on wire, unit conversion here
bool GenfitTrack::addWireMeasurementOnTrack(edm4hep::Track& track,double sigma)
{
- for(int iHit=0;iHit<track.trackerHits_size();iHit++){
+ for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){
edm4hep::ConstTrackerHit hit=track.getTrackerHits(iHit);
double driftVelocity=40;//FIXME, TODO, um/ns
@@ -269,10 +320,10 @@ bool GenfitTrack::addWireMeasurementOnTrack(edm4hep::Track& track,double sigma)
m_gridDriftChamber->cellposition(hit.getCellID(),endPointStart,
endPointEnd);
int lrAmbig=0;
- std::cout<<__FILE__<<" time "<<hit.getTime()<<" driftVelocity "
- <<driftVelocity<<std::endl;
- std::cout<<__FILE__<<" wire pos " <<endPointStart.X()<<" "
- <<endPointStart.Y()<<" " <<endPointStart.Z()<<" "
+ 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);
@@ -294,7 +345,7 @@ bool GenfitTrack::getMOP(int hitID,
if(nullptr == trackRep) trackRep = getRep();
try{
mop = m_track->getFittedState(hitID,trackRep);
- }catch(genfit::Exception e){
+ }catch(genfit::Exception& e){
e.what();
return false;
}
@@ -476,7 +527,7 @@ int GenfitTrack::getFittedState(TLorentzVector& pos, TVector3& mom,
genfit::MeasuredStateOnPlane mop;
try{
mop = m_track->getFittedState(biased);
- }catch(genfit::Exception e){
+ }catch(genfit::Exception& e){
std::cout<<" getNumPointsWithFittedInfo "
<<getNumPointsWithFittedInfo(repID)
<<" no TrackPoint with fitted info "<<std::endl;
@@ -573,7 +624,7 @@ double GenfitTrack::extrapolateToHit( TVector3& poca, TVector3& pocaDir,
// pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise);
extrapoLen = rep->extrapolateToLine(state, end0, end1, poca,
pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise);
- } catch (genfit::Exception& e) {
+ }catch(genfit::Exception& e) {
GenfitMsg::get() << MSG::ERROR
<<"Exception in GenfigFitter::ExtrapolateToHit"<<e.what()<<endmsg;
return extrapoLen;
@@ -593,33 +644,85 @@ double GenfitTrack::extrapolateToHit( TVector3& poca, TVector3& pocaDir,
}//end of ExtrapolateToHit
-///Add space point measurement to track
-int GenfitTrack::addSpacePointMeasurementOnTrack(
- const edm4hep::MCParticle& mcParticle,
- const edm4hep::SimTrackerHitCollection*& simDCHitCol,
- bool smear, double resolution) {
- //Sort sim DC hit by time
- std::vector<edm4hep::SimTrackerHit> sortedSimDCHitCol;
- for(auto simDCHit: *simDCHitCol) sortedSimDCHitCol.push_back(simDCHit);
- std::sort(sortedSimDCHitCol.begin(),sortedSimDCHitCol.end(),sortDCHit);
+///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(auto simDCHit: sortedSimDCHitCol){
- int rand=gRandom->Integer(sortedSimDCHitCol.size());
- if(rand>200){ continue; }
- edm4hep::Vector3d pos=simDCHit.getPosition();
+ 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*dd4hep::mm;
- p[1]=pos.y*dd4hep::mm;
- p[2]=pos.z*dd4hep::mm;
- unsigned long long detID = simDCHit.getCellID();
- if(addSpacePointMeasurementOnTrack(p,smear,resolution,detID,hitID)){
+ 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<<"addSpacePointMeasurementOnTrack "
+ GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurement"
<<detID<<" faieled" <<endmsg;
}
}
+
GenfitMsg::get()<<MSG::DEBUG<<"GenfitTrack nHitAdded "<<hitID<<endmsg;
return hitID;
}
@@ -630,8 +733,10 @@ bool GenfitTrack::storeTrack(edm4hep::ReconstructedParticle& recParticle,
/// Get fit status
const genfit::FitStatus* fitState = m_track->getFitStatus();
- double chi2 = fitState->getChi2();
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();
@@ -659,8 +764,8 @@ bool GenfitTrack::storeTrack(edm4hep::ReconstructedParticle& recParticle,
" ndf "<<ndf
<<endmsg;
}
- float pos[3]={fittedPos.X(),fittedPos.Y(),fittedPos.Z()};
- float mom[3]={fittedMom.X(),fittedMom.Y(),fittedMom.Z()};
+ 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()));
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
index 6e2823e46e57f42b6c19696388381f184b861439..708c33e38ba282c92f5d0f8ef579fde688d150f4 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
@@ -41,7 +41,9 @@ namespace edm4hep{
class MCParticle;
class SimTrackerHitCollection;
class ReconstructedParticle;
+ class MCRecoTrackerAssociationCollection;
class Track;
+ class ConstTrackerHit;
class Vector3d;
class Vector3f;
}
@@ -96,8 +98,11 @@ class GenfitTrack {
// int PrepareHits();//TODO
/// Add a space point measurement, return number of hits on track
- virtual bool addSpacePointMeasurementOnTrack(const TVectorD&, bool, double,
- int detID=-1, int hitID=-1);
+ 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,
@@ -108,13 +113,13 @@ class GenfitTrack {
virtual bool addWireMeasurementOnTrack(edm4hep::Track& track, double sigma);
///Add space point from truth to track
- int addSpacePointMeasurementOnTrack(const edm4hep::MCParticle& mcParticle,
- const edm4hep::SimTrackerHitCollection*& simDCHitCol,bool smear=true,
- double resolution=0.01);
+ 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, double chi2Cut);
+ 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,
diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
index 99d7d3a5f2c273d42c805c877d662f6a0797ce21..9de410346dc9fe0ce43a4fcc1886e24e8f051c47 100644
--- a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
+++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
@@ -36,35 +36,11 @@
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)
+/////////////////////////////////////////////////////////////////////
+RecGenfitAlgDC::RecGenfitAlgDC(const std::string& name, ISvcLocator* pSvcLocator):
+ GaudiAlgorithm(name, pSvcLocator),m_nPDG(5),m_dd4hep(nullptr),
+ m_gridDriftChamber(nullptr),m_decoder(nullptr)
{
- declareProperty("inputType", m_inputType=0);//0:MC particle and hits
- declareProperty("debug", m_debug=0);
- //DAF,DAFRef,KalmanFitter,KalmanFitterRefTrack
- declareProperty("fitterType",m_fitterType="DAFRef");
- declareProperty("correctBremsstrahlung", m_correctBremsstrahlung=false);
- declareProperty("noMaterialEffects", m_noMaterialEffects=false);
- declareProperty("maxIteration", m_maxIteration=20);
- declareProperty("resortHits", m_resortHits=true);
- declareProperty("ndfCut", m_ndfCut=1);
- declareProperty("chi2Cut", m_chi2Cut=1000);
- declareProperty("bStart", m_bStart=100);
- declareProperty("bFinal", m_bFinal=0.01);
- //mdc corner cuts, negtive for no cut
- declareProperty("dcCornerCuts", m_dcCornerCuts=-999);
- //-1,chargedGeantino;0,1,2,3,4:e,mu,pi,K,proton
- declareProperty("debugPid", m_debugPid=-99);
- declareProperty("useTruthSeed", m_useTruthSeed=false);
- declareProperty("useRecLRAmbig", m_useRecLRAmbig=true);
- declareProperty("genfitHistRootName", m_genfitHistRootName="");
- declareProperty("useTruthHit", m_useTruthHit=false);
- declareProperty("showDisplay", m_showDisplay=false);
- declareProperty("initCovResPos", m_initCovResPos=1);
- declareProperty("initCovResMom", m_initCovResMom=0.1);
-
//declareProperty("EventHeaderCol", _headerCol);
declareProperty("MCParticle", m_mcParticleCol,
"Handle of the input MCParticle collection");
@@ -74,8 +50,8 @@ RecGenfitAlgDC::RecGenfitAlgDC(const std::string& name, ISvcLocator* pSvcLocator
"Handle of digi DCHit collection");
declareProperty("DCTrackCollection", m_dcTrackCol,
"Handle of DC track collection");
- //declareProperty("DCHitAssociationCollection", _dcHitAssociationCol,
- // "Handle of association collection");
+ declareProperty("DCHitAssociationCollection", m_dcHitAssociationCol,
+ "Handle of simTrackerHit and TrackerHit association collection");
declareProperty("DCRecParticleCollection", m_dcRecParticleCol,
"Handle of drift chamber reconstructed particle collection");
}
@@ -85,7 +61,7 @@ RecGenfitAlgDC::RecGenfitAlgDC(const std::string& name, ISvcLocator* pSvcLocator
StatusCode RecGenfitAlgDC::initialize()
{
MsgStream log(msgSvc(), name());
- info()<<" RecGenfitAlgDC initialize()"<<endmsg;
+ info()<<"RecGenfitAlgDC initialize()"<<endmsg;
///Get GeomSvc
m_geomSvc=Gaudi::svcLocator()->service("GeomSvc");
@@ -99,7 +75,7 @@ StatusCode RecGenfitAlgDC::initialize()
m_dd4hepField=m_geomSvc->lcdd()->field();
/// New a genfit fitter
- m_genfitFitter=new GenfitFitter(m_fitterType.c_str());
+ 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());
@@ -123,31 +99,27 @@ StatusCode RecGenfitAlgDC::initialize()
m_nDCTrack=0;
///Get Readout
dd4hep::Readout readout=m_dd4hep->readout(m_readout_name);
-std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
///Get Segmentation
m_gridDriftChamber=dynamic_cast<dd4hep::DDSegmentation::GridDriftChamber*>
(readout.segmentation().segmentation());
-std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
if(nullptr==m_gridDriftChamber){
error() << "Failed to get the GridDriftChamber" << endmsg;
return StatusCode::FAILURE;
}
-std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
///Get Decoder
m_decoder = m_geomSvc->getDecoder(m_readout_name);
if (nullptr==m_decoder) {
error() << "Failed to get the decoder" << endmsg;
return StatusCode::FAILURE;
}
-std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
///book tuple
- NTuplePtr nt(ntupleSvc(), "RecGenfitAlgDC/genfit");
+ NTuplePtr nt(ntupleSvc(), "RecGenfitAlgDC/recGenfitAlgDC");
if(nt){
m_tuple=nt;
}else{
- m_tuple=ntupleSvc()->book("RecGenfitAlgDC/genfit",
+ m_tuple=ntupleSvc()->book("RecGenfitAlgDC/recGenfitAlgDC",
CLID_ColumnWiseTuple,"RecGenfitAlgDC");
if(m_tuple){
StatusCode sc;
@@ -155,7 +127,7 @@ std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
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_pocaPosMc,3);
+ 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);
@@ -170,49 +142,51 @@ std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
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("nDofKal",m_pidIndex,m_nDofKal);
- sc=m_tuple->addItem("chi2Kal",m_pidIndex,m_chi2Kal);
- sc=m_tuple->addItem("convergeKal",m_pidIndex,m_convergeKal);
- sc=m_tuple->addItem("isFittedKal",m_pidIndex,m_isFittedKal);
- sc=m_tuple->addItem("nHitFailedKal",m_pidIndex,m_nHitFailedKal);
- sc=m_tuple->addItem("nHitFitted",m_pidIndex,m_nHitFitted);
+ 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("nClusterCgem",m_nClusterCgem);
- sc=m_tuple->addItem("nHitMdc",m_nHitMdc,0,6796);
- sc=m_tuple->addItem("mdcHitDriftT",m_nHitMdc,m_mdcHitDriftT);
- sc=m_tuple->addItem("mdcHitDriftDl",m_nHitMdc,m_mdcHitDriftDl);
- sc=m_tuple->addItem("mdcHitDriftDr",m_nHitMdc,m_mdcHitDriftDr);
- sc=m_tuple->addItem("mdcHitLr",m_nHitMdc,m_mdcHitLr);
- sc=m_tuple->addItem("mdcHitLayer",m_nHitMdc,m_mdcHitLayer);
- sc=m_tuple->addItem("mdcHitWire",m_nHitMdc,m_mdcHitWire);
- sc=m_tuple->addItem("mdcHitExpDoca",m_nHitMdc,m_mdcHitExpDoca);
- sc=m_tuple->addItem("mdcHitExpMcDoca",m_nHitMdc,m_mdcHitExpMcDoca);
- sc=m_tuple->addItem("mdcHitErr",m_nHitMdc,m_mdcHitErr);
- sc=m_tuple->addItem("time",m_pidIndex,m_time);
- sc=m_tuple->addItem("mdcHitMcTkId",m_nHitMdc,m_mdcHitMcTkId);
- sc=m_tuple->addItem("mdcHitMcLr",m_nHitMdc,m_mdcHitMcLr);
- sc=m_tuple->addItem("mdcHitMcDrift",m_nHitMdc,m_mdcHitMcDrift);
- sc=m_tuple->addItem("mdcHitMcX",m_nHitMdc,m_mdcHitMcX);
- sc=m_tuple->addItem("mdcHitMcY",m_nHitMdc,m_mdcHitMcY);
- sc=m_tuple->addItem("mdcHitMcZ",m_nHitMdc,m_mdcHitMcZ);
- sc=m_tuple->addItem("mcPocaX",m_nHitMdc,m_mdcHitExpMcPocaX);
- sc=m_tuple->addItem("mcPocaY",m_nHitMdc,m_mdcHitExpMcPocaY);
- sc=m_tuple->addItem("mcPocaZ",m_nHitMdc,m_mdcHitExpMcPocaZ);
- sc=m_tuple->addItem("mcPocaWireX",m_nHitMdc,m_mdcHitExpMcPocaWireX);
- sc=m_tuple->addItem("mcPocaWireY",m_nHitMdc,m_mdcHitExpMcPocaWireY);
- sc=m_tuple->addItem("mcPocaWireZ",m_nHitMdc,m_mdcHitExpMcPocaWireZ);
- log << MSG::DEBUG<< "Book tuple RecGenfitAlgDC/genfit" << endmsg;
+ 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{
- log << MSG::ERROR << "Cannot book tuple RecGenfitAlgDC/genfit" << endmsg;
+ 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;
}
@@ -220,11 +194,9 @@ std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
StatusCode RecGenfitAlgDC::execute()
{
-
- m_timer=clock();
- m_firstPid=true;
- info()<<" RecGenfitAlgDC in execute()"<<endmsg;
StatusCode sc;
+ m_timer=clock();
+ info()<<"RecGenfitAlgDC in execute()"<<endmsg;
/////retrieve EventHeader
//auto header = _headerCol.get()->at(0);
@@ -234,58 +206,33 @@ StatusCode RecGenfitAlgDC::execute()
//std::cout<<"run "<<header.getEventNumber()
// <<" "<<header.getRunNumber()<<std::endl;
- fitMC();
-
- //if(m_genfitDisplay) while(1){
- // std::cout<<"Press any key to finish..."<<std::endl;
- // //system ("pause");
- //}
-
- if(m_tuple) m_tuple->write();
-
- return StatusCode::SUCCESS;
-}
-
-StatusCode RecGenfitAlgDC::fitMC()
-{
- MsgStream log(msgSvc(), name());
-
- /////retrieve DC hits and MC particle(s)
- //const edm4hep::SimTrackerHitCollection* simDCHitCol=nullptr;
- //simDCHitCol=m_simDCHitCol.get();
- //if(nullptr==simDCHitCol) {
- // log<<MSG::DEBUG<<"DriftChamberHitsCollection not found"<<endmsg;
- // return StatusCode::SUCCESS;
- //}
- //const edm4hep::MCParticleCollection* mcParticleCol = nullptr;
- //mcParticleCol=m_mcParticleCol.get();
- //if(!mcParticleCol) {
- // log<<MSG::DEBUG<<"MCParticleCollection not found"<<endmsg;
- // return StatusCode::SUCCESS;
- //}
///retrieve Track and TrackHits
const edm4hep::TrackCollection* dcTrackCol=nullptr;
dcTrackCol=m_dcTrackCol.get();
if(nullptr==dcTrackCol) {
- log<<MSG::DEBUG<<"TrackCollection not found"<<endmsg;
+ debug()<<"TrackCollection not found"<<endmsg;
return StatusCode::SUCCESS;
}
const edm4hep::TrackerHitCollection* didiDCHitsCol=nullptr;
didiDCHitsCol=m_digiDCHitsCol.get();
if(nullptr==didiDCHitsCol) {
- log<<MSG::DEBUG<<"DigiDCHitsCollection not found"<<endmsg;
+ 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
///----------------------------------------------------
- log<<MSG::DEBUG<<"DCTrackCol size="<<dcTrackCol->size()<<endmsg;
+ 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-1;pidType++){
+ ///-1 for chargedgeantino
+ for(unsigned int pidType=0;pidType<m_nPDG;pidType++){
if((m_debugPid>=0) && (m_debugPid!=pidType)) continue;
///-----------------------------------
///Create a GenFit track
@@ -296,13 +243,21 @@ StatusCode RecGenfitAlgDC::fitMC()
double eventStartTime=0;
if(!genfitTrack->createGenfitTrackFromEDM4HepTrack(pidType,dcTrack,
eventStartTime)){
- log<<MSG::DEBUG<<"createGenfitTrack failed!"<<endmsg;
+ debug()<<"createGenfitTrackFromEDM4HepTrack failed!"<<endmsg;
return StatusCode::SUCCESS;
}
- double sigma=0.02;//cm
- if(0==genfitTrack->addWireMeasurementOnTrack(dcTrack,sigma)){
- log<<MSG::DEBUG<<"addWireMeasurementOnTrack 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();
@@ -318,6 +273,7 @@ StatusCode RecGenfitAlgDC::fitMC()
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) {
@@ -326,11 +282,20 @@ StatusCode RecGenfitAlgDC::fitMC()
}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;
}
@@ -339,11 +304,11 @@ StatusCode RecGenfitAlgDC::fitMC()
StatusCode RecGenfitAlgDC::finalize()
{
MsgStream log(msgSvc(), name());
- log << MSG::INFO << " RecGenfitAlgDC in finalize()" << endmsg;
+ info()<< "RecGenfitAlgDC in finalize()" << endmsg;
m_genfitFitter->writeHist();
delete m_genfitFitter;
- std::cout<<"RecGenfitAlgDC nRecTrack="<<m_nDCTrack<<" success e "
+ 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){
@@ -360,12 +325,15 @@ void RecGenfitAlgDC::debugTrack(int pidType,const GenfitTrack* genfitTrack)
{
/// Get fit status
const genfit::FitStatus* fitState = genfitTrack->getFitStatus();
- double chi2 = fitState->getChi2();
- double ndf = fitState->getNdf();
- double charge = fitState->getCharge();
- int isFitted = fitState->isFitted();
- int isConverged = fitState->isFitConverged();
- int isConvergedFully = fitState->isFitConvergedFully();
+ 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;
@@ -373,11 +341,35 @@ void RecGenfitAlgDC::debugTrack(int pidType,const GenfitTrack* genfitTrack)
TVector3 fittedMom;
int fittedState=genfitTrack->getFittedState(fittedPos,fittedMom,fittedCov);
HelixClass helix;//mm and GeV
- float pos[3]={fittedPos.X()/dd4hep::mm,fittedPos.Y()/dd4hep::mm,fittedPos.Z()/dd4hep::mm};
- float mom[3]={fittedMom.X(),fittedMom.Y(),fittedMom.Z()};
+ 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()
@@ -385,9 +377,11 @@ 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_nHitMdc=simDCHitCol->size();
+ m_nSimDCHit=simDCHitCol->size();
int iMcParticle=0;
int iHit=0;
for(auto mcParticle : *mcParticleCol){
@@ -406,7 +400,7 @@ void RecGenfitAlgDC::debugEvent()
float px=mcPocaMom.x;
float py=mcPocaMom.y;
float pz=mcPocaMom.z;
- std::cout<<__FILE__<<" "<<px<<" "<<py<<" "<<pz<<std::endl;
+ debug()<<" "<<px<<" "<<py<<" "<<pz<<endmsg;
m_pocaMomMcP[iMcParticle]=sqrt(px*px+py*py+pz*pz);
iMcParticle++;
}
diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h
index 8b2754f7226e9ff19d20e55f97bafbb6ba4df8bb..14922bc694d5c6909232a73795d9afbfd2338630 100644
--- a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h
+++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h
@@ -3,7 +3,7 @@
/// This is an algorithm for track fitting for CEPC track with genfit.
///
/// In this file, including:
-/// An algorithm for track fitting with genfit with 5 pid hypothesis
+/// An algorithm for drift chamber track fitting with genfit with 5 hypothesis
///
/// Units are following DD4hepUnits
///
@@ -61,17 +61,14 @@ class RecGenfitAlgDC:public GaudiAlgorithm {
GenfitFitter* m_genfitFitter;//The pointer to a GenfitFitter
const GenfitField* m_genfitField;//The pointer to a GenfitField
- StatusCode fitMC();
void debugTrack(int pidType,const GenfitTrack* genfitTrack);
void debugEvent();
- //void debugDoca(RecMdcTrack* recMdcTrack,const GenfitTrack* genfitTrack);
- //void debugTruthHit();
DataHandle<edm4hep::EventHeaderCollection> _headerCol{
"EventHeaderCol", Gaudi::DataHandle::Reader, this};
//Drift chamber rec hit and trac
DataHandle<edm4hep::TrackerHitCollection> m_digiDCHitsCol{
- "DigiDCHitsCollection", Gaudi::DataHandle::Reader, this};
+ "DigiDCHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::TrackCollection> m_dcTrackCol{
"DCTrackCollection", Gaudi::DataHandle::Reader, this};
//Mc truth
@@ -79,6 +76,10 @@ class RecGenfitAlgDC:public GaudiAlgorithm {
"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};
@@ -91,32 +92,39 @@ class RecGenfitAlgDC:public GaudiAlgorithm {
dd4hep::Detector* m_dd4hep;
dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber;
dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
- Gaudi::Property<std::string> m_readout_name{this, "readout",
- "DriftChamberHitsCollection"};
- int m_inputType;
- int m_debug;
- std::string m_fitterType;//default DAFRef, DAFRef, DAF
- bool m_correctBremsstrahlung;//defalut do not correct bremsstrahlung
- bool m_noMaterialEffects;//defalut false, correct material effects
- int m_maxIteration;//default 20
- bool m_resortHits;//default true
- double m_bStart;
- double m_bFinal;
- double m_dcCornerCuts;
- int m_ndfCut;//default true
- double m_chi2Cut;//default 1000
+ 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;
- int m_debugPid;
- bool m_useTruthSeed;
- bool m_useRecLRAmbig;
+ //bool m_useRecLRAmbig;
- std::string m_genfitHistRootName;
- bool m_firstPid;
- bool m_useTruthHit;
- bool m_showDisplay;
- double m_initCovResPos;
- double m_initCovResMom;
genfit::EventDisplay* m_genfitDisplay;
clock_t m_timer;
@@ -126,6 +134,7 @@ class RecGenfitAlgDC:public GaudiAlgorithm {
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
@@ -140,14 +149,16 @@ class RecGenfitAlgDC:public GaudiAlgorithm {
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_nDofKal;
+ NTuple::Array<int> m_chargeKal;
NTuple::Array<double> m_chi2Kal;
- NTuple::Array<bool> m_convergeKal;
- NTuple::Array<bool> m_isFittedKal;
+ 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_nHitMdc;
- NTuple::Item<int> m_nClusterCgem;
+ 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;
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/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