Add TruthTracerAlg for DC

Examples/options/tut_detsim_digi_truthTracker_SDT.py

+#!/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 NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = ["MyTuples DATAFILE='DCH_digi_ana.root' OPT='NEW' TYP='ROOT'"]
+
+##############################################################################
+# DCHDigiAlg
+##############################################################################
+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
+
+##############################################################################
+# TruthTrackerAlg
+##############################################################################
+from Configurables import TruthTrackerAlg
+truthTrackerAlg = TruthTrackerAlg("TruthTrackerAlg")
+truthTrackerAlg.DCHitAssociationCollection="DCHAssociationCollectio"
+truthTrackerAlg.debug = 1
+
+##############################################################################
+# POD I/O
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "truthRec_DCH.root"
+out.outputCommands = ["keep *"]
+
+##############################################################################
+# ApplicationMgr
+##############################################################################
+
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, truthTrackerAlg, out],
+                EvtSel = 'NONE',
+                EvtMax = 10,
+                ExtSvc = [rndmengine, dsvc, geosvc],
+                HistogramPersistency = "ROOT",
+                OutputLevel=INFO
+)

Reconstruction/Tracking/CMakeLists.txt

@@ -11,17 +11,20 @@ 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
+    DetSegmentation
       -Wl,--no-as-needed
      EDM4HEP::edm4hep EDM4HEP::edm4hepDict
      -Wl,--as-needed

Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp

+#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;
+    }
+
+    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();
+}

Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h

+#ifndef TruthTrackerAlg_h
+#define TruthTrackerAlg_h
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "k4FWCore/DataHandle.h"
+#include "DD4hep/Fields.h"
+
+class IGeomSvc;
+namespace dd4hep {
+    class Detector;
+    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::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"};
+
+        Gaudi::Property<int>  m_debug{ this, "debug", false};
+};
+
+#endif