diff --git a/Examples/options/tut_detsim_digi_truthTracker_SDT.py b/Examples/options/tut_detsim_digi_truthTracker_SDT.py new file mode 100644 index 0000000000000000000000000000000000000000..a5fd439aa790627b439c51e0146b233b5069399a --- /dev/null +++ b/Examples/options/tut_detsim_digi_truthTracker_SDT.py @@ -0,0 +1,200 @@ +#!/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 +) diff --git a/Reconstruction/Tracking/CMakeLists.txt b/Reconstruction/Tracking/CMakeLists.txt index 44b8702ab554e7b705bfc597ae2046ae950c610b..166add75f15af42d8d3bdb8ee087fc3553108192 100644 --- a/Reconstruction/Tracking/CMakeLists.txt +++ b/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 diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp new file mode 100644 index 0000000000000000000000000000000000000000..612fae8e2a0de342361d2ad0527a5c90b9055629 --- /dev/null +++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp @@ -0,0 +1,180 @@ +#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(); +} diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..23e581e1b73a4e1f962ddc3e53cfb8e525258de6 --- /dev/null +++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h @@ -0,0 +1,65 @@ +#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