From 5291b81b80e8517eac2191dc7bd4caa3bb5d4ba5 Mon Sep 17 00:00:00 2001
From: Ke Li <ke.li@cern.ch>
Date: Wed, 15 Sep 2021 19:28:35 +0200
Subject: [PATCH] update the TrackerSPFit

---
 .../TrackerSPFit/CMakeLists.txt               |   3 +-
 .../TrackerSPFit/python/TrackerSPFit.py       |  94 ++++++++
 .../python/TrackerSPFitTI12Data.py            | 102 ++++++++
 .../TrackerSPFit/src/TrackerSPFit.cxx         | 227 +++++++++++++++---
 .../TrackerSPFit/src/TrackerSPFit.h           |  62 ++++-
 .../src/components/TrackerSPFit_entries.cxx   |   2 +-
 6 files changed, 436 insertions(+), 54 deletions(-)
 create mode 100644 Tracker/TrackerRecAlgs/TrackerSPFit/python/TrackerSPFit.py
 create mode 100644 Tracker/TrackerRecAlgs/TrackerSPFit/python/TrackerSPFitTI12Data.py

diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/CMakeLists.txt b/Tracker/TrackerRecAlgs/TrackerSPFit/CMakeLists.txt
index 01df7b9e..665c1a29 100644
--- a/Tracker/TrackerRecAlgs/TrackerSPFit/CMakeLists.txt
+++ b/Tracker/TrackerRecAlgs/TrackerSPFit/CMakeLists.txt
@@ -1,5 +1,5 @@
 ################################################################################
-# Package: TrackerSpacePointFormation
+# Package: TrackerSPFit
 ################################################################################
 
 # Declare the package name:
@@ -16,4 +16,3 @@ atlas_add_component( TrackerSPFit
 		     LINK_LIBRARIES ${EIGEN_LIBRARIES} ${ROOT_LIBRARIES} AthenaBaseComps AthContainers GeoPrimitives Identifier GaudiKernel TrackerReadoutGeometry TrackerPrepRawData TrackerSpacePoint AthenaMonitoringKernelLib FaserDetDescr xAODEventInfo TrkEventPrimitives TrackerIdentifier FaserSiSpacePointToolLib TrkSpacePoint )
 
 atlas_install_python_modules( python/*.py )
-atlas_install_scripts( test/*.py )
diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/python/TrackerSPFit.py b/Tracker/TrackerRecAlgs/TrackerSPFit/python/TrackerSPFit.py
new file mode 100644
index 00000000..cf69e95b
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/TrackerSPFit/python/TrackerSPFit.py
@@ -0,0 +1,94 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+import sys
+from AthenaCommon.Logging import log, logging
+from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+# from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+# from Digitization.DigitizationParametersConfig import writeDigitizationMetadata
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+
+Tracker__TrackerSPFit, THistSvc=CompFactory.getComps("Tracker::TrackerSPFit", "THistSvc")
+
+
+def TrackerSPFitBasicCfg(flags, **kwargs):
+    """Return ComponentAccumulator for TrackerSPFit"""
+    acc = FaserSCT_GeometryCfg(flags)
+    kwargs.setdefault("SpacePointsSCTName", "SCT_SpacePointContainer")
+    acc.addEventAlgo(Tracker__TrackerSPFit(**kwargs))
+   # attach ToolHandles
+    return acc
+
+def TrackerSPFit_OutputCfg(flags):
+    """Return ComponentAccumulator with Output for SCT. Not standalone."""
+    acc = ComponentAccumulator()
+    acc.merge(OutputStreamCfg(flags, "ESD"))
+    ostream = acc.getEventAlgo("OutputStreamESD")
+    ostream.TakeItemsFromInput = True
+    return acc
+
+def TrackerSPFitCfg(flags, **kwargs):
+    acc=TrackerSPFitBasicCfg(flags, **kwargs)
+    histSvc= THistSvc()
+    histSvc.Output += [ "TrackerSPFit DATAFILE='trackerspfit.root' OPT='RECREATE'" ]
+    acc.addService(histSvc)
+    acc.merge(TrackerSPFit_OutputCfg(flags))
+    return acc
+
+if __name__ == "__main__":
+  log.setLevel(DEBUG)
+  Configurable.configurableRun3Behavior = True
+
+  # Configure
+  ConfigFlags.Input.Files = ['my.RDO.pool.root']
+  ConfigFlags.Output.ESDFileName = "mySeeds.ESD.pool.root"
+  ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"             # Always needed; must match FaserVersion
+  ConfigFlags.IOVDb.DatabaseInstance = "OFLP200"               # Use MC conditions for now
+#  ConfigFlags.Input.ProjectName = "data20"                     # Needed to bypass autoconfig
+  ConfigFlags.Input.isMC = True
+  # Needed to bypass autoconfig
+  ConfigFlags.GeoModel.FaserVersion     = "FASER-01"           # FASER cosmic ray geometry (station 2 only)
+  ConfigFlags.Common.isOnline = False
+  ConfigFlags.GeoModel.Align.Dynamic = False
+  ConfigFlags.Beam.NumberOfCollisions = 0.
+
+  ConfigFlags.lock()
+
+  # Core components
+  acc = MainServicesCfg(ConfigFlags)
+  acc.merge(PoolReadCfg(ConfigFlags))
+
+  #acc.merge(writeDigitizationMetadata(ConfigFlags))
+
+  # Inner Detector
+  acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
+  acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+  acc.merge(TrackerSPFitCfg(ConfigFlags))
+
+  # Timing
+  #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+  # Dump config
+  logging.getLogger('forcomps').setLevel(VERBOSE)
+  acc.foreach_component("*").OutputLevel = VERBOSE
+  acc.foreach_component("*ClassID*").OutputLevel = INFO
+  # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+  # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+  acc.getService("StoreGateSvc").Dump = True
+  acc.getService("ConditionStore").Dump = True
+  acc.printConfig(withDetails=True)
+  ConfigFlags.dump()
+
+  # Execute and finish
+  sc = acc.run(maxEvents=-1)
+
+  # Success should be 0
+  sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/python/TrackerSPFitTI12Data.py b/Tracker/TrackerRecAlgs/TrackerSPFit/python/TrackerSPFitTI12Data.py
new file mode 100644
index 00000000..3f21c7c8
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/TrackerSPFit/python/TrackerSPFitTI12Data.py
@@ -0,0 +1,102 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+import sys
+from AthenaCommon.Logging import log, logging
+from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+
+Tracker__TrackerSPFit, THistSvc=CompFactory.getComps("Tracker::TrackerSPFit", "THistSvc")
+
+
+def TrackerSPFitBasicCfg(flags, **kwargs):
+    """Return ComponentAccumulator for TrackerSPFit"""
+    acc = FaserSCT_GeometryCfg(flags)
+    kwargs.setdefault("SpacePointsSCTName", "SCT_SpacePointContainer")
+    acc.addEventAlgo(Tracker__TrackerSPFit(**kwargs))
+   # attach ToolHandles
+    return acc
+
+def TrackerSPFit_OutputCfg(flags):
+    """Return ComponentAccumulator with Output for SCT. Not standalone."""
+    acc = ComponentAccumulator()
+    acc.merge(OutputStreamCfg(flags, "ESD"))
+    ostream = acc.getEventAlgo("OutputStreamESD")
+    ostream.TakeItemsFromInput = True
+    return acc
+
+def TrackerSPFitCfg(flags, **kwargs):
+    acc=TrackerSPFitBasicCfg(flags, **kwargs)
+    histSvc= THistSvc()
+    histSvc.Output += [ "TrackerSPFit DATAFILE='trackerspfit.root' OPT='RECREATE'" ]
+    acc.addService(histSvc)
+    acc.merge(TrackerSPFit_OutputCfg(flags))
+    return acc
+
+if __name__ == "__main__":
+  log.setLevel(DEBUG)
+  Configurable.configurableRun3Behavior = True
+
+  # Configure
+  ConfigFlags.Input.Files = ['/eos/project-f/faser-commissioning/TI12Data/Run-004272/Faser-Physics-004272-00010.raw']
+  ConfigFlags.Output.ESDFileName = "mySeeds.ESD.pool.root"
+  ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"             # Always needed; must match FaserVersion
+  ConfigFlags.IOVDb.DatabaseInstance = "OFLP200"               # Use MC conditions for now
+  ConfigFlags.Input.ProjectName = "data21"                     # Needed to bypass autoconfig
+  ConfigFlags.Input.isMC = False
+  # Needed to bypass autoconfig
+  ConfigFlags.GeoModel.FaserVersion     = "FASER-01"           # FASER cosmic ray geometry (station 2 only)
+  ConfigFlags.Common.isOnline = False
+  ConfigFlags.GeoModel.Align.Dynamic = False
+  ConfigFlags.Beam.NumberOfCollisions = 0.
+
+  ConfigFlags.lock()
+
+  # Core components
+  acc = MainServicesCfg(ConfigFlags)
+  acc.merge(PoolWriteCfg(ConfigFlags))
+
+  from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnvSvcCfg
+  acc.merge(FaserByteStreamCnvSvcCfg(ConfigFlags))
+  acc.merge(WaveformReconstructionCfg(ConfigFlags))
+  acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="SCT_EDGEMODE_RDOs"))
+#  acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
+  acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+  acc.merge(TrackerSPFitCfg(ConfigFlags))
+#  from AthenaConfiguration.ComponentFactory import CompFactory
+#  decoderTool = CompFactory.RawWaveformDecoderTool(name = "RawWaveformDecoderTool", 
+##      CaloChannels = [0, 1, 2,3, 4, 5], 
+##    PreshowerChannels = [6,7], 
+##    TriggerChannels = [8, 9],
+#      VetoChannels=[])
+#  acc.addPublicTool(decoderTool)
+
+  # Timing
+  #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+  # Dump config
+  logging.getLogger('forcomps').setLevel(INFO)
+  acc.foreach_component("*").OutputLevel = INFO
+  acc.foreach_component("*ClassID*").OutputLevel = INFO
+  # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+  # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+  acc.getService("StoreGateSvc").Dump = True
+  acc.getService("ConditionStore").Dump = True
+  acc.printConfig(withDetails=True)
+  ConfigFlags.dump()
+
+  # Execute and finish
+  sc = acc.run(maxEvents=-1)
+
+  # Success should be 0
+  sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx
index 738310f9..8b75e100 100755
--- a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx
+++ b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx
@@ -71,9 +71,9 @@ TrackerSPFit::TrackerSPFit(const std::string& name,
   , m_residual_x_plane6(0)
   , m_residual_x_plane7(0)
   , m_residual_x_plane8(0)
-  , m_chi2_station0(0)
-  , m_chi2_station1(0)
-  , m_chi2_station2(0)
+  , m_chi2(0)
+  , m_edm(0)
+  , m_ndf(0)
   , m_thistSvc("THistSvc", name)
 { 
 }
@@ -89,6 +89,8 @@ StatusCode TrackerSPFit::initialize()
     ATH_MSG_FATAL( "SCTs selected and no name set for SCT clusters");
     return StatusCode::FAILURE;
   }
+  m_tree= new TTree("spfit","tree");
+  TrackerSPFit::initializeTree();
   ATH_CHECK( m_Sct_spcontainerKey.initialize() );
 
   // create containers (requires the Identifier Helpers)
@@ -116,27 +118,28 @@ StatusCode TrackerSPFit::initialize()
   m_hist_x_y_plane6=new TH2D("sp_x_y_plane6","sp_x_y_plane6",100,-200,200,100,-200,200);
   m_hist_x_y_plane7=new TH2D("sp_x_y_plane7","sp_x_y_plane7",100,-200,200,100,-200,200);
   m_hist_x_y_plane8=new TH2D("sp_x_y_plane8","sp_x_y_plane8",100,-200,200,100,-200,200);
-  m_residual_y_plane0=new TH1D("sp_residual_y_plane0","sp_residual_y_plane0",100,-100,100);
-  m_residual_y_plane1=new TH1D("sp_residual_y_plane1","sp_residual_y_plane1",100,-100,100);
-  m_residual_y_plane2=new TH1D("sp_residual_y_plane2","sp_residual_y_plane2",100,-100,100);
-  m_residual_y_plane3=new TH1D("sp_residual_y_plane3","sp_residual_y_plane3",100,-100,100);
-  m_residual_y_plane4=new TH1D("sp_residual_y_plane4","sp_residual_y_plane4",100,-100,100);
-  m_residual_y_plane5=new TH1D("sp_residual_y_plane5","sp_residual_y_plane5",100,-100,100);
-  m_residual_y_plane6=new TH1D("sp_residual_y_plane6","sp_residual_y_plane6",100,-100,100);
-  m_residual_y_plane7=new TH1D("sp_residual_y_plane7","sp_residual_y_plane7",100,-100,100);
-  m_residual_y_plane8=new TH1D("sp_residual_y_plane8","sp_residual_y_plane8",100,-100,100);
-  m_residual_x_plane0=new TH1D("sp_residual_x_plane0","sp_residual_x_plane0",100,-100,100);
-  m_residual_x_plane1=new TH1D("sp_residual_x_plane1","sp_residual_x_plane1",100,-100,100);
-  m_residual_x_plane2=new TH1D("sp_residual_x_plane2","sp_residual_x_plane2",100,-100,100);
-  m_residual_x_plane3=new TH1D("sp_residual_x_plane3","sp_residual_x_plane3",100,-100,100);
-  m_residual_x_plane4=new TH1D("sp_residual_x_plane4","sp_residual_x_plane4",100,-100,100);
-  m_residual_x_plane5=new TH1D("sp_residual_x_plane5","sp_residual_x_plane5",100,-100,100);
-  m_residual_x_plane6=new TH1D("sp_residual_x_plane6","sp_residual_x_plane6",100,-100,100);
-  m_residual_x_plane7=new TH1D("sp_residual_x_plane7","sp_residual_x_plane7",100,-100,100);
-  m_residual_x_plane8=new TH1D("sp_residual_x_plane8","sp_residual_x_plane8",100,-100,100);
-  m_chi2_station0=new TH1D("chi2_station0","chi2_station0",100,0,100);
-  m_chi2_station1=new TH1D("chi2_station1","chi2_station1",100,0,100);
-  m_chi2_station2=new TH1D("chi2_station2","chi2_station2",100,0,100);
+  m_residual_y_plane0=new TH1D("sp_residual_y_plane0","sp_residual_y_plane0",2000,-1,1);
+  m_residual_y_plane1=new TH1D("sp_residual_y_plane1","sp_residual_y_plane1",2000,-1,1);
+  m_residual_y_plane2=new TH1D("sp_residual_y_plane2","sp_residual_y_plane2",2000,-1,1);
+  m_residual_y_plane3=new TH1D("sp_residual_y_plane3","sp_residual_y_plane3",2000,-1,1);
+  m_residual_y_plane4=new TH1D("sp_residual_y_plane4","sp_residual_y_plane4",2000,-1,1);
+  m_residual_y_plane5=new TH1D("sp_residual_y_plane5","sp_residual_y_plane5",2000,-1,1);
+  m_residual_y_plane6=new TH1D("sp_residual_y_plane6","sp_residual_y_plane6",2000,-1,1);
+  m_residual_y_plane7=new TH1D("sp_residual_y_plane7","sp_residual_y_plane7",2000,-1,1);
+  m_residual_y_plane8=new TH1D("sp_residual_y_plane8","sp_residual_y_plane8",2000,-1,1);
+  m_residual_x_plane0=new TH1D("sp_residual_x_plane0","sp_residual_x_plane0",2000,-1,1);
+  m_residual_x_plane1=new TH1D("sp_residual_x_plane1","sp_residual_x_plane1",2000,-1,1);
+  m_residual_x_plane2=new TH1D("sp_residual_x_plane2","sp_residual_x_plane2",2000,-1,1);
+  m_residual_x_plane3=new TH1D("sp_residual_x_plane3","sp_residual_x_plane3",2000,-1,1);
+  m_residual_x_plane4=new TH1D("sp_residual_x_plane4","sp_residual_x_plane4",2000,-1,1);
+  m_residual_x_plane5=new TH1D("sp_residual_x_plane5","sp_residual_x_plane5",2000,-1,1);
+  m_residual_x_plane6=new TH1D("sp_residual_x_plane6","sp_residual_x_plane6",2000,-1,1);
+  m_residual_x_plane7=new TH1D("sp_residual_x_plane7","sp_residual_x_plane7",2000,-1,1);
+  m_residual_x_plane8=new TH1D("sp_residual_x_plane8","sp_residual_x_plane8",2000,-1,1);
+  m_chi2=new TH1D("chi2","chi2",100,0,100);
+  m_edm=new TH1D("edm","edm",100,0,0.01);
+  m_ndf=new TH1D("ndf","ndf",10,0,10);
+  CHECK(m_thistSvc->regTree("/TrackerSPFit/spfit",m_tree));
   CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_x",m_hist_x));
   CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_y",m_hist_y));
   CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_z",m_hist_z));
@@ -173,20 +176,70 @@ StatusCode TrackerSPFit::initialize()
   CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_residual_x_plane6",m_residual_x_plane6));
   CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_residual_x_plane7",m_residual_x_plane7));
   CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_residual_x_plane8",m_residual_x_plane8));
-  CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_chi2_station0",m_chi2_station0));
-  CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_chi2_station1",m_chi2_station1));
-  CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_chi2_station2",m_chi2_station2));
+  CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_chi2",m_chi2));
+  CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_edm",m_edm));
+  CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_ndf",m_ndf));
   ATH_MSG_INFO( "TrackerSPFit::initialized" );
   return StatusCode::SUCCESS;
 }
 
+//-------------------------------------------------------------------------
+void TrackerSPFit::initializeTree(){
+  m_tree->Branch("evtId",&m_eventNumber);
+  m_tree->Branch("sp_x",&m_sp_x);
+  m_tree->Branch("sp_y",&m_sp_y);
+  m_tree->Branch("sp_z",&m_sp_z);
+  m_tree->Branch("sp_x_err",&m_sp_x_err);
+  m_tree->Branch("sp_y_err",&m_sp_y_err);
+  m_tree->Branch("sp_z_err",&m_sp_z_err);
+  m_tree->Branch("sp_x_predicted",&m_sp_x_predicted);
+  m_tree->Branch("sp_y_predicted",&m_sp_y_predicted);
+  m_tree->Branch("sp_z_predicted",&m_sp_z_predicted);
+  m_tree->Branch("sp_x_residual",&m_sp_x_residual);
+  m_tree->Branch("sp_y_residual",&m_sp_y_residual);
+  m_tree->Branch("sp_z_residual",&m_sp_z_residual);
+  m_tree->Branch("trackId",&m_trackId);
+  m_tree->Branch("track_chi2",&m_track_chi2);
+  m_tree->Branch("track_edm",&m_track_edm);
+  m_tree->Branch("track_ndf",&m_track_ndf);
+  m_tree->Branch("track_p0",&m_track_p0);
+  m_tree->Branch("track_p1",&m_track_p1);
+  m_tree->Branch("track_p2",&m_track_p2);
+  m_tree->Branch("track_p3",&m_track_p3);
+  m_tree->Branch("sp_station",&m_sp_station);
+  m_tree->Branch("sp_layer",&m_sp_layer);
+  m_tree->Branch("sp_module",&m_sp_module);
+}
+//-------------------------------------------------------------------------
+void TrackerSPFit::clearVariables() const {
+  m_sp_x.clear();
+  m_sp_y.clear();
+  m_sp_z.clear();
+  m_sp_x_err.clear();
+  m_sp_y_err.clear();
+  m_sp_z_err.clear();
+  m_sp_x_predicted.clear();
+  m_sp_y_predicted.clear();
+  m_sp_z_predicted.clear();
+  m_sp_x_residual.clear();
+  m_sp_y_residual.clear();
+  m_sp_z_residual.clear();
+  m_trackId.clear();
+  m_track_chi2.clear();
+  m_track_edm.clear();
+  m_track_ndf.clear();
+  m_sp_station.clear();
+  m_sp_layer.clear();
+  m_sp_module.clear();
+}
 //-------------------------------------------------------------------------
 
 StatusCode TrackerSPFit::execute (const EventContext& ctx) const
 {
 
-
   ++m_numberOfEvents;
+  m_eventNumber=m_numberOfEvents;
+  clearVariables();
   const TrackerDD::SiDetectorElementCollection* elements = nullptr;
 
 
@@ -227,6 +280,10 @@ StatusCode TrackerSPFit::execute (const EventContext& ctx) const
     Identifier elementID = colNext->identify();
     int istation = m_idHelper->station(elementID);
     int ilayer = m_idHelper->layer(elementID);
+    int ietamodule = m_idHelper->eta_module(elementID);
+    int iphimodule = m_idHelper->phi_module(elementID);
+    int imodule=iphimodule;
+    if(ietamodule<0)imodule+=4;
     m_hist_strip->Fill(m_idHelper->strip(elementID));
     m_hist_station->Fill(istation);
     m_hist_layer->Fill(ilayer);
@@ -246,13 +303,13 @@ StatusCode TrackerSPFit::execute (const EventContext& ctx) const
 	m_hist_x->Fill(gloPos.x());
 	m_hist_y->Fill(gloPos.y());
 	m_hist_z->Fill(gloPos.z());
-	  struct SP_Seed tmp{gloPos,sp->globCovariance(),ilayer};
+	struct SP_Seed tmp{gloPos,sp->globCovariance(),ilayer,imodule};
 	if(istation==1)
-	sp_sta0.push_back(tmp);
+	  sp_sta0.push_back(tmp);
 	if(istation==2)
-	sp_sta1.push_back(tmp);
+	  sp_sta1.push_back(tmp);
 	if(istation==3)
-	sp_sta2.push_back(tmp);
+	  sp_sta2.push_back(tmp);
 
 	if ( ((istation-1)*3+ilayer) == 0 ) m_hist_x_y_plane0->Fill(gloPos.x(),gloPos.y());
 	if ( ((istation-1)*3+ilayer) == 1 ) m_hist_x_y_plane1->Fill(gloPos.x(),gloPos.y());
@@ -267,14 +324,43 @@ StatusCode TrackerSPFit::execute (const EventContext& ctx) const
       ATH_MSG_VERBOSE( size << " SpacePoints successfully added to Container !" );
     }
   }
+  ATH_MSG_DEBUG( "TrackerSPFit number of spacepoints in each stations: "<<sp_sta0.size()<<" "<<sp_sta1.size()<<" "<<sp_sta2.size() );
 
   double maxchi2=10.;
+  bool ifTrack=false;
   if(sp_sta0.size()>2) {
     auto tracks_sta0=makeTrackSeg(sp_sta0,maxchi2);
+  ATH_MSG_DEBUG( "TrackerSPFit found "<<tracks_sta0.size()<<" track segments in station 0 " );
     if(tracks_sta0.size()>0){
+      ifTrack=true;
       for(auto track:tracks_sta0){
 	auto pre=track.predicted;
 	auto pos=track.pos;
+	++m_numberOfTrack;
+	m_trackId.push_back(int(m_numberOfTrack/3));
+	m_sp_x.push_back(pos.x());
+	m_sp_y.push_back(pos.y());
+	m_sp_z.push_back(pos.z());
+	m_sp_x_err.push_back(track.err.x());
+	m_sp_y_err.push_back(track.err.y());
+	m_sp_z_err.push_back(track.err.z());
+	m_sp_x_predicted.push_back(pre.x());
+	m_sp_y_predicted.push_back(pre.y());
+	m_sp_z_predicted.push_back(pre.z());
+	m_sp_x_residual.push_back(pre.x()-pos.x());
+	m_sp_y_residual.push_back(pre.y()-pos.y());
+	m_sp_z_residual.push_back(pre.z()-pos.z());
+	m_sp_station.push_back(0);
+	m_sp_layer.push_back(track.layer);
+	m_sp_module.push_back(track.module);
+	m_track_chi2.push_back(track.chi2);
+	m_track_edm.push_back(track.edm);
+	m_track_ndf.push_back(track.ndf);
+	double* par=const_cast<double*>(track.p);
+	m_track_p0.push_back(par[0]);
+	m_track_p1.push_back(par[1]);
+	m_track_p2.push_back(par[2]);
+	m_track_p3.push_back(par[3]);
 	if(track.layer==0){ m_residual_x_plane0->Fill(pre.x()-pos.x());m_residual_y_plane0->Fill(pre.y()-pos.y());}
 	if(track.layer==1){ m_residual_x_plane1->Fill(pre.x()-pos.x());m_residual_y_plane1->Fill(pre.y()-pos.y());}
 	if(track.layer==2){ m_residual_x_plane2->Fill(pre.x()-pos.x());m_residual_y_plane2->Fill(pre.y()-pos.y());}
@@ -283,10 +369,37 @@ StatusCode TrackerSPFit::execute (const EventContext& ctx) const
   }
   if(sp_sta1.size()>2) {
   auto tracks_sta1=makeTrackSeg(sp_sta1,maxchi2);
+  ATH_MSG_DEBUG( "TrackerSPFit found "<<tracks_sta1.size()<<" track segments in station 1 " );
   if(tracks_sta1.size()>0){
+      ifTrack=true;
     for(auto track:tracks_sta1){
       auto pre=track.predicted;
       auto pos=track.pos;
+	++m_numberOfTrack;
+	m_trackId.push_back(int(m_numberOfTrack/3));
+	m_sp_x.push_back(pos.x());
+	m_sp_y.push_back(pos.y());
+	m_sp_z.push_back(pos.z());
+	m_sp_x_err.push_back(track.err.x());
+	m_sp_y_err.push_back(track.err.y());
+	m_sp_z_err.push_back(track.err.z());
+	m_sp_x_predicted.push_back(pre.x());
+	m_sp_y_predicted.push_back(pre.y());
+	m_sp_z_predicted.push_back(pre.z());
+	m_sp_x_residual.push_back(pre.x()-pos.x());
+	m_sp_y_residual.push_back(pre.y()-pos.y());
+	m_sp_z_residual.push_back(pre.z()-pos.z());
+	m_sp_station.push_back(1);
+	m_sp_layer.push_back(track.layer);
+	m_sp_module.push_back(track.module);
+	m_track_chi2.push_back(track.chi2);
+	m_track_edm.push_back(track.edm);
+	m_track_ndf.push_back(track.ndf);
+	double* par=const_cast<double*>(track.p);
+	m_track_p0.push_back(par[0]);
+	m_track_p1.push_back(par[1]);
+	m_track_p2.push_back(par[2]);
+	m_track_p3.push_back(par[3]);
       if(track.layer==0){ m_residual_x_plane3->Fill(pre.x()-pos.x());m_residual_y_plane3->Fill(pre.y()-pos.y());}
       if(track.layer==1){ m_residual_x_plane4->Fill(pre.x()-pos.x());m_residual_y_plane4->Fill(pre.y()-pos.y());}
       if(track.layer==2){ m_residual_x_plane5->Fill(pre.x()-pos.x());m_residual_y_plane5->Fill(pre.y()-pos.y());}
@@ -295,10 +408,37 @@ StatusCode TrackerSPFit::execute (const EventContext& ctx) const
   }
   if(sp_sta2.size()>2) {
   auto tracks_sta2=makeTrackSeg(sp_sta2,maxchi2);
+  ATH_MSG_DEBUG( "TrackerSPFit found "<<tracks_sta2.size()<<" track segments in station 2 " );
   if(tracks_sta2.size()>0){
+      ifTrack=true;
     for(auto track:tracks_sta2){
       auto pre=track.predicted;
       auto pos=track.pos;
+	++m_numberOfTrack;
+	m_trackId.push_back(int(m_numberOfTrack/3));
+	m_sp_x.push_back(pos.x());
+	m_sp_y.push_back(pos.y());
+	m_sp_z.push_back(pos.z());
+	m_sp_x_err.push_back(track.err.x());
+	m_sp_y_err.push_back(track.err.y());
+	m_sp_z_err.push_back(track.err.z());
+	m_sp_x_predicted.push_back(pre.x());
+	m_sp_y_predicted.push_back(pre.y());
+	m_sp_z_predicted.push_back(pre.z());
+	m_sp_x_residual.push_back(pre.x()-pos.x());
+	m_sp_y_residual.push_back(pre.y()-pos.y());
+	m_sp_z_residual.push_back(pre.z()-pos.z());
+	m_sp_station.push_back(2);
+	m_sp_layer.push_back(track.layer);
+	m_sp_module.push_back(track.module);
+	m_track_chi2.push_back(track.chi2);
+	m_track_edm.push_back(track.edm);
+	m_track_ndf.push_back(track.ndf);
+	double* par=const_cast<double*>(track.p);
+	m_track_p0.push_back(par[0]);
+	m_track_p1.push_back(par[1]);
+	m_track_p2.push_back(par[2]);
+	m_track_p3.push_back(par[3]);
       if(track.layer==0){ m_residual_x_plane6->Fill(pre.x()-pos.x());m_residual_y_plane6->Fill(pre.y()-pos.y());}
       if(track.layer==1){ m_residual_x_plane7->Fill(pre.x()-pos.x());m_residual_y_plane7->Fill(pre.y()-pos.y());}
       if(track.layer==2){ m_residual_x_plane8->Fill(pre.x()-pos.x());m_residual_y_plane8->Fill(pre.y()-pos.y());}
@@ -306,7 +446,8 @@ StatusCode TrackerSPFit::execute (const EventContext& ctx) const
   }
   }
 
-
+  if(ifTrack)
+  m_tree->Fill();
   return StatusCode::SUCCESS;
 }
 
@@ -365,13 +506,23 @@ std::vector<TrackerSPFit::SP_TSOS> TrackerSPFit::makeTrackSeg(SP_Seed sp0, SP_Se
   //   for (int i = 0; i < 4; ++i) fitter.Config().ParSettings(i).SetStepSize(0.01);
   if(fitter.FitFCN()){
       const ROOT::Fit::FitResult & result =fitter.Result();
-    m_chi2_station0->Fill(result.Chi2());
-    if(result.Chi2()<maxchi2){
+      double chi2=result.MinFcnValue();
+      double edm=result.Edm();
+      int ndf=result.Ndf();
+  ATH_MSG_DEBUG( "TrackerSPFit::makeTrackSeg(), track chi2 = "<<chi2<<"  ; edm = "<<edm<<" ; ndf = "<<ndf );
+    m_chi2->Fill(chi2);
+    m_edm->Fill(edm);
+    m_ndf->Fill(ndf);
+
+    if(chi2<maxchi2){
       //result.Print(std::cout);
       const double *fitParam=result.GetParams();
-      struct SP_TSOS tsos0{sp0.pos,predicted(sp0.pos.z(),fitParam),sp0.layer};
-      struct SP_TSOS tsos1{sp1.pos,predicted(sp1.pos.z(),fitParam),sp1.layer};
-      struct SP_TSOS tsos2{sp2.pos,predicted(sp2.pos.z(),fitParam),sp2.layer};
+      Amg::Vector3D err0(sqrt((sp0.cov)(0,0)),sqrt((sp0.cov)(1,1)),sqrt((sp0.cov)(2,2)));
+      Amg::Vector3D err1(sqrt((sp1.cov)(0,0)),sqrt((sp1.cov)(1,1)),sqrt((sp1.cov)(2,2)));
+      Amg::Vector3D err2(sqrt((sp2.cov)(0,0)),sqrt((sp2.cov)(1,1)),sqrt((sp2.cov)(2,2)));
+      struct SP_TSOS tsos0{sp0.pos,predicted(sp0.pos.z(),fitParam),err0,sp0.layer,sp0.module,chi2,edm,ndf,fitParam};
+      struct SP_TSOS tsos1{sp1.pos,predicted(sp1.pos.z(),fitParam),err1,sp1.layer,sp0.module,chi2,edm,ndf,fitParam};
+      struct SP_TSOS tsos2{sp2.pos,predicted(sp2.pos.z(),fitParam),err2,sp2.layer,sp0.module,chi2,edm,ndf,fitParam};
       spt.push_back(tsos0);
       spt.push_back(tsos1);
       spt.push_back(tsos2);
diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h
index 9ff147ef..8ae562fb 100755
--- a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h
+++ b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h
@@ -27,6 +27,7 @@
 #include <vector>
 #include "TH1.h"
 #include "TH2.h"
+#include "TTree.h"
 #include "TGraph2DErrors.h"
 #include "TMath.h"
 #include "Math/Functor.h"
@@ -58,8 +59,8 @@ namespace Tracker
 
     virtual StatusCode finalize() override;
 
-    struct SP_Seed{Amg::Vector3D pos; Amg::MatrixX cov; int layer;};
-    struct SP_TSOS{Amg::Vector3D pos; Amg::Vector3D predicted; int layer;};
+    struct SP_Seed{Amg::Vector3D pos; Amg::MatrixX cov; int layer; int module;};
+    struct SP_TSOS{Amg::Vector3D pos; Amg::Vector3D predicted; Amg::Vector3D err; int layer;int module; double chi2; double edm; int ndf; const double *p;};
     std::vector<SP_TSOS> makeTrackSeg(std::vector<SP_Seed> sps, double maxchi2) const;
     std::vector<SP_TSOS> makeTrackSeg(SP_Seed sp0, SP_Seed sp1, SP_Seed sp2, double maxchi2) const;
     Amg::Vector3D predicted(double z, const double *p) const;
@@ -70,14 +71,19 @@ namespace Tracker
       TGraph2DErrors *fGraph;
       SumDistance2(TGraph2DErrors *g) : fGraph(g) {}
       // calculate distance line-point
-      double distance2(double x,double y,double z, const double *p) {
+      double distance2(double x,double y,double z, double xe, double ye, const double *p) {
+	double xpred=p[0]+p[1]*z;
+	double ypred=p[2]+p[3]*z;
+	double xchi2= (xpred-x)*(xpred-x)/xe/xe;
+	double ychi2= (ypred-y)*(ypred-y)/ye/ye;
+	double d2=xchi2+ychi2;
 	// distance line point is D= | (xp-x0) cross  ux |
 	// where ux is direction of line and x0 is a point in the line (like t = 0)
-	ROOT::Math::XYZVector xp(x,y,z);
-	ROOT::Math::XYZVector x0(p[0], p[2], 0. );
-	ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. );
-	ROOT::Math::XYZVector u = (x1-x0).Unit();
-	double d2 = ((xp-x0).Cross(u)).Mag2();
+//	ROOT::Math::XYZVector xp(x,y,z);
+//	ROOT::Math::XYZVector x0(p[0], p[2], 0. );
+//	ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. );
+//	ROOT::Math::XYZVector u = (x1-x0).Unit();
+//	double d2 = ((xp-x0).Cross(u)).Mag2();
 	return d2;
       }
       // implementation of the function to be minimized
@@ -86,10 +92,12 @@ namespace Tracker
 	double * x = fGraph->GetX();
 	double * y = fGraph->GetY();
 	double * z = fGraph->GetZ();
+	double * xe = fGraph->GetEX();
+	double * ye = fGraph->GetEY();
 	int npoints = fGraph->GetN();
 	double sum = 0;
 	for (int i  = 0; i < npoints; ++i) {
-	  double d = distance2(x[i],y[i],z[i],par);
+	  double d = distance2(x[i],y[i],z[i],xe[i],ye[i],par);
 	  sum += d;
 	}
 	return sum;
@@ -97,12 +105,15 @@ namespace Tracker
     };
 
   private:
+    void initializeTree();
+    void clearVariables() const;
+    mutable TTree* m_tree;
     TrackerSPFit() = delete;
     TrackerSPFit(const TrackerSPFit&) =delete;
     TrackerSPFit &operator=(const TrackerSPFit&) = delete;
     //@}
 
-    SG::ReadHandleKey<FaserSCT_SpacePointContainer>  m_Sct_spcontainerKey{this, "SpacePointsSCTName", "SCT spContainer"};
+    SG::ReadHandleKey<FaserSCT_SpacePointContainer>  m_Sct_spcontainerKey{this, "SpacePointsSCTName", "SCT_SpacePointContainer"};
     //@}
 
     SG::ReadCondHandleKey<TrackerDD::SiDetectorElementCollection> m_SCTDetEleCollKey{this, "SCTDetEleCollKey", "SCT_DetectorElementCollection", "Key of SiDetectorElementCollection for SCT"};
@@ -111,6 +122,7 @@ namespace Tracker
     const FaserSCT_ID* m_idHelper{nullptr};
     mutable std::atomic<int> m_numberOfEvents{0};
     mutable std::atomic<int> m_numberOfSCT{0};
+    mutable std::atomic<int> m_numberOfTrack{0};
     mutable std::atomic<int> m_numberOfSP{0};
     //@}
    TH1* m_hist_x;
@@ -149,11 +161,35 @@ namespace Tracker
    TH1* m_residual_x_plane6;
    TH1* m_residual_x_plane7;
    TH1* m_residual_x_plane8;
-   TH1* m_chi2_station0;
-   TH1* m_chi2_station1;
-   TH1* m_chi2_station2;
+   TH1* m_chi2;
+   TH1* m_edm;
+   TH1* m_ndf;
    ServiceHandle<ITHistSvc>  m_thistSvc;
 
+   mutable int m_eventNumber;
+   mutable std::vector<double> m_sp_x;
+   mutable std::vector<double> m_sp_y;
+   mutable std::vector<double> m_sp_z;
+   mutable std::vector<double> m_sp_x_err;
+   mutable std::vector<double> m_sp_y_err;
+   mutable std::vector<double> m_sp_z_err;
+   mutable std::vector<double> m_sp_x_predicted;
+   mutable std::vector<double> m_sp_y_predicted;
+   mutable std::vector<double> m_sp_z_predicted;
+   mutable std::vector<double> m_sp_x_residual;
+   mutable std::vector<double> m_sp_y_residual;
+   mutable std::vector<double> m_sp_z_residual;
+   mutable std::vector<int> m_trackId;
+   mutable std::vector<double> m_track_chi2;
+   mutable std::vector<double> m_track_edm;
+   mutable std::vector<int> m_track_ndf;
+   mutable std::vector<double> m_track_p0;
+   mutable std::vector<double> m_track_p1;
+   mutable std::vector<double> m_track_p2;
+   mutable std::vector<double> m_track_p3;
+   mutable std::vector<int> m_sp_station;
+   mutable std::vector<int> m_sp_layer;
+   mutable std::vector<int> m_sp_module;
 
   };
 
diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/src/components/TrackerSPFit_entries.cxx b/Tracker/TrackerRecAlgs/TrackerSPFit/src/components/TrackerSPFit_entries.cxx
index 9505e880..dae3ad64 100644
--- a/Tracker/TrackerRecAlgs/TrackerSPFit/src/components/TrackerSPFit_entries.cxx
+++ b/Tracker/TrackerRecAlgs/TrackerSPFit/src/components/TrackerSPFit_entries.cxx
@@ -1,4 +1,4 @@
 #include "../TrackerSPFit.h"
 
-DECLARE_COMPONENT( Tracker::TrackerSPFit)
+DECLARE_COMPONENT( Tracker::TrackerSPFit )
 
-- 
GitLab