From 18e0929e510d92b62b8b9519d97045d881aaa40a Mon Sep 17 00:00:00 2001
From: Zhang Yao <zhangyao@ihep.ac.cn>
Date: Tue, 1 Dec 2020 14:27:16 +0800
Subject: [PATCH] Add RecGenfitAlg with wire measurements
---
Examples/options/tut_detsim_digi_fit_STD.py | 210 +++++
Reconstruction/RecGenfitAlg/CMakeLists.txt | 39 +
Reconstruction/RecGenfitAlg/README.md | 12 +
.../RecGenfitAlg/src/GenfitField.cpp | 59 ++
Reconstruction/RecGenfitAlg/src/GenfitField.h | 38 +
.../RecGenfitAlg/src/GenfitFitter.cpp | 614 +++++++++++++++
.../RecGenfitAlg/src/GenfitFitter.h | 200 +++++
.../src/GenfitMaterialInterface.cpp | 353 +++++++++
.../src/GenfitMaterialInterface.h | 74 ++
Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp | 19 +
Reconstruction/RecGenfitAlg/src/GenfitMsg.h | 26 +
.../RecGenfitAlg/src/GenfitTrack.cpp | 733 ++++++++++++++++++
Reconstruction/RecGenfitAlg/src/GenfitTrack.h | 199 +++++
.../RecGenfitAlg/src/RecGenfitAlgDC.cpp | 416 ++++++++++
.../RecGenfitAlg/src/RecGenfitAlgDC.h | 179 +++++
Reconstruction/Tracking/CMakeLists.txt | 12 +-
.../src/TruthTracker/TruthTrackerAlg.cpp | 203 +++++
.../src/TruthTracker/TruthTrackerAlg.h | 79 ++
Utilities/DataHelper/src/HelixClass.cc | 2 +
19 files changed, 3463 insertions(+), 4 deletions(-)
create mode 100644 Examples/options/tut_detsim_digi_fit_STD.py
create mode 100644 Reconstruction/RecGenfitAlg/CMakeLists.txt
create mode 100644 Reconstruction/RecGenfitAlg/README.md
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitField.cpp
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitField.h
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitFitter.h
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitMsg.h
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
create mode 100644 Reconstruction/RecGenfitAlg/src/GenfitTrack.h
create mode 100644 Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
create mode 100644 Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h
create mode 100644 Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
create mode 100644 Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
diff --git a/Examples/options/tut_detsim_digi_fit_STD.py b/Examples/options/tut_detsim_digi_fit_STD.py
new file mode 100644
index 00000000..9dd7f4bb
--- /dev/null
+++ b/Examples/options/tut_detsim_digi_fit_STD.py
@@ -0,0 +1,210 @@
+#!/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 DETCEPCV4ROOT." )
+ 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 = [90] # rad; 0deg
+gun.PhiMaxs = [90] # 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")
+dCHDigiAlg.readout = "DriftChamberHitsCollection"
+dCHDigiAlg.drift_velocity = 40#um/ns
+dCHDigiAlg.mom_threshold = 0 #GeV
+dCHDigiAlg.SimDCHitCollection = "DriftChamberHitsCollection"
+dCHDigiAlg.DigiDCHitCollection = "DigiDCHitsCollection"
+dCHDigiAlg.AssociationCollection = "DCHAssociationCollectio"
+dCHDigiAlg.WriteAna = True
+
+##############################################################################
+# POD I/O
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "fit_DCH.root"
+out.outputCommands = ["keep *"]
+
+##############################################################################
+# TruthTrackerAlg
+##############################################################################
+from Configurables import TruthTrackerAlg
+truthTrackerAlg = TruthTrackerAlg("TruthTrackerAlg")
+truthTrackerAlg.DCHitAssociationCollection="DCHAssociationCollectio"
+truthTrackerAlg.debug = 1
+
+##############################################################################
+# GenfitAlg
+##############################################################################
+from Configurables import RecGenfitAlgDC
+recGenfitAlgDC = RecGenfitAlgDC("RecGenfitAlgDC")
+recGenfitAlgDC.debug = 11
+
+##############################################################################
+# NTupleSvc
+##############################################################################
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = ["MyTuples DATAFILE='DCH_digi_ana.root' OPT='NEW' TYP='ROOT'",
+ "TruthTrackerAlg DATAFILE='result_truthTracker_STD.root' OPT='NEW' TYP='ROOT'",
+ "RecGenfitAlgDC DATAFILE='result_fit_SDT.root' OPT='NEW' TYP='ROOT'"]
+
+##############################################################################
+# ApplicationMgr
+##############################################################################
+
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, truthTrackerAlg,
+ recGenfitAlgDC, out],
+ EvtSel = 'NONE',
+ EvtMax = 1,
+ ExtSvc = [rndmengine, dsvc, geosvc, ntsvc],
+ HistogramPersistency = "ROOT",
+ OutputLevel=DEBUG
+)
diff --git a/Reconstruction/RecGenfitAlg/CMakeLists.txt b/Reconstruction/RecGenfitAlg/CMakeLists.txt
new file mode 100644
index 00000000..920ae7b9
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/CMakeLists.txt
@@ -0,0 +1,39 @@
+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
+)
+
+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
+ )
+
+# 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
+ )
diff --git a/Reconstruction/RecGenfitAlg/README.md b/Reconstruction/RecGenfitAlg/README.md
new file mode 100644
index 00000000..3b08f667
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/README.md
@@ -0,0 +1,12 @@
+# RecGenfitAlg package
+
+This package is the track fitting algorithms to do track fitting with genfit.
+It's incudingg the interface to genfit fitter, track and measurement:
+* GenfitAlg: The algorithm for drift chamber fitting with genfit
+* GenfitFitter: fitter implementation
+* GenfitTrack: class for create and access of genfit Track and Measurement
+* GenfitField: implementation of genfit AbsField
+* GenfitGeoMaterialInterface: implementation of genfit AbsMaterialInterface
+* GenfitMsg: a class to get Gaudi message instance in Genfit classes.
+* HelixTrackModel: a class to get Gaudi message instance in Genfit classes.
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
new file mode 100644
index 00000000..de6df542
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
@@ -0,0 +1,59 @@
+#include "GenfitField.h"
+#include "GenfitMsg.h"
+
+//External
+#include "DD4hep/DD4hepUnits.h"
+
+//genfit
+#include "FieldManager.h"
+#include "ConstField.h"
+
+///-------------------------------------------------------------
+/// GenfitField constructor
+///-------------------------------------------------------------
+GenfitField::GenfitField(dd4hep::OverlayedField dd4hepField)
+:m_dd4hepField(dd4hepField)
+{
+ genfit::FieldManager::getInstance()->init(this);
+}
+
+
+///-------------------------------------------------------------
+/// GenfitField get field value by TVector3
+///-------------------------------------------------------------
+/// unit in genfit is cm and kGauss
+TVector3 GenfitField::get(const TVector3& pos) const
+{
+ double B[3]={1e9,1e9,1e9};
+ get(pos.X(),pos.Y(),pos.Z(),B[0],B[1],B[2]);
+ return TVector3(B[0],B[1],B[2]);
+}
+
+///-------------------------------------------------------------
+/// GenfitField get field value by double
+///-------------------------------------------------------------
+/// unit in genfit for position is cm and Bxyz is kGauss
+/// unit in DD4hepUnits kilogaus and cm? FIXME
+void
+GenfitField::get(const double& posX, const double& posY, const double& posZ,
+ double& Bx, double& By, double& Bz) const
+{
+ /// get field from dd4hepField
+ const dd4hep::Direction& field=m_dd4hepField.magneticField(
+ dd4hep::Position(posX/dd4hep::cm,posY/dd4hep::cm,posZ/dd4hep::cm));
+ //m_dd4hepField.magneticField(pos,B);
+
+ double factor = 1./dd4hep::kilogauss;
+ Bx=field.X()/dd4hep::kilogauss;
+ By=field.Y()/dd4hep::kilogauss;
+ Bz=field.Z()/dd4hep::kilogauss;
+
+// std::cout<<"GenfitField "
+// <<Form("xyz(%f,%f,%f)cm B(%f,%f,%f)kilogauss",posX,posY,posZ,Bx,By,Bz)
+// <<std::endl;
+}
+
+double GenfitField::getBz(const TVector3& pos) const
+{
+ return get(pos).Z();
+}
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitField.h b/Reconstruction/RecGenfitAlg/src/GenfitField.h
new file mode 100644
index 00000000..406d0e6b
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitField.h
@@ -0,0 +1,38 @@
+///////////////////////////////////////////////////////////
+//// An implementation of genfit AbsBField - GenfitField
+//// Authors:
+//// Yao ZHANG (zhangyao@ihep.ac.cn)
+///////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_GENFITFIELD_H
+#define RECGENFITALG_GENFITFIELD_H
+
+#include "AbsBField.h"
+#include "DD4hep/Fields.h"
+
+class TVector3;
+
+
+/// Calss for translating Field into genfit::AbsBField
+class GenfitField : public genfit::AbsBField{
+ public:
+ GenfitField(dd4hep::OverlayedField dd4hepField);
+ virtual ~GenfitField(){;}
+
+ //Get field value from TVector3
+ TVector3 get(const TVector3& pos) const override;
+
+ //Get field value from doubles
+ void get(const double& posX, const double& posY, const double& posZ,
+ double& Bx, double& By, double& Bz) const override;
+
+ //Get Bz, kilogauss
+ double getBz(const TVector3& pos) const;
+
+ //Get dd4hep field
+ const dd4hep::OverlayedField field() const {return m_dd4hepField;}
+
+ private:
+ dd4hep::OverlayedField m_dd4hepField;
+};
+#endif
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
new file mode 100644
index 00000000..2a7745b6
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
@@ -0,0 +1,614 @@
+#include "GenfitFitter.h"
+#include "GenfitTrack.h"
+#include "GenfitMsg.h"
+#include "GenfitField.h"
+#include "GenfitMaterialInterface.h"
+
+//Gaudi
+#include "GaudiKernel/StatusCode.h"
+#include "GaudiKernel/ISvcLocator.h"
+
+//External
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/Detector.h"
+
+//genfit
+#include <Track.h>
+#include <Exception.h>
+#include <FieldManager.h>
+#include <TGeoMaterialInterface.h>
+#include <TGeoManager.h>
+#include <MaterialEffects.h>
+#include <MeasuredStateOnPlane.h>
+#include <KalmanFitter.h>
+#include <KalmanFitterRefTrack.h>
+#include <StateOnPlane.h>
+#include <DAF.h>
+#include <AbsKalmanFitter.h>
+#include <KalmanFitterInfo.h>
+
+//ROOT
+#include <TVector3.h>
+#include <TGeoManager.h>
+
+//STL
+#include <iostream>
+#include <string>
+
+
+GenfitFitter::~GenfitFitter(){
+ delete m_absKalman;
+}
+
+GenfitFitter::GenfitFitter(const char* type, const char* name):
+ m_absKalman(nullptr)
+ ,m_genfitField(nullptr)
+ ,m_geoMaterial(nullptr)
+ ,m_fitterType(type)
+ ,m_name(name)
+ ,m_minIterations(4)
+ ,m_maxIterations(10)
+ ,m_deltaPval(1e-3)
+ ,m_relChi2Change(0.2)
+ ,m_blowUpFactor(1e3)
+ ,m_resetOffDiagonals(true)
+ ,m_blowUpMaxVal(1.e6)
+ ,m_multipleMeasurementHandling(genfit::unweightedClosestToPredictionWire)
+ ,m_maxFailedHits(-1)
+ ,m_deltaWeight(1e-3)
+ ,m_annealingBetaStart(100)
+ ,m_annealingBetaStop(0.01)
+ ,m_annealingNSteps(0.01)
+ ,m_noEffects(false)
+ ,m_energyLossBetheBloch(true)
+ ,m_noiseBetheBloch(true)
+ ,m_noiseCoulomb(true)
+ ,m_energyLossBrems(false)
+ ,m_noiseBrems(false)
+ ,m_ignoreBoundariesBetweenEqualMaterials(true)
+ ,m_mscModelName("GEANE")
+ ,m_debug(0)
+ ,m_hist(0)
+{
+ /// Initialize genfit fitter
+ init();
+}
+
+void GenfitFitter::setField(const GenfitField* field)
+{
+ if(nullptr==m_genfitField) m_genfitField=field;
+}
+
+/// Set geometry for material, use geometry from IOADatabase
+void GenfitFitter::setGeoMaterial(const dd4hep::Detector* dd4hepGeo)
+{
+ if(nullptr==m_geoMaterial){
+ m_geoMaterial=GenfitMaterialInterface::getInstance(dd4hepGeo);
+ }
+}
+
+/// initialize genfit fitter
+int GenfitFitter::init(bool deleteOldFitter)
+{
+ if(deleteOldFitter && m_absKalman) delete m_absKalman;
+
+ GenfitMsg::get()<<MSG::DEBUG<<"Initialize GenfitFitter with "
+ <<m_fitterType<<endmsg;
+
+ if (m_fitterType=="DAFRef") {
+ m_absKalman = new genfit::DAF(true,getDeltaPval(),
+ getConvergenceDeltaWeight());
+ }
+ else if (m_fitterType=="DAF") {
+ m_absKalman = new genfit::DAF(false,getDeltaPval(),
+ getConvergenceDeltaWeight());
+ }
+ else if (m_fitterType=="KalmanFitter") {
+ m_absKalman = new genfit::KalmanFitter(getMaxIterations());
+ }
+ else if (m_fitterType=="KalmanFitterRefTrack") {
+ m_absKalman = new genfit::KalmanFitterRefTrack(getMaxIterations());
+ }
+ else {
+ m_absKalman = nullptr;
+ GenfitMsg::get()<< MSG::DEBUG<<"Fitter type is invalid:"
+ <<m_fitterType<<endmsg;
+ return -1;
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<"Fitter type is "<<m_fitterType<<endmsg;
+std::cout<<__FILE__<<" "<<__LINE__<<m_absKalman<<std::endl;
+ m_absKalman->setDebugLvl(m_debug);
+std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
+
+ return 0;
+}
+
+/// Fit a track from a candidate track
+int GenfitFitter::processTrackWithRep(GenfitTrack* track,int repID,bool resort)
+{
+ GenfitMsg::get()<<MSG::DEBUG<< "In ProcessTrackWithRep rep "<<repID<<endmsg;
+ if(getDebug()>2) print("");
+
+ if(getDebug()>0){
+ GenfitMsg::get()<<MSG::DEBUG<<"Print track seed "<<endmsg;
+ track->getTrack()->getStateSeed().Print();
+ }
+ /// Do the fitting
+ try{
+ m_absKalman->processTrackWithRep(track->getTrack(), track->getRep(repID),
+ resort);
+ }catch(genfit::Exception& e){
+ GenfitMsg::get()<<MSG::DEBUG<<"Genfit exception caught "<<endmsg;
+ e.what();
+ return false;
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<"End of ProcessTrackWithRep"<<endmsg;
+ return true;
+} // End of ProcessTrackWithRep
+
+/// Fit a track from a candidate track
+int GenfitFitter::processTrack(GenfitTrack* track, bool resort)
+{
+ GenfitMsg::get()<<MSG::DEBUG<<"In ProcessTrack"<<endmsg;
+ if(getDebug()>2) print("");
+
+ /// Do the fitting
+ try{
+ m_absKalman->processTrack(track->getTrack(),resort);
+ }catch(genfit::Exception& e){
+ GenfitMsg::get()<<MSG::DEBUG<<"Genfit exception caught "<<endmsg;
+ e.what();
+ return false;
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<"End of ProcessTrack"<<endmsg;
+ return true;
+} // End of ProcessTrack
+
+
+void GenfitFitter::setFitterType(const char* val)
+{
+ std::string oldFitterType=m_fitterType;
+ m_fitterType = val;
+ GenfitMsg::get()<<MSG::DEBUG<<"Fitter type is "<<m_fitterType<<endmsg;
+ init(oldFitterType==val);
+}
+
+/// Extrapolate track to the cloest point of approach(POCA) to the wire of hit,
+/// return StateOnPlane of this POCA
+/// inputs
+/// pos,mom ... position & momentum at starting point, unit= [mm]&[GeV/c]
+/// (converted to [cm]&[GeV/c] inside this function)
+/// hit ... destination
+/// outputs poca [mm] (converted from [cm] in this function) ,global
+double GenfitFitter::extrapolateToHit( TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, const GenfitTrack* track,
+ TVector3 pos, TVector3 mom,
+ TVector3 end0,//one end of the hit wire
+ TVector3 end1,//the orhter end of the hit wire
+ int debug,
+ int repID,
+ bool stopAtBoundary,
+ bool calcJacobianNoise)
+{
+
+ return track->extrapolateToHit(poca,pocaDir,pocaOnWire,doca,pos,mom,end0,
+ end1,debug,repID,stopAtBoundary,calcJacobianNoise);
+}//end of ExtrapolateToHit
+
+
+/// Extrapolate the track to the cyliner at fixed raidus
+/// position & momentum as starting point
+/// position and momentum at global coordinate in dd4hepUnit
+/// return trackLength in dd4hepUnit
+ double
+GenfitFitter::extrapolateToCylinder(TVector3& pos, TVector3& mom,
+ GenfitTrack* track, double radius, const TVector3 linePoint,
+ const TVector3 lineDirection, int hitID, int repID,
+ bool stopAtBoundary, bool calcJacobianNoise)
+{
+ double trackLength(1e9*dd4hep::cm);
+ if(!track->getFitStatus(repID)->isFitted()) return trackLength;
+ try{
+ // get track rep
+ genfit::AbsTrackRep* rep = track->getRep(repID);
+ if(nullptr == rep) {
+ GenfitMsg::get()<<MSG::DEBUG<<"In ExtrapolateToCylinder rep is null"
+ <<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+
+ // get track point with fitter info
+ genfit::TrackPoint* tp =
+ track->getTrack()->getPointWithFitterInfo(hitID,rep);
+ if(nullptr == tp) {
+ GenfitMsg::get()<<MSG::DEBUG<<
+ "In ExtrapolateToCylinder TrackPoint is null"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+
+ // get fitted state on plane of this track point
+ genfit::KalmanFittedStateOnPlane* state =
+ static_cast<genfit::KalmanFitterInfo*>(
+ tp->getFitterInfo(rep))->getBackwardUpdate();
+
+ if(nullptr == state){
+ GenfitMsg::get()<<MSG::DEBUG<<"In ExtrapolateToCylinder "
+ <<"no KalmanFittedStateOnPlane in backwardUpdate"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+ rep->setPosMom(*state, pos*(1/dd4hep::cm), mom*(1/dd4hep::GeV));
+
+ /// extrapolate
+ trackLength = rep->extrapolateToCylinder(*state,
+ radius/dd4hep::cm, linePoint*(1/dd4hep::cm), lineDirection,
+ stopAtBoundary, calcJacobianNoise);
+ // get pos&mom at extrapolated point on the cylinder
+ rep->getPosMom(*state,pos,mom);//FIXME exception exist
+ pos = pos*dd4hep::cm;
+ mom = mom*dd4hep::GeV;
+ } catch(genfit::Exception& e){
+ GenfitMsg::get() << MSG::ERROR
+ <<"Exception in GenfitFitter::extrapolateToCylinder "
+ << e.what()<<endmsg;
+ trackLength = 1e9*dd4hep::cm;
+ }
+ return trackLength*dd4hep::cm;
+}
+
+double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom,
+ const GenfitTrack* track,
+ const TVector3& point,
+ int repID,// same with pidType
+ bool stopAtBoundary,
+ bool calcJacobianNoise) const
+{
+ double trackLength(1e9*dd4hep::cm);
+ if(!track->getFitStatus(repID)->isFitted()) return trackLength;
+ try{
+ // get track rep
+ genfit::AbsTrackRep* rep = track->getRep(repID);
+ if(nullptr == rep) {
+ GenfitMsg::get()<<MSG::DEBUG<<"In ExtrapolateToPoint rep "
+ <<repID<<" not exist!"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+
+ /// extrapolate to point
+ //genfit::StateOnPlane state(*(&(track->getTrack()->getFittedState(0,rep))));
+
+ // get track point with fitter info
+ genfit::TrackPoint* tp =
+ track->getTrack()->getPointWithFitterInfo(0,rep);
+ if(nullptr == tp) {
+ GenfitMsg::get()<<MSG::DEBUG<<
+ "In ExtrapolateToPoint TrackPoint is null"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+
+ // get fitted state on plane of this track point
+ genfit::KalmanFittedStateOnPlane* state =
+ static_cast<genfit::KalmanFitterInfo*>(
+ tp->getFitterInfo(rep))->getBackwardUpdate();
+
+ if(nullptr == state) {
+ GenfitMsg::get()<<MSG::DEBUG<<
+ "In ExtrapolateToPoint KalmanFittedStateOnPlane is null"<<endmsg;
+ return trackLength*dd4hep::cm;
+ }
+ trackLength = rep->extrapolateToPoint(*state,
+ point*(1/dd4hep::cm),stopAtBoundary, calcJacobianNoise);
+ rep->getPosMom(*state,pos,mom);//FIXME exception exist
+ pos = pos*dd4hep::cm;
+ mom = mom*dd4hep::GeV;
+ } catch(genfit::Exception& e){
+ GenfitMsg::get() << MSG::ERROR
+ <<"Exception in GenfitFitter::extrapolateToPoint"
+ << e.what()<<endmsg;
+ trackLength = 1e9*dd4hep::cm;
+ }
+ return trackLength*dd4hep::cm;
+}//end of ExtrapolateToPoint
+
+GenfitFitter& GenfitFitter::operator=(
+ const GenfitFitter& r)
+{
+ m_fitterType = r.m_fitterType;
+ m_minIterations = r.m_minIterations;
+ m_maxIterations = r.m_maxIterations;
+ m_deltaPval = r.m_deltaPval;
+ m_relChi2Change = r.m_relChi2Change;
+ m_blowUpFactor = r.m_blowUpFactor;
+ m_resetOffDiagonals = r.m_resetOffDiagonals;
+ m_blowUpMaxVal = r.m_blowUpMaxVal;
+ m_multipleMeasurementHandling = r.m_multipleMeasurementHandling;
+ m_maxFailedHits = r.m_maxFailedHits;
+ m_annealingNSteps = r.m_annealingNSteps;
+ m_deltaWeight = r.m_deltaWeight;
+ m_annealingBetaStart = r.m_annealingBetaStart;
+ m_annealingBetaStop = r.m_annealingBetaStop;
+ m_noEffects = r.m_noEffects;
+ m_energyLossBetheBloch= r.m_energyLossBetheBloch;
+ m_noiseBetheBloch = r.m_noiseBetheBloch;
+ m_noiseCoulomb = r.m_noiseCoulomb;
+ m_energyLossBrems = r.m_energyLossBrems;
+ m_noiseBrems = r.m_noiseBrems;
+ m_ignoreBoundariesBetweenEqualMaterials
+ = r.m_ignoreBoundariesBetweenEqualMaterials;
+ m_mscModelName = r.m_mscModelName;
+ m_debug = r.m_debug;
+ m_hist = r.m_hist;
+ return *this;
+}
+
+void GenfitFitter::print(const char* name)
+{
+
+ if(0==strlen(name)) name=m_name;
+ //TODO print all fitting parameters
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" GenfitFitter Global Parameters:"<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" Fitter type = " << m_fitterType<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" Fitting Parameters:"<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" MinIteration = " << m_minIterations<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" MaxIteration = " << m_maxIterations<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_deltaPval = " << m_deltaPval<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" Material Effect Parameters:"<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_noEffects = " << m_noEffects<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_energyLossBetheBloch= " << m_energyLossBetheBloch<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_noiseBetheBloch = " << m_noiseBetheBloch<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_noiseCoulomb = " << m_noiseCoulomb<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_energyLossBrems = " << m_energyLossBrems<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_noiseBrems = " << m_noiseBrems<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_ignoreBoundariesBetweenEqualMaterials= "
+ << m_ignoreBoundariesBetweenEqualMaterials<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_mscModelName = " << m_mscModelName<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_debug = " << m_debug<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<name
+ <<" m_hist = " << m_hist<<endmsg;
+
+ if(m_fitterType=="DAF"||m_fitterType=="DAFRef"){
+ genfit::DAF* daf = getDAF();
+
+ if(nullptr != daf){
+ std::ostringstream sBetas;
+ std::vector<double> betas = daf->getBetas();
+ for (unsigned int i=0; i<betas.size(); ++i) {
+ sBetas<<betas.at(i)<<",";
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<" print "<<betas.size()<<
+ " betas: %s "<<sBetas.str()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<m_name
+ << " m_deltaWeight = " << m_deltaWeight<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF maxIterations= "
+ << daf->getMaxIterations()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF minIterations= "
+ << daf->getMinIterations()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF deltaPval= "
+ << daf->getDeltaPval()<<endmsg;
+ }
+ }
+ genfit::KalmanFitterRefTrack* ref = getKalRef();
+ if(nullptr != ref){
+ std::ostringstream sBetas;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF maxIterations= "
+ << ref->getMaxIterations()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF minIterations= "
+ << ref->getMinIterations()<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<" DAF deltaPval= "
+ << ref->getDeltaPval()<<endmsg;
+ }
+}
+
+///Setters of AbsKalmanFitter
+void GenfitFitter::setMinIterations(unsigned int val)
+{
+ m_absKalman->setMinIterations(val);
+ m_minIterations = val;
+}
+
+void GenfitFitter::setMaxIterations(unsigned int val)
+{
+ if(val<=4) {
+ GenfitMsg::get()<<MSG::ERROR<<"Could NOT set MaxIteration<=4!"<<endmsg;
+ }
+ m_absKalman->setMaxIterations(val);
+ m_maxIterations = val;
+}
+
+void GenfitFitter::setMaxIterationsBetas(double bStart,double bFinal,
+ unsigned int val)
+{
+ m_absKalman->setMaxIterations(val);
+ m_maxIterations = val;
+ if(val<=4) {
+ GenfitMsg::get()<<MSG::ERROR<<"Could NOT set MaxIteration<=4!"<<endmsg;
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<"bStart "<<bStart<<" bFinal "<<bFinal
+ <<" MaxIteration "<<val<<endmsg;
+ // genfit version less than 2.2.0
+ genfit::DAF* daf = dynamic_cast<genfit::DAF*> (m_absKalman);
+ daf->setAnnealingScheme(bStart,bFinal,val);
+}
+
+void GenfitFitter::setDeltaPval(double val)
+{
+ m_absKalman->setDeltaPval(val);
+ m_deltaPval = val;
+}
+
+void GenfitFitter::setRelChi2Change(double val)
+{
+ m_absKalman->setRelChi2Change(val);
+ m_relChi2Change = val;
+}
+
+void GenfitFitter::setBlowUpFactor(double val)
+{
+ m_absKalman->setBlowUpFactor(val);
+ m_blowUpFactor = val;
+}
+
+void GenfitFitter::setResetOffDiagonals(bool val)
+{
+ m_absKalman->setResetOffDiagonals(val);
+ m_resetOffDiagonals = val;
+}
+
+void GenfitFitter::setBlowUpMaxVal(double val)
+{
+ m_absKalman->setBlowUpMaxVal(val);
+ m_blowUpMaxVal = val;
+}
+
+void GenfitFitter::setMultipleMeasurementHandling(
+ genfit::eMultipleMeasurementHandling val)
+{
+ m_absKalman->setMultipleMeasurementHandling(val);
+ m_multipleMeasurementHandling = val;
+}
+
+void GenfitFitter::setMaxFailedHits(int val)
+{
+ m_absKalman->setMaxFailedHits(val);
+ m_maxFailedHits = val;
+}
+
+///DAF setters
+void GenfitFitter::setConvergenceDeltaWeight(double val)
+{
+ genfit::DAF* daf = getDAF();
+ if(nullptr != daf) daf->setConvergenceDeltaWeight(val);
+ m_deltaWeight = val;
+}
+void GenfitFitter::setAnnealingScheme(
+ double bStart, double bFinal, unsigned int nSteps)
+{
+ genfit::DAF* daf = getDAF();
+ if(nullptr != daf) daf->setAnnealingScheme(bStart, bFinal, nSteps);
+ m_annealingBetaStart = bStart;
+ m_annealingBetaStop = bFinal;
+ m_annealingNSteps = nSteps;
+}
+
+//TODO chi2Cuts?
+///Material effects
+void GenfitFitter::setNoEffects(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setNoEffects();
+ m_noEffects = val;
+}
+
+void GenfitFitter::setEnergyLossBetheBloch(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setEnergyLossBetheBloch(val);
+ m_energyLossBetheBloch = val;
+}
+
+void GenfitFitter::setNoiseBetheBloch(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setNoiseBetheBloch(val);
+ m_noiseBetheBloch = val;
+}
+
+void GenfitFitter::setNoiseCoulomb(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setNoiseCoulomb(val);
+ m_noiseCoulomb = val;
+}
+
+void GenfitFitter::setEnergyLossBrems(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setEnergyLossBrems(val);
+ m_energyLossBrems = val;
+}
+
+void GenfitFitter::setNoiseBrems(bool val)
+{
+ genfit::MaterialEffects::getInstance()->setNoiseBrems(val);
+ m_noiseBrems = val;
+}
+
+void GenfitFitter::setIgnoreBoundariesBetweenEqualMaterials(bool val)
+{
+ genfit::MaterialEffects::getInstance()->
+ ignoreBoundariesBetweenEqualMaterials(val);
+ m_ignoreBoundariesBetweenEqualMaterials = val;
+}
+
+void GenfitFitter::setMscModelName(std::string val)
+{
+ genfit::MaterialEffects::getInstance()->setMscModel(val);
+ m_mscModelName = val;
+}
+
+void GenfitFitter::setMaterialDebugLvl(unsigned int val)
+{
+ genfit::MaterialEffects::getInstance()->setDebugLvl(val);
+}
+
+///Set GenfitFitter parameters
+void GenfitFitter::setDebug(unsigned int val)
+{
+ GenfitMsg::get()<<MSG::DEBUG<<"set fitter debugLvl "<<val<<endmsg;
+ m_absKalman->setDebugLvl(val);
+ if(nullptr!=getDAF()) getDAF()->setDebugLvl(val);
+ if(val>10) genfit::MaterialEffects::getInstance()->setDebugLvl(val);
+ m_debug = val;
+}
+
+void GenfitFitter::setHist(unsigned int val) {m_hist = val;}
+
+genfit::DAF* GenfitFitter::getDAF()
+{
+ genfit::DAF* daf = nullptr;
+ try{
+ daf = dynamic_cast<genfit::DAF*> (m_absKalman);
+ }catch(...){
+ GenfitMsg::get()<< MSG::ERROR
+ << "dynamic_cast m_rom AbsFitter to DAF m_ailed!"<<endmsg;
+ return nullptr;
+ }
+ return daf;
+}
+
+genfit::KalmanFitterRefTrack* GenfitFitter::getKalRef()
+{
+ genfit::KalmanFitterRefTrack* ref=nullptr;
+ try{
+ ref = dynamic_cast<genfit::KalmanFitterRefTrack*> (m_absKalman);
+ }catch(...){
+ GenfitMsg::get()<< MSG::ERROR
+ << "dynamic_cast m_rom AbsFitter to KalmanFitterRefTrack m_ailed!"
+ <<endmsg;
+ }
+ return ref;
+}
+
+void GenfitFitter::initHist(std::string name)
+{
+ GenfitMsg::get()<<MSG::DEBUG<<"GenfitFitter::initHist "<<name<<endmsg;
+ //getDAF()->initHist(name);
+}
+
+void GenfitFitter::writeHist()
+{
+ GenfitMsg::get()<<MSG::DEBUG<<"GenfitFitter::writeHist "<<endmsg;
+ //getDAF()->writeHist();
+}
+
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitFitter.h b/Reconstruction/RecGenfitAlg/src/GenfitFitter.h
new file mode 100644
index 00000000..9ca46644
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitFitter.h
@@ -0,0 +1,200 @@
+//////////////////////////////////////////////////////////////////
+///
+/// This is an interface of to call genfit fitter
+///
+/// In this file, including:
+/// a genfit fitter class
+/// Fitter type is fixed after creation
+/// Magnetic field and geometry for material effect are singletons in genfit.
+/// They will be set before fitting and before each time call track rep.
+///
+/// Authors:
+/// Yao ZHANG (zhangyao@ihep.ac.cn)
+///
+//////////////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_GENFITFITTER_H
+#define RECGENFITALG_GENFITFITTER_H
+
+#include "AbsKalmanFitter.h"
+#include <string>
+
+class GenfitField;
+class GenfitMaterialInterface;
+class GenfitTrack;
+class TVector3;
+class TGeoManager;
+class ISvcLocator;
+namespace genfit{
+ class AbsKalmanFitter;
+ class KalmanFitterRefTrack;
+ class DAF;
+}
+namespace dd4hep{
+ class OverlayedField;
+ class Detector;
+}
+
+
+/// The GenfitTrack class
+class GenfitFitter{
+ public:
+ /// Type of fitters are :DAFRef,DAF,KalmanFitter,KalmanFitterRefTrack
+ GenfitFitter(const char* type="DAFRef",const char* name="GenfitFitter");
+ virtual ~GenfitFitter();
+
+ /// Magnetic field and geometry for material effect in genfit
+ /// please SET before use !!!!
+ void setField(const GenfitField* field);
+ /// please SET before use !!!!
+ void setGeoMaterial(const dd4hep::Detector* dd4hepGeo);
+
+ /// Main fitting function
+ int processTrack(GenfitTrack* track, bool resort=false);
+
+ /// fitting with rep
+ int processTrackWithRep(GenfitTrack* track,int repID=0,
+ bool resort=false);
+
+ /// Extrapolate the track to the CDC hit
+ /// Output: poca pos and dir and poca distance to the hit wire
+ /// Input: genfit track, pos and mom, two ends of a wire
+ /// pos, and mom are position & momentum at starting point
+ double extrapolateToHit(TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, const GenfitTrack* track,
+ TVector3 pos, TVector3 mom, TVector3 end0, TVector3 end1,
+ int debug=0, int repID=0, bool stopAtBoundary=false,
+ bool calcJacobianNoise=false);
+
+ /// Extrapolate the track to the cyliner at fixed raidus
+ /// Output: pos and mom at fixed radius
+ /// Input: genfitTrack, radius of cylinder at center of the origin,
+ /// repID, stopAtBoundary and calcAverageState
+ double extrapolateToCylinder(TVector3& pos, TVector3& mom,
+ GenfitTrack* track, double radius, const TVector3 linePoint,
+ const TVector3 lineDirection, int hitID =0, int repID=0,
+ bool stopAtBoundary=false, bool calcJacobianNoise=false);
+
+ /// Extrapolate the track to the point
+ /// Output: pos and mom of POCA point to point
+ /// Input: genfitTrack,point,repID,stopAtBoundary and calcAverageState
+ /// repID same with pidType
+ double extrapolateToPoint(TVector3& pos, TVector3& mom,
+ const GenfitTrack* genfitTrack, const TVector3& point,
+ int repID=0, bool stopAtBoundary = false,
+ bool calcJacobianNoise = false) const;
+
+ /// setters of fitter properties
+ void setFitterType(const char* val);
+ void setMinIterations(unsigned int val);
+ void setMaxIterations(unsigned int val);
+ void setMaxIterationsBetas(double bStart,double bFinal,unsigned int val);
+ void setDeltaPval(double val);
+ void setRelChi2Change(double val);
+ void setBlowUpFactor(double val);
+ void setResetOffDiagonals(bool val);
+ void setBlowUpMaxVal(double val);
+ void setMultipleMeasurementHandling(
+ genfit::eMultipleMeasurementHandling val);
+ void setMaxFailedHits(int val);
+ void setConvergenceDeltaWeight(double val);
+ void setAnnealingScheme(double bStart,double bFinal,unsigned int nSteps);
+ //TODO chi2cut?
+ void setNoEffects(bool val);
+ void setEnergyLossBetheBloch(bool val);
+ void setNoiseBetheBloch(bool val);
+ void setNoiseCoulomb(bool val);
+ void setEnergyLossBrems(bool val);
+ void setNoiseBrems(bool val);
+ void setIgnoreBoundariesBetweenEqualMaterials(bool val);
+ void setMscModelName(std::string val);
+ void setMaterialDebugLvl(unsigned int val);
+ void setDebug(unsigned int val);
+ void setHist(unsigned int val);
+
+ /// getters of fitter properties
+ std::string getFitterType() const {return m_fitterType;}
+ unsigned int getMinIterations() const { return m_minIterations; }
+ unsigned int getMaxIterations() const { return m_maxIterations; }
+ double getDeltaPval() const { return m_deltaPval; }
+ double getRelChi2Change() const { return m_relChi2Change; }
+ double getBlowUpFactor() const { return m_blowUpFactor; }
+ double getBlowUpMaxVal() const { return m_blowUpMaxVal; }
+ bool getResetOffDiagonals() const { return m_resetOffDiagonals; }
+ genfit::eMultipleMeasurementHandling getMultipleMeasurementHandling()
+ const { return m_multipleMeasurementHandling; }
+ int getMaxFailedHits() const { return m_maxFailedHits; }
+ double getConvergenceDeltaWeight() const { return m_deltaWeight; }
+ float getAnnealingBetaStart() const {return m_annealingBetaStart;}
+ float getAnnealingBetaStop() const {return m_annealingBetaStop;}
+ bool getNoEffects(){return m_noEffects;}
+ bool getEnergyLossBetheBloch(){return m_energyLossBetheBloch;}
+ bool getNoiseBetheBloch(){return m_noiseBetheBloch;}
+ bool getNoiseCoulomb(){return m_noiseCoulomb;}
+ bool getEnergyLossBrems(){return m_energyLossBrems;}
+ bool getNoiseBrems(){return m_noiseBrems;}
+ bool getIgnoreBoundariesBetweenEqualMaterials()
+ {return m_ignoreBoundariesBetweenEqualMaterials;}
+ std::string getMscModelName(){return m_mscModelName;}
+ int getDebug() const {return m_debug;}
+ int getHist() const {return m_hist;}
+
+ /// Printer
+ void print(const char* name="");
+ void initHist(std::string name);
+ void writeHist();
+
+ const char* getName(void) const {return m_name;}
+
+ private:
+ GenfitFitter(const GenfitFitter&){};
+ GenfitFitter& operator=(const GenfitFitter&);
+
+ /// Initialze fitter and setting fitter parameter
+ int init(bool deleteOldFitter=false);
+
+ /// Get DAF fitter
+ genfit::DAF* getDAF();
+ /// Get KalmanFitterRefTrack
+ genfit::KalmanFitterRefTrack* getKalRef();
+
+ genfit::AbsKalmanFitter* m_absKalman;/// kalman fitter object
+
+ const GenfitField* m_genfitField;//pointer to genfit field
+ GenfitMaterialInterface* m_geoMaterial;//pointer to genfit geo
+
+ ///fitting method: DAFRef,DAF,KalmanFitter,KalmanFitterRefTrack
+ std::string m_fitterType;
+
+ const char* m_name;
+
+ /// control parameters of fitter
+ unsigned int m_minIterations; /// minimum number of iterations
+ unsigned int m_maxIterations; /// maximum number of iterations
+ double m_deltaPval; /// delta pVal that converged
+ double m_relChi2Change; ///
+ double m_blowUpFactor; ///
+ bool m_resetOffDiagonals; ///
+ double m_blowUpMaxVal; ///
+ genfit::eMultipleMeasurementHandling m_multipleMeasurementHandling;///
+ int m_maxFailedHits; ///
+ double m_deltaWeight; /// delta weight that converged
+ double m_annealingBetaStart;/// start beta for annealing
+ double m_annealingBetaStop; /// stop beta for annealing
+ unsigned int m_annealingNSteps; /// n steps
+
+ /// control parAmeters of material effects
+ bool m_noEffects;
+ bool m_energyLossBetheBloch;
+ bool m_noiseBetheBloch;
+ bool m_noiseCoulomb;
+ bool m_energyLossBrems;
+ bool m_noiseBrems;
+ bool m_ignoreBoundariesBetweenEqualMaterials;
+ std::string m_mscModelName;
+
+ int m_debug; /// debug
+ int m_hist; /// hist
+};
+
+#endif
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp
new file mode 100644
index 00000000..73a59316
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp
@@ -0,0 +1,353 @@
+///
+/// Authors: ZHANG Yao (zhangyao@ihep.ac.cn)
+///
+
+#include "GenfitMaterialInterface.h"
+//Gaudi
+#include "GaudiKernel/Bootstrap.h"
+#include "GaudiKernel/SmartIF.h"
+#include "GaudiKernel/ISvcLocator.h"
+//CEPCSW
+#include "DetInterface/IGeomSvc.h"
+//Externals
+#include "DD4hep/Detector.h"
+//#include "DD4hep/DD4hepUnits.h"
+//genfit
+#include "Exception.h"
+//ROOT
+#include "TGeoManager.h"
+#include "TGeoNode.h"
+#include "TGeoMedium.h"
+#include "TGeoMaterial.h"
+#include "TGeoManager.h"
+#include "MaterialEffects.h"
+#include "Material.h"
+//STL
+#include <assert.h>
+#include <math.h>
+
+GenfitMaterialInterface* GenfitMaterialInterface::m_instance = nullptr;
+
+//#define DEBUG_GENFITGEO 1
+
+double MeanExcEnergy_get(int Z);
+double MeanExcEnergy_get(TGeoMaterial*);
+
+GenfitMaterialInterface::GenfitMaterialInterface(
+ const dd4hep::Detector* dd4hepGeo)
+{
+ assert(nullptr!=dd4hepGeo);
+ m_geoManager=&(dd4hepGeo->manager());
+ assert(nullptr!=m_geoManager);
+}
+
+GenfitMaterialInterface* GenfitMaterialInterface::getInstance(
+ const dd4hep::Detector* dd4hepGeo)
+{
+ if(nullptr == m_instance){
+ m_instance = new GenfitMaterialInterface(dd4hepGeo);
+ }
+ genfit::MaterialEffects::getInstance()->init(m_instance);
+ return m_instance;
+}
+
+void GenfitMaterialInterface::destruct(){
+ delete m_instance;
+}
+
+TGeoManager* GenfitMaterialInterface::getGeoManager()
+{
+ return m_geoManager;
+}
+
+ bool
+GenfitMaterialInterface::initTrack(double posX, double posY,
+ double posZ, double dirX, double dirY, double dirZ)
+{
+ //debugLvl_ = 1;
+#ifdef DEBUG_GENFITGEO
+ std::cout << "GenfitMaterialInterface::initTrack. \n";
+ std::cout << "Pos "; TVector3(posX, posY, posZ).Print();
+ std::cout << "Dir "; TVector3(dirX, dirY, dirZ).Print();
+#endif
+
+ // Move to the new point.
+ bool result = !isSameLocation(posX, posY, posZ, kTRUE);
+
+ // Set the intended direction.
+ setCurrentDirection(dirX, dirY, dirZ);
+
+ if (debugLvl_ > 0) {
+ std::cout << " GenfitMaterialInterface::initTrack at \n";
+ std::cout << " position: "; TVector3(posX, posY, posZ).Print();
+ std::cout << " direction: "; TVector3(dirX, dirY, dirZ).Print();
+ }
+
+ return result;
+}
+
+
+genfit::Material GenfitMaterialInterface::getMaterialParameters()
+{
+ TGeoMaterial* mat =
+ getGeoManager()->GetCurrentVolume()->GetMedium()->GetMaterial();
+ //Scalar density; /// Density in g / cm^3
+ //Scalar Z; /// Atomic number
+ //Scalar A; /// Mass number in g / mol
+ //Scalar radiationLength; /// Radiation Length in cm
+ //Scalar mEE; /// Mean excitaiton energy in eV
+ //Material from CEPCSW is NOT follow the dd4hep?FIXME
+
+ //std::cout<<__FILE__<<" "<<__LINE__<<" yzhang debug material "<<std::endl;
+ //mat->Print();
+ return genfit::Material(mat->GetDensity(), mat->GetZ(), mat->GetA(),
+ mat->GetRadLen(), MeanExcEnergy_get(mat));
+}
+
+ double
+GenfitMaterialInterface::findNextBoundary(const genfit::RKTrackRep* rep,
+ const genfit::M1x7& stateOrig,
+ double sMax, // signed
+ bool varField)
+{
+ //debugLvl_ = 1;
+ // cm, distance limit beneath which straight-line steps are taken.
+ const double delta(1.E-2);
+ const double epsilon(1.E-1); // cm, allowed upper bound on arch
+ // deviation from straight line
+
+ genfit::M1x3 SA;
+ genfit::M1x7 state7, oldState7;
+ oldState7 = stateOrig;
+
+ int stepSign(sMax < 0 ? -1 : 1);
+
+ double s = 0; // trajectory length to boundary
+
+ const unsigned maxIt = 300;
+ unsigned it = 0;
+
+ // Initialize the geometry to the current location (set by caller).
+ findNextBoundary(fabs(sMax) - s);
+ double safety = getSafeDistance(); // >= 0
+ double slDist = getStep();
+
+ // this should not happen, but it happens sometimes.
+ // The reason are probably overlaps in the geometry.
+ // Without this "hack" many small steps with size of minStep will be made,
+ // until eventually the maximum number of iterations are exceeded
+ // and the extrapolation fails
+ if (slDist < safety)
+ slDist = safety;
+ double step = slDist;
+
+ if (debugLvl_ > 0)
+ std::cout << " safety = " << safety << "; step, slDist = " << step << "\n";
+
+ while (1) {
+ if (++it > maxIt){
+ genfit::Exception exc("GenfitMaterialInterface::findNextBoundary ==> maximum number of iterations exceeded",__LINE__,__FILE__);
+ exc.setFatal();
+ throw exc;
+ }
+
+ // No boundary in sight?
+ if (s + safety > fabs(sMax)) {
+ if (debugLvl_ > 0)
+ std::cout << " next boundary is further away than sMax \n";
+ return stepSign*(s + safety); //sMax;
+ }
+
+ // Are we at the boundary?
+ if (slDist < delta) {
+ if (debugLvl_ > 0)
+ std::cout << " very close to the boundary -> return"
+ << " stepSign*(s + slDist) = "
+ << stepSign << "*(" << s + slDist << ")\n";
+ return stepSign*(s + slDist);
+ }
+
+ // We have to find whether there's any boundary on our path.
+
+ // Follow curved arch, then see if we may have missed a boundary.
+ // Always propagate complete way from original start to avoid
+ // inconsistent extrapolations.
+ state7 = stateOrig;
+ rep->RKPropagate(state7, nullptr, SA, stepSign*(s + step), varField);
+
+ // Straight line distance² between extrapolation finish and
+ // the end of the previously determined safe segment.
+ double dist2 = (pow(state7[0] - oldState7[0], 2)
+ + pow(state7[1] - oldState7[1], 2)
+ + pow(state7[2] - oldState7[2], 2));
+ // Maximal lateral deviation².
+ double maxDeviation2 = 0.25*(step*step - dist2);
+
+ if (step > safety
+ && maxDeviation2 > epsilon*epsilon) {
+ // Need to take a shorter step to reliably estimate material,
+ // but only if we didn't move by safety. In order to avoid
+ // issues with roundoff we check 'step > safety' instead of
+ // 'step != safety'. If we moved by safety, there couldn't have
+ // been a boundary that we missed along the path, but we could
+ // be on the boundary.
+
+ // Take a shorter step, but never shorter than safety.
+ step = std::max(step / 2, safety);
+ } else {
+ getGeoManager()->PushPoint();
+ bool volChanged = initTrack(state7[0], state7[1], state7[2],
+ stepSign*state7[3], stepSign*state7[4],
+ stepSign*state7[5]);
+
+ if (volChanged) {
+ if (debugLvl_ > 0)
+ std::cout << " volChanged\n";
+ // Move back to start.
+ getGeoManager()->PopPoint();
+
+ // Extrapolation may not take the exact step length we asked
+ // for, so it can happen that a requested step < safety takes
+ // us across the boundary. This is then the best estimate we
+ // can get of the distance to the boundary with the stepper.
+ if (step <= safety) {
+ if (debugLvl_ > 0)
+ std::cout <<
+ " step <= safety, return stepSign*(s + step) = "
+ << stepSign*(s + step) << "\n";
+ return stepSign*(s + step);
+ }
+
+ // Volume changed during the extrapolation. Take a shorter
+ // step, but never shorter than safety.
+ step = std::max(step / 2, safety);
+ } else {
+ // we're in the new place, the step was safe, advance
+ s += step;
+
+ oldState7 = state7;
+ getGeoManager()->PopDummy(); // Pop stack, but stay in place.
+
+ findNextBoundary(fabs(sMax) - s);
+ safety = getSafeDistance();
+ slDist = getStep();
+
+ // this should not happen, but it happens sometimes.
+ // The reason are probably overlaps in the geometry.
+ // Without this "hack" many small steps with size of minStep will be made,
+ // until eventually the maximum number of iterations are exceeded and the extrapolation fails
+ if (slDist < safety)
+ slDist = safety;
+ step = slDist;
+ if (debugLvl_ > 0){
+ std::cout << " safety = " << safety
+ << "; step, slDist = " << step << "\n";
+ }
+ }
+ }
+ }
+}
+
+
+/*
+ Reference for elemental mean excitation energies at:
+http://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html
+
+Code ported from GEANT 3
+*/
+
+const int MeanExcEnergy_NELEMENTS = 93; // 0 = vacuum, 1 = hydrogen, 92 = uranium
+const double MeanExcEnergy_vals[MeanExcEnergy_NELEMENTS] = {
+ 1.e15, 19.2, 41.8, 40.0, 63.7, 76.0, 78., 82.0,
+ 95.0, 115.0, 137.0, 149.0, 156.0, 166.0, 173.0, 173.0,
+ 180.0, 174.0, 188.0, 190.0, 191.0, 216.0, 233.0, 245.0,
+ 257.0, 272.0, 286.0, 297.0, 311.0, 322.0, 330.0, 334.0,
+ 350.0, 347.0, 348.0, 343.0, 352.0, 363.0, 366.0, 379.0,
+ 393.0, 417.0, 424.0, 428.0, 441.0, 449.0, 470.0, 470.0,
+ 469.0, 488.0, 488.0, 487.0, 485.0, 491.0, 482.0, 488.0,
+ 491.0, 501.0, 523.0, 535.0, 546.0, 560.0, 574.0, 580.0,
+ 591.0, 614.0, 628.0, 650.0, 658.0, 674.0, 684.0, 694.0,
+ 705.0, 718.0, 727.0, 736.0, 746.0, 757.0, 790.0, 790.0,
+ 800.0, 810.0, 823.0, 823.0, 830.0, 825.0, 794.0, 827.0,
+ 826.0, 841.0, 847.0, 878.0, 890.0, };
+// Logarithms of the previous table, used to calculate mixtures.
+const double MeanExcEnergy_logs[MeanExcEnergy_NELEMENTS] = {
+ 34.5388, 2.95491, 3.7329, 3.68888, 4.15418, 4.33073, 4.35671, 4.40672,
+ 4.55388, 4.74493, 4.91998, 5.00395, 5.04986, 5.11199, 5.15329, 5.15329,
+ 5.19296, 5.15906, 5.23644, 5.24702, 5.25227, 5.37528, 5.45104, 5.50126,
+ 5.54908, 5.6058, 5.65599, 5.69373, 5.73979, 5.77455, 5.79909, 5.81114,
+ 5.85793, 5.84932, 5.8522, 5.83773, 5.86363, 5.8944, 5.90263, 5.93754,
+ 5.97381, 6.03309, 6.04973, 6.05912, 6.08904, 6.10702, 6.15273, 6.15273,
+ 6.1506, 6.19032, 6.19032, 6.18826, 6.18415, 6.19644, 6.17794, 6.19032,
+ 6.19644, 6.21661, 6.25958, 6.28227, 6.30262, 6.32794, 6.35263, 6.36303,
+ 6.38182, 6.41999, 6.44254, 6.47697, 6.4892, 6.51323, 6.52796, 6.54247,
+ 6.5582, 6.57647, 6.58893, 6.60123, 6.61473, 6.62936, 6.67203, 6.67203,
+ 6.68461, 6.69703, 6.71296, 6.71296, 6.72143, 6.71538, 6.67708, 6.7178,
+ 6.71659, 6.73459, 6.7417, 6.77765, 6.79122, };
+
+
+double
+MeanExcEnergy_get(int Z, bool logs = false) {
+ assert(Z >= 0 && Z < MeanExcEnergy_NELEMENTS);
+ if (logs)
+ return MeanExcEnergy_logs[Z];
+ else
+ return MeanExcEnergy_vals[Z];
+}
+
+
+double
+MeanExcEnergy_get(TGeoMaterial* mat) {
+ if (mat->IsMixture()) {
+ // The mean excitation energy is calculated as the geometric mean
+ // of the mEEs of the components, weighted for abundance.
+ double logMEE = 0.;
+ double denom = 0.;
+ TGeoMixture* mix = (TGeoMixture*)mat;
+ for (int i = 0; i < mix->GetNelements(); ++i) {
+ int index = int(mix->GetZmixt()[i]);
+ double weight = mix->GetWmixt()[i] * mix->GetZmixt()[i] / mix->GetAmixt()[i];
+ logMEE += weight * MeanExcEnergy_get(index, true);
+ denom += weight;
+ }
+ logMEE /= denom;
+ return exp(logMEE);
+ }
+
+ // not a mixture
+ int index = int(mat->GetZ());
+ return MeanExcEnergy_get(index, false);
+}
+
+double GenfitMaterialInterface::getSafeDistance()
+{
+ return getGeoManager()->GetSafeDistance();
+}
+
+double GenfitMaterialInterface::getStep()
+{
+ return getGeoManager()->GetSafeDistance();
+}
+
+TGeoNode* GenfitMaterialInterface::findNextBoundary(double stepmax,
+ const char* path,bool frombdr)
+{
+ return getGeoManager()->FindNextBoundary(stepmax,path,frombdr);
+}
+
+bool GenfitMaterialInterface::isSameLocation(double posX, double posY,
+ double posZ, bool change)
+{
+ //std::cout<<" MatInt "<<__LINE__<<" posXYZ*dd4hep::cm "<<posX*dd4hep::cm
+ //<<" posY "<<posY*dd4hep::cm<<" posZ "<<posZ*dd4hep::cm<<std::endl;
+ return getGeoManager()->IsSameLocation(posX,posY,
+ posZ,change);
+}
+
+void GenfitMaterialInterface::setCurrentDirection(double nx, double ny,
+ double nz)
+{
+ return getGeoManager()->SetCurrentDirection(nx, ny, nz);
+
+}
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h
new file mode 100644
index 00000000..83da9fc1
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h
@@ -0,0 +1,74 @@
+/////////////////////////////////////////////////////////////////////////
+// An implementation of genfit AbsMaterialInterface
+//
+// Author:
+// Yao ZHANG (zhangyao@ihep.ac.cn)
+//
+//
+/////////////////////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_GENFITMATERIALINTERFACE_H
+#define RECGENFITALG_GENFITMATERIALINTERFACE_H
+
+#include "AbsMaterialInterface.h"
+
+class RKTrackRep;
+class TGeoManager;
+class TGeoNode;
+namespace dd4hep{
+ class Detector;
+}
+
+/**
+ * @brief MaterialInterface implementation for use with ROOT's TGeoManager.
+ */
+class GenfitMaterialInterface : public genfit::AbsMaterialInterface{
+ public:
+ GenfitMaterialInterface(const dd4hep::Detector* dd4hepGeo);
+ virtual ~GenfitMaterialInterface(){};
+ static GenfitMaterialInterface* getInstance(
+ const dd4hep::Detector* dd4hepGeo);
+ void destruct();
+
+ //Set Detector pointer of dd4hep
+ void setDetector(dd4hep::Detector*);
+
+ /** @brief Initialize the navigator at given position and with given
+ direction. Returns true if the volume changed.
+ */
+ bool initTrack(double posX, double posY, double posZ,
+ double dirX, double dirY, double dirZ) override;
+
+ genfit::Material getMaterialParameters() override;
+
+
+ /** @brief Make a step (following the curvature) until step length
+ * sMax or the next boundary is reached. After making a step to a
+ * boundary, the position has to be beyond the boundary, i.e. the
+ * current material has to be that beyond the boundary. The actual
+ * step made is returned.
+ */
+ double findNextBoundary(const genfit::RKTrackRep* rep,
+ const genfit::M1x7& state7,
+ double sMax,
+ bool varField = true) override;
+
+ // ClassDef(GenfitMaterialInterface, 1);
+
+ private:
+ static GenfitMaterialInterface* m_instance;
+ TGeoManager* getGeoManager();
+ TGeoManager* m_geoManager;
+ double getSafeDistance();
+ double getStep();
+ TGeoNode* findNextBoundary(double stepmax, const char* path="",
+ bool frombdr=false);
+ bool isSameLocation(double posX, double posY, double posZ,
+ bool change=false);
+ void setCurrentDirection(double nx, double ny, double nz);
+
+};
+
+/** @} */
+
+#endif // GenfitMaterialInterface
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp b/Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp
new file mode 100644
index 00000000..6799d68e
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp
@@ -0,0 +1,19 @@
+#include "GenfitMsg.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/IMessageSvc.h"
+
+MsgStream* GenfitMsg::m_log = nullptr;
+
+MsgStream GenfitMsg::get(){
+ if(nullptr == m_log){
+ SmartIF<IMessageSvc> msgSvc(
+ Gaudi::svcLocator()->service<IMessageSvc>("MessageSvc"));
+ m_log = new MsgStream(msgSvc, "GenfitMsg");
+ (*m_log) << MSG::DEBUG << "initialize GenfitMsg" << endmsg;
+ }
+ return (*m_log);
+}
+
+GenfitMsg::~GenfitMsg(){
+ delete m_log;
+}
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMsg.h b/Reconstruction/RecGenfitAlg/src/GenfitMsg.h
new file mode 100644
index 00000000..43808b42
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMsg.h
@@ -0,0 +1,26 @@
+//////////////////////////////////////////////////////////////////////
+///
+/// This is a class for genfit interfaces classes to access Gaudi log
+/// system
+///
+/// Authors:
+/// Yao ZHANG(zhangyao@ihep.ac.cn)
+///
+//////////////////////////////////////////////////////////////////////
+#ifndef RECGENFITALG_GENFITMSG_H
+#define RECGENFITALG_GENFITMSG_H
+
+#include "GaudiKernel/MsgStream.h"
+
+class GenfitMsg {
+ public:
+ static MsgStream get();
+ private:
+ GenfitMsg(){};
+ virtual ~GenfitMsg();
+ static MsgStream* m_log;
+
+};
+
+#endif
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
new file mode 100644
index 00000000..a73f6108
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
@@ -0,0 +1,733 @@
+#include "GenfitTrack.h"
+#include "GenfitMsg.h"
+#include "GenfitField.h"
+
+//CEPCSW
+#include "DataHelper/HelixClass.h"
+
+//Externals
+#include "DD4hep/DD4hepUnits.h"
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/Track.h"
+#include "edm4hep/ReconstructedParticle.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/Vector3d.h"
+#include "DetSegmentation/GridDriftChamber.h"
+
+//genfit
+#include "Track.h"
+#include "MeasuredStateOnPlane.h"
+#include "RKTrackRep.h"
+#include "TrackPoint.h"
+#include "StateOnPlane.h"
+#include "KalmanFitterInfo.h"
+#include "KalmanFittedStateOnPlane.h"
+#include "AbsTrackRep.h"
+#include "FitStatus.h"
+#include "SpacepointMeasurement.h"
+#include "WireMeasurementNew.h"
+
+//ROOT
+#include "TRandom.h"
+#include "TVector3.h"
+#include "TLorentzVector.h"
+#include "TMatrixDSym.h"
+
+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)
+{
+ return (hit1.getTime()<hit2.getTime());
+}
+
+ 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)
+{
+
+}
+
+GenfitTrack::~GenfitTrack()
+{
+ ///Note: track reps and points will be deleted in the destructor of track
+ ///implemented in genfit::Track::Clear()
+ delete m_track;
+}
+
+/// create a Genfit track from track state, without trackRep
+/// Initialize track with seed states
+/// NO unit conversion here
+bool GenfitTrack::createGenfitTrack(int pdgType,int charge,
+ TLorentzVector posInit, TVector3 momInit, TMatrixDSym covMInit)
+{
+ TVectorD seedState(6);
+ TMatrixDSym seedCov(6);
+
+ //for(int i = 0; i < 6; ++i) {
+ // for(int j = 0; j < 6; ++j) {
+ // seedCov(i,j)=covMInit(i,j);
+ // }
+ //}
+ //yzhang FIXME
+ //seed position
+ for(int i = 0; i < 3; ++i) {
+ seedState(i)=posInit[i];
+ //yzhang TODO from covMInit to seedCov
+ double resolution = 0.1;//*dd4hep::mm/dd4hep::cm;
+ seedCov(i,i)=resolution*resolution;
+ if(i==2) seedCov(i,i)=0.5*0.5;
+ }
+ //seed momentum
+ for(int i = 3; i < 6; ++i){
+ //seedState(i)=momInit[i-3]*(dd4hep::GeV);
+ seedState(i)=momInit[i-3];
+ //yzhang TODO from covMInit to seedCov
+ seedCov(i,i)=0.01;//pow(resolution / sqrt(3),2);
+ }
+
+ if(nullptr==m_track) m_track=new genfit::Track();
+ m_track->setStateSeed(seedState);
+ m_track->setCovSeed(seedCov);
+
+ /// new a track representation and add to the track
+ int chargeId=0;
+ charge<0 ? chargeId=0 : chargeId=1;
+
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" CreateGenfitTrack seed pos("
+ <<seedState[0]<<" " <<seedState[1]<<" " <<seedState[2]<<")cm ("
+ <<seedState[3]<<" " <<seedState[4]<<" " <<seedState[5]<<")GeV pdg "
+ <<s_PDG[chargeId][pdgType]<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<"seedCov "<<endmsg;
+
+ addTrackRep(s_PDG[chargeId][pdgType]);
+
+ if(m_debug>0) seedCov.Print();
+ return true;
+}
+
+///Create a Genfit track with MCParticle, unit conversion here
+bool GenfitTrack::createGenfitTrackFromMCParticle(int pidType,
+ const edm4hep::MCParticle& mcParticle, double eventStartTime)
+{
+ ///get track parameters from McParticle
+ edm4hep::Vector3d mcPocaPos = mcParticle.getVertex();//mm
+ edm4hep::Vector3f mcPocaMom = mcParticle.getMomentum();//GeV
+ GenfitMsg::get()<<MSG::DEBUG<<"seedPos poca "<< mcPocaPos.x
+ <<" "<<mcPocaPos.y<<" "<<mcPocaPos.z<<" mm "<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<"seedMom poca "<< mcPocaMom.x
+ <<" "<<mcPocaMom.y<<" "<<mcPocaMom.z<<" GeV "<<endmsg;
+
+ ///Pivot to first layer to avoid correction of beam pipe
+ edm4hep::Vector3d firstLayerPos(1e9,1e9,1e9);
+ edm4hep::Vector3f firstLayerMom(1e9,1e9,1e9);
+ pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom);
+
+ //TODO convert unit
+ ///Get seed position and momentum
+ TLorentzVector seedPos(firstLayerPos.x,firstLayerPos.y,firstLayerPos.z,
+ eventStartTime);
+ TVector3 seedMom(firstLayerMom.x,firstLayerMom.y,firstLayerMom.z);
+ GenfitMsg::get()<<MSG::DEBUG<<"seedPos "<< firstLayerPos.x
+ <<" "<<firstLayerPos.y<<" "<<firstLayerPos.z<<endmsg;
+ GenfitMsg::get()<<MSG::DEBUG<<"seedMom "<< firstLayerMom.x
+ <<" "<<firstLayerMom.y<<" "<<firstLayerMom.z<<endmsg;
+
+ ///Get error matrix of seed track
+ TMatrixDSym covM(5);//FIXME, TODO
+
+ ///Create a genfit track with seed
+ GenfitMsg::get()<<MSG::DEBUG<<"createGenfitTrack " ;
+ if(!GenfitTrack::createGenfitTrack(pidType,mcParticle.getCharge(),
+ seedPos,seedMom,covM)){
+ GenfitMsg::get()<<MSG::ERROR<<"GenfitTrack"
+ <<" Error in createGenfitTrack" <<endmsg;
+ return false;
+ }else{
+ GenfitMsg::get()<<MSG::DEBUG<<"GenfitTrack "
+ <<"createGenfitTrackFromMCParticle track created" <<endmsg;
+ }
+ return true;
+}//end of createGenfitTrackFromMCParticle
+
+///Create a Genfit track with MCParticle, unit conversion here
+bool GenfitTrack::createGenfitTrackFromEDM4HepTrack(int pidType,
+ const edm4hep::Track& track, double eventStartTime)
+{
+ std::cout<<__FILE__<<" "<<__LINE__<<" bz kilogauss "<<m_genfitField->getBz({0.,0.,0.})/dd4hep::kilogauss<<std::endl;
+ std::cout<<__FILE__<<" "<<__LINE__<<" bz tesla "<<m_genfitField->getBz({0.,0.,0.})/dd4hep::tesla<<std::endl;
+ std::cout<<__FILE__<<" "<<__LINE__<<" bz "<<m_genfitField->getBz({0.,0.,0.})<< dd4hep::kilogauss <<" "<<dd4hep::tesla<<std::endl;
+ //TODO
+ //pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom);
+ //Get track parameters
+ edm4hep::TrackState trackStat=track.getTrackStates(0);//FIXME?
+ HelixClass helixClass;
+ helixClass.Initialize_Canonical(trackStat.phi,trackStat.D0,
+ trackStat.Z0,trackStat.omega,trackStat.tanLambda,
+ m_genfitField->getBz({0.,0.,0.})*dd4hep::kilogauss/dd4hep::tesla);
+ TLorentzVector posInit(helixClass.getReferencePoint()[0],
+ helixClass.getReferencePoint()[1],
+ helixClass.getReferencePoint()[2],eventStartTime);
+ posInit.SetX(posInit.X()*dd4hep::mm);
+ posInit.SetY(posInit.Y()*dd4hep::mm);
+ posInit.SetZ(posInit.Z()*dd4hep::mm);
+ TVector3 momInit(helixClass.getMomentum()[0],
+ helixClass.getMomentum()[1],helixClass.getMomentum()[2]);
+ momInit.SetX(momInit.x()*dd4hep::GeV);
+ momInit.SetY(momInit.y()*dd4hep::GeV);
+ momInit.SetZ(momInit.z()*dd4hep::GeV);
+ TMatrixDSym covMInit;
+ if(!createGenfitTrack(pidType,
+ int(trackStat.omega/fabs(trackStat.omega)),//charge,FIXME?
+ posInit,momInit,covMInit)){
+ return false;
+ }
+ return true;
+}
+
+/// Add a 3d SpacepointMeasurement with MC truth position smeared by sigma
+bool GenfitTrack::addSpacePointMeasurementOnTrack(const TVectorD& pos,
+ bool smear, double resolution, int detID, int hitID)
+{
+
+ /// Convert from CEPCSW unit to genfit unit, cm
+ TVectorD pos_t(3);
+ pos_t(0) = pos(0)*dd4hep::mm;
+ pos_t(1) = pos(1)*dd4hep::mm;
+ pos_t(2) = pos(2)*dd4hep::mm;
+ /// smear hit position with same weight
+ TVectorD pos_smeared(pos_t);
+ if (smear) {
+ for (int ii=0;ii<3;ii++){
+ pos_smeared[ii] += gRandom->Gaus(0, resolution/TMath::Sqrt(3.));
+ }
+ }
+
+ /// New a SpacepointMeasurement
+ TMatrixDSym hitCov(3);
+ hitCov(0,0) = resolution*resolution;
+ hitCov(1,1) = resolution*resolution;
+ hitCov(2,2) = resolution*resolution;
+
+ 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;
+
+ genfit::SpacepointMeasurement* sMeas =
+ new genfit::SpacepointMeasurement(pos_smeared,hitCov,detID,hitID,nullptr);
+ genfit::TrackPoint* trackPoint = new genfit::TrackPoint(sMeas,m_track);
+ m_track->insertPoint(trackPoint);
+
+ return true;
+}
+
+
+/// Add a WireMeasurement, no Unit conversion here
+void GenfitTrack::addWireMeasurement(double driftDistance,
+ double driftDistanceError, const TVector3& endPoint1,
+ const TVector3& endPoint2, int lrAmbig, int detID, int hitID)
+{
+ try{
+ /// New a WireMeasurement
+ genfit::WireMeasurementNew* wireMeas = new genfit::WireMeasurementNew(
+ driftDistance, driftDistanceError, endPoint1, endPoint2, detID,
+ hitID, nullptr);
+ wireMeas->setMaxDistance(0.5);//0.5 cm FIXME
+ wireMeas->setLeftRightResolution(lrAmbig);
+
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" Add wire measurement(cm) "<<hitID
+ <<" ep1("<<endPoint1[0]<<" "<<endPoint1[1]<<" "<<endPoint1[2]
+ <<") ep2("<<endPoint2[0]<<" "<<endPoint2[1]<<" "<<endPoint2[2]
+ <<") drift "<<driftDistance<<" driftErr "<<driftDistanceError
+ <<" lr "<<lrAmbig<<" detId "<<detID << " hitId "<< hitID
+ <<endmsg;
+
+ ///New a TrackPoint,create connection between meas. and trackPoint
+ genfit::TrackPoint* trackPoint=new genfit::TrackPoint(wireMeas,m_track);
+ wireMeas->setTrackPoint(trackPoint);
+
+ m_track->insertPoint(trackPoint);
+
+ }catch(genfit::Exception e){
+ GenfitMsg::get()<<MSG::DEBUG<<m_name
+ <<"Add wire measurement exception"<<endmsg;
+ e.what();
+ }
+}//end of addWireMeasurementOnTrack
+
+//Add wire measurement on wire, unit conversion here
+bool GenfitTrack::addWireMeasurementOnTrack(edm4hep::Track& track,double sigma)
+{
+ for(int iHit=0;iHit<track.trackerHits_size();iHit++){
+ edm4hep::ConstTrackerHit hit=track.getTrackerHits(iHit);
+
+ double driftVelocity=40;//FIXME, TODO, um/ns
+ double driftDistance=hit.getTime()*driftVelocity*dd4hep::um;
+ TVector3 endPointStart(0,0,0);
+ TVector3 endPointEnd(0,0,0);
+ m_gridDriftChamber->cellposition(hit.getCellID(),endPointStart,
+ endPointEnd);
+ int lrAmbig=0;
+ std::cout<<__FILE__<<" time "<<hit.getTime()<<" driftVelocity "
+ <<driftVelocity<<std::endl;
+ std::cout<<__FILE__<<" wire pos " <<endPointStart.X()<<" "
+ <<endPointStart.Y()<<" " <<endPointStart.Z()<<" "
+ <<endPointEnd.X()<<" " <<endPointEnd.Y()<<" "
+ <<endPointEnd.Z()<<" " <<std::endl;
+ endPointStart.SetX(endPointStart.x()*dd4hep::cm);
+ endPointStart.SetY(endPointStart.y()*dd4hep::cm);
+ endPointStart.SetZ(endPointStart.z()*dd4hep::cm);
+ endPointEnd.SetX(endPointEnd.x()*dd4hep::cm);
+ endPointEnd.SetY(endPointEnd.y()*dd4hep::cm);
+ endPointEnd.SetZ(endPointEnd.z()*dd4hep::cm);
+ addWireMeasurement(driftDistance,sigma*dd4hep::cm,endPointStart,
+ endPointEnd,lrAmbig,hit.getCellID(),iHit);
+ }
+ return true;
+}//end of addWireMeasurementOnTrack of Track
+
+/// Get MOP
+bool GenfitTrack::getMOP(int hitID,
+ genfit::MeasuredStateOnPlane& mop, genfit::AbsTrackRep* trackRep) const
+{
+ if(nullptr == trackRep) trackRep = getRep();
+ try{
+ mop = m_track->getFittedState(hitID,trackRep);
+ }catch(genfit::Exception e){
+ e.what();
+ return false;
+ }
+ return true;
+}
+
+/// New and add a track representation to track
+genfit::RKTrackRep* GenfitTrack::addTrackRep(int pdg)
+{
+ /// create a new track representation
+ genfit::RKTrackRep* rep = new genfit::RKTrackRep(pdg);
+ m_reps.push_back(rep);
+ m_track->addTrackRep(rep);
+ //try{
+ // genfit::MeasuredStateOnPlane stateInit(rep);
+ // rep->setPosMomCov(stateInit, pos, mom, covM);
+ //}catch(genfit::Exception e){
+ // GenfitMsg::get() << MSG::DEBUG<<m_name<<" Exception in set track status"
+ // <<endmsg ;
+ // std::cout<<e.what()<<std::endl;
+ // return false;
+ //}
+ return rep;
+}
+
+/// Get the position from genfit::Track::getStateSeed
+const TLorentzVector GenfitTrack::getSeedStatePos()const
+{
+ TVectorD seedStat(6);
+ seedStat = m_track->getStateSeed();
+ TVector3 p(seedStat[0],seedStat[1],seedStat[2]);
+ p = p*dd4hep::cm;
+ TLorentzVector pos(p[0],p[1],p[2],9999);//FIXME
+ return pos;
+}
+
+/// Get the momentum from genfit::Track::getStateSeed
+const TVector3 GenfitTrack::getSeedStateMom() const
+{
+ TVectorD seedStat(6); seedStat = m_track->getStateSeed();
+ TVector3 mom(seedStat[3],seedStat[4],seedStat[5]);
+ return mom*dd4hep::GeV;
+}
+
+/// Get the seed states of momentum and position
+void GenfitTrack::getSeedStateMom(TLorentzVector& pos, TVector3& mom) const
+{
+ TVectorD seedStat(6); seedStat = m_track->getStateSeed();
+ mom = TVector3(seedStat[3],seedStat[4],seedStat[5])*dd4hep::GeV;
+ seedStat = m_track->getStateSeed();
+ TVector3 p = TVector3(seedStat[0],seedStat[1],seedStat[2])*dd4hep::cm;
+ pos.SetXYZT(p[0],p[1],p[2],9999);//FIXME time
+}
+
+unsigned int GenfitTrack::getNumPoints() const
+{
+ return m_track->getNumPoints();
+}
+
+/// Test the fit result FIXME
+bool GenfitTrack::fitSuccess(int repID) const
+{
+
+ genfit::FitStatus* fitStatus = m_track->getFitStatus(getRep(repID));
+
+ /// Test fitting converged
+ if (!fitStatus->isFitted()||!fitStatus->isFitConverged()
+ ||fitStatus->isFitConvergedFully()) {
+ GenfitMsg::get() << MSG::DEBUG<<m_name<< "Fitting is failed... isFitted"
+ <<fitStatus->isFitted()<<" , isFitConverged "
+ <<fitStatus->isFitConverged()<<", isFitConvergedFully "
+ <<fitStatus->isFitConvergedFully()<<endmsg;
+ return false;
+ }
+
+ double chi2 = fitStatus->getChi2();
+ double ndf = fitStatus->getNdf();
+ GenfitMsg::get() << MSG::INFO << "Fit result: chi2 "<<chi2 <<" ndf "<<ndf
+ << " chi2/ndf = " << chi2/ndf<<endmsg;
+
+ /// Test fitting chi2
+ if (chi2<= 0) {
+ GenfitMsg::get() << MSG::DEBUG<<m_name<< "Fit chi2<0 (chi2,ndf) = (" <<
+ chi2 << "," << ndf << ")"<<endmsg;
+ return false;
+ }
+ return true;
+}
+
+void GenfitTrack::setDebug(int debug)
+{
+ m_debug = debug;
+ for(unsigned int i=0;i<m_reps.size();i++){ m_reps[i]->setDebugLvl(debug); }
+}
+
+void GenfitTrack::printSeed() const
+{
+ TLorentzVector pos = getSeedStatePos();
+ TVector3 mom = getSeedStateMom();
+ print(pos,mom);
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" NumPoints "<<getNumPoints()<<endmsg;
+}
+
+void GenfitTrack::printFitted(int repID) const
+{
+ TLorentzVector fittedPos;
+ TVector3 fittedMom;
+ TMatrixDSym cov;
+
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<< "printFitted nHit="
+ <<m_track->getNumPoints()<<endmsg;
+ for(unsigned int iHit=0; iHit<m_track->getNumPoints(); iHit++){
+ if (getPosMomCovMOP((int) iHit, fittedPos, fittedMom, cov, repID)){
+ //print(fittedPos,fittedMom,to_string(iHit).c_str());//TODO
+ }else{
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<"Hit "<<iHit
+ <<" have no fitter info"<<endmsg;
+ }
+ }
+}
+
+/// Print track information
+void GenfitTrack::print( TLorentzVector pos, TVector3 mom,
+ const char* comment) const
+{
+ TVector3 pglo = pos.Vect();
+ TVector3 mglo = mom;
+
+ // TODO
+ GenfitMsg::get()<<MSG::DEBUG<<m_name<<" "<<comment<<endmsg;
+
+ if(m_debug>1){
+ for(unsigned int i=0;i<m_reps.size();i++){ m_reps[i]->Print(); }
+ }
+ //for(unsigned int i=0; i<m_track->getNumPoints(); i++){
+ // m_track->getPoint(i)->print();
+ //}
+}
+
+/// Get position, momentum, cov on plane of hitID-th hit
+bool GenfitTrack::getPosMomCovMOP(int hitID, TLorentzVector& pos,
+ TVector3& mom, TMatrixDSym& cov, int repID) const
+{
+ TVector3 p;
+ genfit::MeasuredStateOnPlane mop;
+ if(!getMOP(hitID,mop,getRep(repID))) return false;
+ mop.getPosMomCov(p,mom,cov);
+ pos.SetVect(p*dd4hep::cm);
+ pos.SetT(9999);//FIXME
+ mom = mom*(dd4hep::GeV);
+ //FIXME
+ //TrackingUtils::CovConvertUnit(cov, dd4hep::cm, dd4hep::GeV);
+ return true;
+}
+
+int GenfitTrack::getNumPointsWithFittedInfo(int repID) const
+{
+ int nHitWithFittedInfo = 0;
+ int nHit = m_track->getNumPointsWithMeasurement();
+ for(int i=0; i<nHit; i++){
+ if(nullptr != m_track->getPointWithFitterInfo(i,getRep(repID))){
+ nHitWithFittedInfo++;
+ }
+ }
+ return nHitWithFittedInfo;
+}
+
+int GenfitTrack::getFittedState(TLorentzVector& pos, TVector3& mom,
+ TMatrixDSym& cov, int repID, bool biased) const
+{
+ //check number of hit with fitter info
+ if(getNumPointsWithFittedInfo(repID)<=2) return 1;
+
+ //get track rep
+ genfit::AbsTrackRep* rep = getRep(repID);
+ if(nullptr == rep) return 2;
+
+ //get first or last measured state on plane
+ genfit::MeasuredStateOnPlane mop;
+ try{
+ mop = m_track->getFittedState(biased);
+ }catch(genfit::Exception e){
+ std::cout<<" getNumPointsWithFittedInfo "
+ <<getNumPointsWithFittedInfo(repID)
+ <<" no TrackPoint with fitted info "<<std::endl;
+ GenfitMsg::get()<<MSG::DEBUG<<m_name
+ <<"Exception in getFittedState"<<endmsg;
+ std::cout<<e.what()<<std::endl;
+ return 3;
+ }
+
+ //get state
+ TVector3 p;
+ mop.getPosMomCov(p,mom,cov);
+ pos.SetVect(p*dd4hep::cm);
+ pos.SetT(9999);//FIXME
+ mom = mom*(dd4hep::GeV);
+
+ return 0;//success
+}
+
+// Get point with fitter info
+int GenfitTrack::getDetIDWithFitterInfo(int hitID, int idRaw) const
+{
+ return m_track->getPointWithFitterInfo(hitID)->
+ getRawMeasurement(idRaw)->getDetId();
+}
+
+int GenfitTrack::getPDG(int id) const
+{
+ return m_reps[id]->getPDG();
+}
+
+int GenfitTrack::getPDGCharge(int id) const
+{
+ return m_reps[id]->getPDGCharge();
+}
+
+const genfit::FitStatus*
+GenfitTrack::getFitStatus(int repID) const
+{
+ return m_track->getFitStatus(getRep(repID));
+}
+
+/// Extrapolate track to the cloest point of approach(POCA) to the wire of hit,
+/// return StateOnPlane of this POCA
+/// inputs
+/// pos,mom ... position & momentum at starting point, unit= [mm]&[GeV/c]
+/// (converted to [cm]&[GeV/c] inside this function)
+/// hit ... destination
+/// outputs poca [mm] (converted from [cm] in this function) ,global
+double GenfitTrack::extrapolateToHit( TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, TVector3 pos, TVector3 mom,
+ TVector3 end0,//one end of the hit wire
+ TVector3 end1,//the orhter end of the hit wire
+ int debug,
+ int repID,
+ bool stopAtBoundary,
+ bool calcJacobianNoise)const
+{
+
+ //genfit::MeasuredStateOnPlane state = getMOP(iPoint); // better?
+ //genfit::MeasuredStateOnPlane state = getMOP(0); // better?
+ //To do the extrapolation without IHitSelection,above 2 lines are not used.
+ pos = pos*dd4hep::cm;//FIXME
+ mom = mom*dd4hep::GeV;
+
+ //std::cout<<__LINE__<<" extrapolate to Hit pos and mom"<<std::endl;
+ //pos.Print();
+ //mom.Print();
+
+ genfit::AbsTrackRep* rep = new genfit::RKTrackRep(getRep(repID)->getPDG());
+ rep->setDebugLvl(debug);
+ genfit::MeasuredStateOnPlane state(rep);
+ rep->setPosMom(state, pos, mom);
+
+
+ //genfit::MeasuredStateOnPlane state(m_track->getRep(repID));
+ //m_track->getRep(repID)->setPosMom(state, pos, mom);
+
+ //m_track->PrintSeed();
+ double extrapoLen(0);
+ //std::cout<<" wire1 "<<std::endl;
+ //end0.Print();
+ //std::cout<<" wire0 "<<std::endl;
+ //end1.Print();
+ //std::cout<<" state "<<std::endl;
+ //state.Print();
+
+
+ // forth
+ try {
+ //genfit::RKTrackRep* rkRep =
+ // dynamic_cast<genfit::RKTrackRep*> (m_track->getRep(repID));
+ //extrapoLen = rkRep->extrapolateToLine(state, end0, end1, poca,
+ // pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise);
+ extrapoLen = rep->extrapolateToLine(state, end0, end1, poca,
+ pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise);
+ } catch (genfit::Exception& e) {
+ GenfitMsg::get() << MSG::ERROR
+ <<"Exception in GenfigFitter::ExtrapolateToHit"<<e.what()<<endmsg;
+ return extrapoLen;
+ }
+
+ //poca = poca*(dd4hep::cm);
+ //pocaOnWire = pocaOnWire*(dd4hep::cm);
+ pocaOnWire = pocaOnWire;
+ doca = (pocaOnWire-poca).Mag();
+ //TODO: debug pocaOnWire
+ GenfitMsg::get() << MSG::DEBUG << " poca "<<poca.x()<<","<<poca.y()
+ <<" "<<poca.z()<<" doca "<<doca<<endmsg;
+ GenfitMsg::get() << MSG::DEBUG << " pocaOnWire "<<pocaOnWire.x()
+ <<","<<pocaOnWire.y()<<" "<<pocaOnWire.z()<<" doca "<<doca<<endmsg;
+
+ return extrapoLen*(dd4hep::cm);
+}//end of ExtrapolateToHit
+
+
+///Add space point measurement 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);
+
+ int hitID=0;
+ for(auto simDCHit: sortedSimDCHitCol){
+ int rand=gRandom->Integer(sortedSimDCHitCol.size());
+ if(rand>200){ continue; }
+ edm4hep::Vector3d pos=simDCHit.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)){
+ hitID++;
+ }else{
+ GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurementOnTrack "
+ <<detID<<" faieled" <<endmsg;
+ }
+ }
+ GenfitMsg::get()<<MSG::DEBUG<<"GenfitTrack nHitAdded "<<hitID<<endmsg;
+ return hitID;
+}
+
+bool GenfitTrack::storeTrack(edm4hep::ReconstructedParticle& recParticle,
+ int pidType, int ndfCut, double chi2Cut)
+{
+
+ /// Get fit status
+ const genfit::FitStatus* fitState = m_track->getFitStatus();
+ double chi2 = fitState->getChi2();
+ double ndf = fitState->getNdf();
+ double charge = fitState->getCharge();
+ int isFitted = fitState->isFitted();
+ int isConverged = fitState->isFitConverged();
+ int isConvergedFully = fitState->isFitConvergedFully();
+ TMatrixDSym fittedCov;
+ TLorentzVector fittedPos;
+ TVector3 fittedMom;
+ int fittedState=getFittedState(fittedPos,fittedMom,fittedCov);
+ GenfitMsg::get()<<MSG::DEBUG<<"fit result: get status OK? pidType "
+ <<pidType<<" fittedState==0 " <<(0==fittedState)<<" isFitted "<<isFitted
+ <<" isConverged "<<isConverged<<" ndf "<<ndf<<endmsg;
+ if((0!=fittedState) || (!isFitted) || !isConvergedFully || (ndf<ndfCut)){
+ GenfitMsg::get()<<MSG::DEBUG<<" fitting failed"<<endmsg;
+ }else{
+ GenfitMsg::get()<<MSG::DEBUG<<"fit result: Pos("<<
+ fittedPos.X()<<" "<<
+ fittedPos.Y()<<" "<<
+ fittedPos.Z()<<") mom("<<
+ fittedMom.X()<<" "<<
+ fittedMom.Y()<<" "<<
+ fittedMom.Z()<<") p_tot "<<
+ fittedMom.Mag()<<" pt "<<
+ fittedMom.Perp()<<
+ " chi2 "<<chi2<<
+ " ndf "<<ndf
+ <<endmsg;
+ }
+ float pos[3]={fittedPos.X(),fittedPos.Y(),fittedPos.Z()};
+ float mom[3]={fittedMom.X(),fittedMom.Y(),fittedMom.Z()};
+ HelixClass helix;
+ helix.Initialize_VP(pos,mom,charge,m_genfitField->getBz(fittedPos.Vect()));
+
+
+ /////track status at POCA to origin
+ //TVector3 origin(pivot);
+ //TVector3 pocaToOrigin_pos(1e9*dd4hep::cm,1e9*dd4hep::cm,1e9*dd4hep::cm);
+ //TVector3 pocaToOrigin_mom(1e9*dd4hep::GeV,1e9*dd4hep::GeV,1e9*dd4hep::GeV);
+
+ //if(ExtrapolateToPoint(pocaToOrigin_pos,pocaToOrigin_mom,
+ // m_track,origin) > 1e6*dd4hep::cm){
+ // log<<"extrapolate to origin failed";
+ // return false;
+ //}
+
+
+ //new TrackState
+ edm4hep::TrackState* trackState = new edm4hep::TrackState();
+ trackState->location=0;//edm4hep::AtIP;
+ trackState->D0=helix.getD0();
+ trackState->phi=helix.getPhi0();
+ trackState->omega=helix.getOmega();
+ trackState->Z0=helix.getZ0();
+ trackState->tanLambda=helix.getTanLambda();
+ trackState->referencePoint=helix.getReferencePoint();
+
+ // std::array<float,15> covMatrix;
+ // int k=0;
+ // for(int i=0;i<5;i++){
+ // for(int j=0;j<5;j++){
+ // if(i<=j) covMatrix[k]=;//FIXME
+ // }
+ // }
+ // trackState.covMatrix=
+
+ //new Track
+ edm4hep::Track* track = new edm4hep::Track();
+ //track->setType();
+ track->setChi2(fitState->getChi2());
+ track->setNdf(fitState->getNdf());
+ //track->setDEdx();
+ //track->setRadiusOfInnermostHit();//FIXME
+ //track->addToTrackerHits();
+
+ //new ReconstructedParticle
+ //recParticle->setType();
+ //dcRecParticle->setEnergy();
+
+ edm4hep::Vector3f momVec3(helix.getMomentum()[0],
+ helix.getMomentum()[1],helix.getMomentum()[2]);
+ recParticle.setMomentum(momVec3);
+ //recParticle->setReferencePoint(referencePoint);
+ recParticle.setCharge(helix.getCharge());
+ // recParticle->setMass();
+ // recParticle->setCovMatrix();
+ // rcRecParticle->setStartVertex();
+ //recParticle->addToTracks(track);
+
+ return true;
+}
+
+void GenfitTrack::pivotToFirstLayer(edm4hep::Vector3d& pos,
+ edm4hep::Vector3f& mom, edm4hep::Vector3d& firstPos,
+ edm4hep::Vector3f& firstMom)
+{
+ //FIXME, TODO
+ firstPos=pos;
+ firstMom=mom;
+}
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
new file mode 100644
index 00000000..6e2823e4
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
@@ -0,0 +1,199 @@
+//////////////////////////////////////////////////////////////////
+///
+/// This is an interface of to call genfit
+/// A genfit track can be created, fitted and extrapolated
+/// Track with only track representation(s)
+///
+/// In this file, including:
+/// a genfit track class
+///
+/// Units are following DD4hepUnits
+///
+/// Authors:
+/// Zhang Yao (zhangyao@ihep.ac.cn)
+/// Y.Fujii (yfujii@ihep.ac.cn)
+/// Yohei Nakatsugawa (yohei@ihep.ac.cn)
+///
+//////////////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_GENFITTRACK_H
+#define RECGENFITALG_GENFITTRACK_H
+
+#include "GenfitFitter.h"
+
+//ROOT
+#include "TVectorD.h"
+#include "TMatrixDSym.h"
+
+//STL
+#include <vector>
+
+class TLorentzVector;
+
+namespace genfit{
+ class Track;
+ class FitStatus;
+ class AbsTrackRep;
+ class RKTrackRep;
+ class KalmanFittedStateOnPlane;
+}
+namespace edm4hep{
+ class MCParticle;
+ class SimTrackerHitCollection;
+ class ReconstructedParticle;
+ class Track;
+ class Vector3d;
+ class Vector3f;
+}
+namespace dd4hep {
+ namespace DDSegmentation{
+ class GridDriftChamber;
+ }
+}
+
+class GenfitTrack {
+ friend int GenfitFitter::processTrack(
+ GenfitTrack* track, bool resort=false);
+
+ friend int GenfitFitter::processTrackWithRep(
+ GenfitTrack* track, int repID=0, bool resort=false);
+
+ friend double GenfitFitter::extrapolateToHit(TVector3& poca,
+ TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, const GenfitTrack* track,
+ TVector3 pos, TVector3 mom, TVector3 end0, TVector3 end1, int debug,
+ int repID=0, bool stopAtBoundary=false, bool calcJacobianNoise=false);
+
+ friend double GenfitFitter::extrapolateToCylinder(TVector3& pos,
+ TVector3& mom,
+ GenfitTrack* track, double radius, const TVector3 linePoint,
+ const TVector3 lineDirection, int hitID =0, int repID=0,
+ bool stopAtBoundary=false, bool calcJacobianNoise=false);
+
+ friend double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom,
+ const GenfitTrack* genfitTrack, const TVector3& point, int repID=0,
+ bool stopAtBoundary = false, bool calcJacobianNoise = false) const;
+
+ public:
+ GenfitTrack(const GenfitField* field,
+ const dd4hep::DDSegmentation::GridDriftChamber* seg,
+ const char* name="GenfitTrack");
+ virtual ~GenfitTrack();
+
+ /// Add a Genfit track
+ virtual bool createGenfitTrack(int pdgType,int charge,TLorentzVector pos, TVector3 mom,
+ TMatrixDSym covM);
+ //virtual bool createGenfitTrack(TLorentzVector posInit,TVector3 momInit,
+ //TMatrixDSym covMInit);
+
+ ///Create genfit track from MCParticle
+ bool createGenfitTrackFromMCParticle(int pidTyep,const edm4hep::MCParticle&
+ mcParticle, double eventStartTime=0.);
+ bool createGenfitTrackFromEDM4HepTrack(int pidType,const edm4hep::Track& track,
+ double eventStartTime);
+
+ // /// Prepare a hit list, return number of hits on track
+ // int PrepareHits();//TODO
+
+ /// Add a space point measurement, return number of hits on track
+ virtual bool addSpacePointMeasurementOnTrack(const TVectorD&, bool, double,
+ int detID=-1, int hitID=-1);
+
+ /// Add a WireMeasurement with MC truth position smeared by sigma
+ virtual void addWireMeasurement(double driftDistance,
+ double driftDistanceError, const TVector3& endPoint1,
+ const TVector3& endPoint2, int lrAmbig, int detID, int hitID);
+
+ /// Add a WireMeasurement with DC digi
+ virtual bool addWireMeasurementOnTrack(edm4hep::Track& track, double sigma);
+
+ ///Add space point from truth to track
+ int addSpacePointMeasurementOnTrack(const edm4hep::MCParticle& mcParticle,
+ const edm4hep::SimTrackerHitCollection*& simDCHitCol,bool smear=true,
+ double resolution=0.01);
+
+ ///Store track to ReconstructedParticle
+ bool storeTrack(edm4hep::ReconstructedParticle& dcRecParticle,int pidType,
+ int ndfCut, double chi2Cut);
+
+ ///A tool to convert track to the first layer of DC
+ void pivotToFirstLayer(edm4hep::Vector3d& pos,edm4hep::Vector3f& mom,
+ edm4hep::Vector3d& firstPos, edm4hep::Vector3f& firstMom);
+
+ /// Copy a track to event
+ //void CopyATrack()const;
+
+ ///Extrapolate to Hit
+ double extrapolateToHit( TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca, TVector3 pos, TVector3 mom,
+ TVector3 end0,//one end of the hit wire
+ TVector3 end1,//the orhter end of the hit wire
+ int debug,
+ int repID,
+ bool stopAtBoundary,
+ bool calcJacobianNoise) const;
+
+ bool getPosMomCovMOP(int hitID, TLorentzVector& pos, TVector3& mom,
+ TMatrixDSym& cov, int repID=0) const;
+
+ /// get the seed position and momentum from track
+ const TLorentzVector getSeedStatePos() const;
+ const TVector3 getSeedStateMom() const;
+ void getSeedStateMom(TLorentzVector& pos, TVector3& mom) const;
+ unsigned int getNumPoints() const;
+
+ /// get the fitted track status
+ const genfit::FitStatus* getFitStatus(int repID=0) const;
+ int getFittedState(TLorentzVector& pos, TVector3& mom, TMatrixDSym& cov,
+ int repID=0, bool biased=true) const;
+ int getNumPointsWithFittedInfo(int repID=0) const;
+ bool getFirstPointWithFittedInfo(int repID=0) const;
+ bool fitSuccess(int repID) const;
+
+ /// get the wire infomation
+ int getDetIDWithFitterInfo(int hitID, int idRaw=0) const;
+
+ int getPDG(int id=0) const;
+ int getPDGCharge(int id=0) const;
+
+ /// print genfit track
+ void printSeed() const;//print seed
+ void printFitted(int repID=0) const;//print seed
+ void print(TLorentzVector pos, TVector3 mom, const char* comment="") const;
+
+ /// set and get debug level
+ void setDebug(int debug);
+ int getDebug(void) const { return m_debug;}
+
+ /// get name of this object
+ const char* getName() const {return m_name;}
+ genfit::Track* getTrack() const{return m_track;}
+
+ /// Add a track representation
+ genfit::RKTrackRep* addTrackRep(int pdg);
+
+ protected:
+ //genfit::Track* getTrack() {return m_track;}
+
+ private:
+
+ const char* m_name;
+
+ ///Note! private functions are using genfit unit, cm and MeV
+
+ genfit::AbsTrackRep* getRep(int id=0) const {return m_reps[id];}
+ bool getMOP(int hitID, genfit::MeasuredStateOnPlane& mop,
+ genfit::AbsTrackRep* trackRep=nullptr) const;
+
+ genfit::Track* m_track;/// track
+ std::vector<genfit::AbsTrackRep*> m_reps;/// track repesentations
+ int m_debug;/// debug level
+
+ const GenfitField* m_genfitField;//pointer to genfit field
+ const dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber;
+
+ static const int s_PDG[2][5];
+
+};
+
+#endif
diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
new file mode 100644
index 00000000..99d7d3a5
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
@@ -0,0 +1,416 @@
+#include "RecGenfitAlgDC.h"
+#include "GenfitTrack.h"
+#include "GenfitFitter.h"
+#include "GenfitField.h"
+
+//genfit
+#include "EventDisplay.h"
+
+//cepcsw
+#include "DetInterface/IGeomSvc.h"
+#include "DataHelper/HelixClass.h"
+#include "DetSegmentation/GridDriftChamber.h"
+
+//externals
+#include "edm4hep/EventHeaderCollection.h"
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/ReconstructedParticle.h"
+#include "edm4hep/ReconstructedParticleCollection.h"
+#include "edm4hep/Track.h"
+#include "edm4hep/TrackCollection.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/Detector.h"
+#include "UTIL/BitField64.h"
+#include "DDSegmentation/Segmentation.h"
+#include "TRandom.h"
+#include "TLorentzVector.h"
+
+//stl
+#include "time.h"
+
+
+DECLARE_COMPONENT( RecGenfitAlgDC )
+
+///////////////////////////////////////////////////////////////////////////
+
+RecGenfitAlgDC::RecGenfitAlgDC(const std::string& name, ISvcLocator* pSvcLocator) :
+ GaudiAlgorithm(name, pSvcLocator),m_nPDG(5),m_dd4hep(nullptr),m_gridDriftChamber(nullptr),m_decoder(nullptr)
+{
+ declareProperty("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");
+ declareProperty("DriftChamberHitsCollection", m_simDCHitCol,
+ "Handle of the input SimHit collection");
+ declareProperty("DigiDCHitCollection", m_digiDCHitsCol,
+ "Handle of digi DCHit collection");
+ declareProperty("DCTrackCollection", m_dcTrackCol,
+ "Handle of DC track collection");
+ //declareProperty("DCHitAssociationCollection", _dcHitAssociationCol,
+ // "Handle of association collection");
+ declareProperty("DCRecParticleCollection", m_dcRecParticleCol,
+ "Handle of drift chamber reconstructed particle collection");
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode RecGenfitAlgDC::initialize()
+{
+ MsgStream log(msgSvc(), name());
+ info()<<" RecGenfitAlgDC initialize()"<<endmsg;
+
+ ///Get GeomSvc
+ m_geomSvc=Gaudi::svcLocator()->service("GeomSvc");
+ if (nullptr==m_geomSvc) {
+ std::cout<<"Failed to find GeomSvc"<<std::endl;
+ return StatusCode::FAILURE;
+ }
+ ///Get Detector
+ m_dd4hep = m_geomSvc->lcdd();
+ ///Get Field
+ m_dd4hepField=m_geomSvc->lcdd()->field();
+
+ /// New a genfit fitter
+ m_genfitFitter=new GenfitFitter(m_fitterType.c_str());
+ m_genfitField=new GenfitField(m_dd4hepField);
+ m_genfitFitter->setField(m_genfitField);
+ m_genfitFitter->setGeoMaterial(m_geomSvc->lcdd());
+ m_genfitFitter->setEnergyLossBrems(m_correctBremsstrahlung);
+ m_genfitFitter->setNoiseBrems(m_correctBremsstrahlung);
+ if(m_debug>10) m_genfitFitter->setDebug(m_debug-10);
+ if(m_noMaterialEffects) m_genfitFitter->setNoEffects(true);
+ if(-1==m_debugPid) m_genfitFitter->setNoEffects(true);
+ if(-1==m_debugPid) m_debugPid=0;//charged geantino with electron pid
+ if(m_fitterType=="DAF"||m_fitterType=="DafRef"){
+ m_genfitFitter->setMaxIterationsBetas(m_bStart,m_bFinal,m_maxIteration);
+ } else {
+ m_genfitFitter->setMaxIterations(m_maxIteration);
+ }
+ //print genfit parameters
+ if(m_debug) m_genfitFitter->print();
+ if(""!=m_genfitHistRootName) m_genfitFitter->initHist(m_genfitHistRootName);
+
+ //initialize member vairables
+ for(int i=0;i<5;i++) m_fitSuccess[i]=0;
+ m_nDCTrack=0;
+ ///Get Readout
+ dd4hep::Readout readout=m_dd4hep->readout(m_readout_name);
+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");
+ if(nt){
+ m_tuple=nt;
+ }else{
+ m_tuple=ntupleSvc()->book("RecGenfitAlgDC/genfit",
+ CLID_ColumnWiseTuple,"RecGenfitAlgDC");
+ if(m_tuple){
+ StatusCode sc;
+ sc=m_tuple->addItem("run",m_run);
+ sc=m_tuple->addItem("evt",m_evt);
+ sc=m_tuple->addItem("tkId",m_tkId);
+ sc=m_tuple->addItem("mcIndex",m_mcIndex,0,100);//max. 100 particles
+ sc=m_tuple->addItem("truthPocaMc",m_mcIndex,m_pocaPosMc,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("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("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;
+ }else{
+ log << MSG::ERROR << "Cannot book tuple RecGenfitAlgDC/genfit" << endmsg;
+ }
+ }//end of book tuple
+
+ //init genfit event display
+ //if(m_showDisplay) m_genfitDisplay = genfit::EventDisplay::getInstance();
+
+
+ return StatusCode::SUCCESS;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode RecGenfitAlgDC::execute()
+{
+
+ m_timer=clock();
+ m_firstPid=true;
+ info()<<" RecGenfitAlgDC in execute()"<<endmsg;
+ StatusCode sc;
+
+ /////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;
+
+ 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;
+ return StatusCode::SUCCESS;
+ }
+ const edm4hep::TrackerHitCollection* didiDCHitsCol=nullptr;
+ didiDCHitsCol=m_digiDCHitsCol.get();
+ if(nullptr==didiDCHitsCol) {
+ log<<MSG::DEBUG<<"DigiDCHitsCollection not found"<<endmsg;
+ return StatusCode::SUCCESS;
+ }
+ edm4hep::ReconstructedParticleCollection* dcRecParticleCol=
+ m_dcRecParticleCol.createAndPut();
+ ///----------------------------------------------------
+ ///Loop over Track and do fitting for each track
+ ///----------------------------------------------------
+ log<<MSG::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++){
+ if((m_debugPid>=0) && (m_debugPid!=pidType)) continue;
+ ///-----------------------------------
+ ///Create a GenFit track
+ ///-----------------------------------
+ GenfitTrack* genfitTrack=new GenfitTrack(m_genfitField,
+ m_gridDriftChamber);
+ genfitTrack->setDebug(m_debug);
+ double eventStartTime=0;
+ if(!genfitTrack->createGenfitTrackFromEDM4HepTrack(pidType,dcTrack,
+ eventStartTime)){
+ log<<MSG::DEBUG<<"createGenfitTrack 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_debug) genfitTrack->printSeed();
+
+ ///-----------------------------------
+ ///call genfit fitting procedure
+ ///-----------------------------------
+ m_genfitFitter->processTrack(genfitTrack,m_resortHits);
+ m_genfitFitter->setDebug(m_debug);
+
+ ///-----------------------------------
+ ///Store track
+ ///-----------------------------------
+ auto dcRecParticle=dcRecParticleCol->create();
+ genfitTrack->storeTrack(dcRecParticle,pidType,m_ndfCut,
+ m_chi2Cut);
+
+ if(m_tuple) debugTrack(pidType,genfitTrack);
+ if(m_showDisplay) {
+ //m_genfitDisplay->addEvent(genfitTrack->getTrack());
+ //m_genfitDisplay->open();
+ }else{
+ delete genfitTrack;
+ }
+ }//end loop over particle type
+ }//end loop over a track
+
+ if(m_tuple) debugEvent();
+
+ return StatusCode::SUCCESS;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode RecGenfitAlgDC::finalize()
+{
+ MsgStream log(msgSvc(), name());
+ log << MSG::INFO << " RecGenfitAlgDC in finalize()" << endmsg;
+
+ m_genfitFitter->writeHist();
+ delete m_genfitFitter;
+ std::cout<<"RecGenfitAlgDC nRecTrack="<<m_nDCTrack<<" success e "
+ <<m_fitSuccess[0]<<" mu "<<m_fitSuccess[1]<<" pi "<<m_fitSuccess[2]
+ <<" K "<<m_fitSuccess[3]<<" p "<<m_fitSuccess[4]<<std::endl;
+ if(m_nDCTrack>0){
+ std::cout<<"RecGenfitAlgDC Success rate = "<<std::endl;
+ for (int i=0;i<5;i++){
+ std::cout<<Form("%d %2.2f",i,((float) m_fitSuccess[i])/m_nDCTrack)
+ <<std::endl;
+ }
+ }
+ return StatusCode::SUCCESS;
+}
+
+void RecGenfitAlgDC::debugTrack(int pidType,const GenfitTrack* genfitTrack)
+{
+ /// Get fit status
+ const genfit::FitStatus* fitState = genfitTrack->getFitStatus();
+ double chi2 = fitState->getChi2();
+ double ndf = fitState->getNdf();
+ double charge = fitState->getCharge();
+ int isFitted = fitState->isFitted();
+ int isConverged = fitState->isFitConverged();
+ int isConvergedFully = 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]={fittedPos.X()/dd4hep::mm,fittedPos.Y()/dd4hep::mm,fittedPos.Z()/dd4hep::mm};
+ float mom[3]={fittedMom.X(),fittedMom.Y(),fittedMom.Z()};
+ helix.Initialize_VP(pos,mom,charge,m_genfitField->getBz(fittedPos.Vect()));
+ m_pocaMomKalP[pidType]=fittedMom.Mag();
+
+}
+
+void RecGenfitAlgDC::debugEvent()
+{
+ const edm4hep::MCParticleCollection* mcParticleCol = nullptr;
+ const edm4hep::SimTrackerHitCollection* simDCHitCol=nullptr;
+
+ mcParticleCol=m_mcParticleCol.get();
+ simDCHitCol=m_simDCHitCol.get();
+ m_nHitMdc=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;
+ std::cout<<__FILE__<<" "<<px<<" "<<py<<" "<<pz<<std::endl;
+ m_pocaMomMcP[iMcParticle]=sqrt(px*px+py*py+pz*pz);
+ iMcParticle++;
+ }
+ m_mcIndex=iHit;
+
+ if(m_tuple) m_tuple->write();
+}
diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h
new file mode 100644
index 00000000..8b2754f7
--- /dev/null
+++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h
@@ -0,0 +1,179 @@
+//////////////////////////////////////////////////////////////////////
+///
+/// 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
+///
+/// Units are following DD4hepUnits
+///
+/// Authors:
+/// Yao ZHANG(zhangyao@ihep.ac.cn)
+///
+/////////////////////////////////////////////////////////////////////
+
+#ifndef RECGENFITALG_RECGENFITALGDC_H
+#define RECGENFITALG_RECGENFITALGDC_H
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/NTuple.h"
+#include "k4FWCore/DataHandle.h"
+#include "DD4hep/Fields.h"
+#include <string>
+
+class GenfitFitter;
+class GenfitField;
+class GenfitTrack;
+class IGeomSvc;
+class time;
+namespace genfit{
+ class EventDisplay;
+}
+namespace dd4hep {
+ class Detector;
+ //class rec::CellIDPositionConverter;
+ namespace DDSegmentation{
+ class GridDriftChamber;
+ class BitFieldCoder;
+ }
+}
+namespace edm4hep{
+ class EventHeaderCollection;
+ class MCParticleCollection;
+ class SimTrackerHitCollection;
+ class TrackCollection;
+ class TrackerHitCollection;
+ class MCRecoTrackerAssociationCollection;
+ class ReconstructedParticle;
+ class ReconstructedParticleCollection;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+class RecGenfitAlgDC:public GaudiAlgorithm {
+ public:
+ RecGenfitAlgDC (const std::string& name, ISvcLocator* pSvcLocator);
+ StatusCode initialize() override;
+ StatusCode execute() override;
+ StatusCode finalize() override;
+
+ private:
+ GenfitFitter* m_genfitFitter;//The pointer to a GenfitFitter
+ const GenfitField* m_genfitField;//The pointer to a GenfitField
+
+ 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};
+ 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};
+ //output hits and particles
+ DataHandle<edm4hep::TrackerHitCollection> m_dcFitRecHitCol{
+ "DCFitRecHitsCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::ReconstructedParticleCollection> m_dcRecParticleCol{
+ "DCRecParticleCollection", Gaudi::DataHandle::Writer, this};
+
+ const unsigned int m_nPDG;//5:e,mu,pi,K,proton
+ SmartIF<IGeomSvc> m_geomSvc;
+ dd4hep::OverlayedField m_dd4hepField;
+ dd4hep::Detector* m_dd4hep;
+ dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+ Gaudi::Property<std::string> m_readout_name{this, "readout",
+ "DriftChamberHitsCollection"};
+ 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
+ int m_fitSuccess[5];
+ int m_nDCTrack;
+ int m_debugPid;
+ bool m_useTruthSeed;
+ 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;
+
+ /// 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_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_nDofKal;
+ NTuple::Array<double> m_chi2Kal;
+ NTuple::Array<bool> m_convergeKal;
+ NTuple::Array<bool> m_isFittedKal;
+ NTuple::Item<int> m_nDigi;
+ NTuple::Item<int> m_nHitMc;
+ NTuple::Item<int> m_nHitMdc;
+ NTuple::Item<int> m_nClusterCgem;
+ NTuple::Item<int> m_nHitKalInput;
+ NTuple::Array<double> m_mdcHitDriftT;
+ NTuple::Array<double> m_mdcHitDriftDl;
+ NTuple::Array<double> m_mdcHitDriftDr;
+ NTuple::Array<int> m_mdcHitLr;
+ NTuple::Array<int> m_mdcHitLayer;
+ NTuple::Array<int> m_mdcHitWire;
+ NTuple::Array<double> m_mdcHitExpDoca;
+ NTuple::Array<double> m_mdcHitExpMcDoca;
+ NTuple::Array<double> m_mdcHitErr;
+ NTuple::Array<int> m_nHitFailedKal;
+ NTuple::Array<int> m_nHitFitted;
+ NTuple::Array<double> m_time;
+ //truth
+ NTuple::Array<int> m_mdcHitMcLr;
+ NTuple::Array<int> m_mdcHitMcTkId;
+ NTuple::Array<double> m_mdcHitMcDrift;
+ NTuple::Array<double> m_mdcHitMcX;
+ NTuple::Array<double> m_mdcHitMcY;
+ NTuple::Array<double> m_mdcHitMcZ;
+ NTuple::Array<double> m_mdcHitExpMcPocaX;
+ NTuple::Array<double> m_mdcHitExpMcPocaY;
+ NTuple::Array<double> m_mdcHitExpMcPocaZ;
+ NTuple::Array<double> m_mdcHitExpMcPocaWireX;
+ NTuple::Array<double> m_mdcHitExpMcPocaWireY;
+ NTuple::Array<double> m_mdcHitExpMcPocaWireZ;
+
+};
+#endif
diff --git a/Reconstruction/Tracking/CMakeLists.txt b/Reconstruction/Tracking/CMakeLists.txt
index 44b8702a..0a4452cc 100644
--- a/Reconstruction/Tracking/CMakeLists.txt
+++ b/Reconstruction/Tracking/CMakeLists.txt
@@ -1,9 +1,9 @@
gaudi_subdir(Tracking v0r0)
find_package(GEAR REQUIRED)
-find_package(GSL REQUIRED )
-find_package(LCIO REQUIRED )
-find_package(EDM4HEP REQUIRED )
+find_package(GSL REQUIRED )
+find_package(LCIO REQUIRED )
+find_package(EDM4HEP REQUIRED )
find_package(DD4hep COMPONENTS DDCore DDRec REQUIRED)
@@ -11,17 +11,21 @@ gaudi_depends_on_subdirs(
Service/GearSvc
Service/EventSeeder
Service/TrackSystemSvc
+ Detector/DetSegmentation
)
set(Tracking_srcs
src/Clupatra/*.cpp
src/FullLDCTracking/*.cpp
+ src/TruthTracker/*.cpp
)
# Modules
gaudi_add_module(Tracking ${Tracking_srcs}
INCLUDE_DIRS GaudiKernel gear ${GSL_INCLUDE_DIRS} ${LCIO_INCLUDE_DIRS}
- LINK_LIBRARIES GaudiAlgLib GaudiKernel ${GEAR_LIBRARIES} ${GSL_LIBRARIES} ${LCIO_LIBRARIES} TrackSystemSvcLib
+ LINK_LIBRARIES GaudiAlgLib GaudiKernel ${GEAR_LIBRARIES} ${GSL_LIBRARIES}
+ ${LCIO_LIBRARIES} TrackSystemSvcLib
+ DetSegmentation
-Wl,--no-as-needed
EDM4HEP::edm4hep EDM4HEP::edm4hepDict
-Wl,--as-needed
diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
new file mode 100644
index 00000000..b6da32ed
--- /dev/null
+++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
@@ -0,0 +1,203 @@
+#include "TruthTrackerAlg.h"
+#include "DataHelper/HelixClass.h"
+#include "edm4hep/EventHeaderCollection.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/TrackCollection.h"
+#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)
+{
+ 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,
+ "Handle of drift chamber reconstructed particle collection");
+ declareProperty("DCRecParticleAssociationCollection", m_dcRecParticleAssociationCol,
+ "Handle of drift chamber reconstructed particle collection");
+ declareProperty("debug", m_debug=false);
+}
+
+StatusCode TruthTrackerAlg::initialize()
+{
+ ///Get geometry
+ m_geomSvc=service<IGeomSvc>("GeomSvc");
+ if (!m_geomSvc) {
+ error() << "Failed to get GeomSvc." << endmsg;
+ return StatusCode::FAILURE;
+ }
+ ///Get Detector
+ m_dd4hep=m_geomSvc->lcdd();
+ if (nullptr==m_dd4hep) {
+ error() << "Failed to get dd4hep::Detector." << endmsg;
+ return StatusCode::FAILURE;
+ }
+ //Get Field
+ m_dd4hepField=m_geomSvc->lcdd()->field();
+ ///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 N-tuple 1
+ NTuplePtr nt1(ntupleSvc(), "TruthTrackerAlg/1");
+ if ( nt1 ) {
+ m_tuple_id = nt1;
+ } else {
+ m_tuple_id = ntupleSvc()->book( "TruthTrackerAlg/1", CLID_RowWiseTuple, "Row-wise N-Tuple example" );
+ if ( m_tuple_id ) {
+ m_tuple_id->addItem( "system", m_id_system ).ignore();
+ m_tuple_id->addItem( "module", m_id_module ).ignore();
+ m_tuple_id->addItem( "stave", m_id_stave ).ignore();
+ m_tuple_id->addItem( "tower", m_id_tower ).ignore();
+ m_tuple_id->addItem( "layer", m_id_layer ).ignore();
+ m_tuple_id->addItem( "wafer", m_id_wafer ).ignore();
+ m_tuple_id->addItem( "cellX", m_id_cellX ).ignore();
+ m_tuple_id->addItem( "cellY", m_id_cellY ).ignore();
+
+ } else { // did not manage to book the N tuple....
+ error() << " Cannot book N-tuple:" << long( m_tuple_id ) << endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+
+
+ return GaudiAlgorithm::initialize();
+}
+
+StatusCode TruthTrackerAlg::execute()
+{
+ ///Retrieve MC particle(s)
+ const edm4hep::MCParticleCollection* mcParticleCol=nullptr;
+ mcParticleCol=m_mcParticleCol.get();
+ if(nullptr==mcParticleCol){
+ debug()<<"MCParticleCollection not found"<<endmsg;
+ return StatusCode::SUCCESS;
+ }
+ ///Retrieve DC digi
+ const edm4hep::TrackerHitCollection* digiDCHitsCol=nullptr;
+ digiDCHitsCol=m_digiDCHitsCol.get();//FIXME DEBUG
+ if(nullptr==digiDCHitsCol){
+ debug()<<"TrackerHitCollection not found"<<endmsg;
+ //return StatusCode::SUCCESS;
+ }
+
+ ///Output Track collection
+ edm4hep::TrackCollection* dcTrackCol=m_dcTrackCol.createAndPut();
+ ////TODO
+ /////Output MCRecoTrackerAssociationCollection collection
+ //const edm4hep::MCRecoTrackerAssociationCollection*
+ // mcRecoTrackerAssociationCol=nullptr;
+ //if(nullptr==mcRecoTrackerAssociationCol){
+ // log<<MSG::DEBUG<<"MCRecoTrackerAssociationCollection not found"
+ // <<endmsg;
+ // return StatusCode::SUCCESS;
+ //}
+ //mcRecoTrackerAssociationCol=m_mcRecoParticleAssociation.get();
+
+ ///Convert MCParticle to Track and ReconstructedParticle
+ if(m_debug){
+ 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);
+ 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);
+ edm4hep::TrackState trackState;
+ trackState.D0=helix.getD0();
+ trackState.phi=helix.getPhi0();
+ trackState.omega=helix.getOmega();
+ trackState.Z0=helix.getZ0();
+ trackState.tanLambda=helix.getTanLambda();
+ trackState.referencePoint=helix.getReferencePoint();
+ std::array<float,15> covMatrix;
+ 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);
+ }
+ }
+ //TODO
+ //new ReconstructedParticle
+ //recParticle->setType();
+ //dcRecParticle->setEnergy();
+ //edm4hep::Vector3f momVec3(helix.getMomentum()[0],
+ // helix.getMomentum()[1],helix.getMomentum()[2]);
+ //recParticle->setMomentum(momVec3);
+ }
+
+ if(m_debug) std::cout<<"Output DCTrack size="<<dcTrackCol->size()<<std::endl;
+ return StatusCode::SUCCESS;
+}
+
+StatusCode TruthTrackerAlg::finalize()
+{
+ return GaudiAlgorithm::finalize();
+}
diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
new file mode 100644
index 00000000..f8f31c13
--- /dev/null
+++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
@@ -0,0 +1,79 @@
+#ifndef TruthTrackerAlg_h
+#define TruthTrackerAlg_h
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/NTuple.h"
+#include "DD4hep/Fields.h"
+
+class IGeomSvc;
+namespace dd4hep {
+ class Detector;
+ //class rec::CellIDPositionConverter;
+ namespace DDSegmentation{
+ class GridDriftChamber;
+ class BitFieldCoder;
+ }
+}
+namespace edm4hep {
+ class MCParticleCollection;
+ class TrackerHitCollection;
+ class TrackCollection;
+ class MCRecoTrackerAssociationCollection;
+ class ReconstructedParticleCollection;
+ class MCRecoParticleAssociationCollection;
+}
+
+class TruthTrackerAlg: public GaudiAlgorithm
+{
+ public:
+ TruthTrackerAlg(const std::string& name, ISvcLocator* svcLoc);
+
+ virtual StatusCode initialize();
+ virtual StatusCode execute();
+ virtual StatusCode finalize();
+
+ private:
+ SmartIF<IGeomSvc> m_geomSvc;
+ dd4hep::Detector* m_dd4hep;
+ dd4hep::OverlayedField m_dd4hepField;
+ //dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
+ dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+
+ //reader
+ DataHandle<edm4hep::MCParticleCollection> m_mcParticleCol{
+ "MCParticle", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_digiDCHitsCol{
+ "DigiDCHitsCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection>
+ m_dcHitAssociationCol{ "DCHitAssociationCollection",
+ Gaudi::DataHandle::Reader, this};
+ //writer
+ DataHandle<edm4hep::TrackCollection> m_dcTrackCol{"DCTrackCollection",
+ Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::ReconstructedParticleCollection> m_dcRecParticleCol{
+ "DCRecParticleCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoParticleAssociationCollection>
+ m_dcRecParticleAssociationCol{"DCRecMCRecoParticleAssociationCollection",
+ Gaudi::DataHandle::Writer, this};
+
+ //readout for getting segmentation
+ Gaudi::Property<std::string> m_readout_name{this, "readout",
+ "DriftChamberHitsCollection"};
+
+ int m_debug;
+ //strore all the id for later analysis
+ NTuple::Tuple* m_tuple_id = nullptr;
+
+ NTuple::Item<int> m_id_system;
+ NTuple::Item<int> m_id_module;
+ NTuple::Item<int> m_id_stave;
+ NTuple::Item<int> m_id_tower;
+ NTuple::Item<int> m_id_layer;
+ NTuple::Item<int> m_id_wafer;
+ NTuple::Item<int> m_id_cellX;
+ NTuple::Item<int> m_id_cellY;
+};
+
+#endif
diff --git a/Utilities/DataHelper/src/HelixClass.cc b/Utilities/DataHelper/src/HelixClass.cc
index 591acc51..8075c062 100644
--- a/Utilities/DataHelper/src/HelixClass.cc
+++ b/Utilities/DataHelper/src/HelixClass.cc
@@ -49,6 +49,8 @@ void HelixClass::Initialize_VP(float * pos, float * mom, float q, float B) {
else {
d0 = double(q)*radius + double(sqrt(xCentre*xCentre+yCentre*yCentre));
}
+ std::cout<<__FILE__<<" "<<__LINE__<<" d0 "<<d0<<" r "<<radius<<" q "<<q<<" xCentre "<<xCentre<<" yCentre "<<yCentre<<std::endl;
+ std::cout<<__FILE__<<" "<<__LINE__<<" pxy "<<_pxy<<" FCT "<<_FCT<<" B "<<B<<std::endl;
_d0 = float(d0);
--
GitLab