diff --git a/Digitisers/DCHDigi/CMakeLists.txt b/Digitisers/DCHDigi/CMakeLists.txt index a2d07c1b835859cf4c3e4deff65fcb055cad2a09..c777aa6ff3b4b2a5b344661505e245f7132bc759 100644 --- a/Digitisers/DCHDigi/CMakeLists.txt +++ b/Digitisers/DCHDigi/CMakeLists.txt @@ -3,11 +3,10 @@ gaudi_subdir(DCHDigi v0r0) find_package(CLHEP REQUIRED;CONFIG) find_package(DD4hep COMPONENTS DDG4 REQUIRED) find_package(EDM4HEP REQUIRED ) -message("EDM4HEP_INCLUDE_DIRS: ${EDM4HEP_INCLUDE_DIR}") -message("EDM4HEP_LIB: ${EDM4HEP_LIBRARIES}") include_directories(${EDM4HEP_INCLUDE_DIR}) find_package(podio REQUIRED ) +find_package(ROOT COMPONENTS MathCore Physics GenVector Geom REQUIRED) set(srcs src/*.cpp @@ -15,11 +14,12 @@ set(srcs gaudi_depends_on_subdirs( Detector/DetInterface + Detector/DetSegmentation ) ## Modules gaudi_add_module(DCHDigi ${srcs} - INCLUDE_DIRS FWCore GaudiKernel GaudiAlgLib ${CLHEP_INCLUDE_DIR} DD4hep - LINK_LIBRARIES FWCore GaudiKernel GaudiAlgLib ${CLHEP_LIBRARIES} DD4hep ${DD4hep_COMPONENT_LIBRARIES} DDRec + INCLUDE_DIRS FWCore GaudiKernel GaudiAlgLib ${CLHEP_INCLUDE_DIR} DD4hep ROOT + LINK_LIBRARIES FWCore GaudiKernel GaudiAlgLib ${CLHEP_LIBRARIES} DD4hep ${DD4hep_COMPONENT_LIBRARIES} DDRec ROOT DetSegmentation -Wl,--no-as-needed EDM4HEP::edm4hep EDM4HEP::edm4hepDict ) diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp index 4f782b71cd79c852d974daa9d8d001be237398c6..d7c2949b4416a5d5e62809a620c3cb33f1b4d274 100644 --- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp +++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp @@ -10,6 +10,9 @@ #include <DD4hep/Objects.h> #include "DDRec/Vector3D.h" +#include "GaudiKernel/INTupleSvc.h" +#include "GaudiKernel/MsgStream.h" + #include <array> #include <math.h> @@ -29,39 +32,74 @@ DCHDigiAlg::DCHDigiAlg(const std::string& name, ISvcLocator* svcLoc) // Output collections declareProperty("DigiDCHitCollection", w_DigiDCHCol, "Handle of Digi DCHit collection"); - declareProperty("CaloAssociationCollection", w_AssociationCol, "Handle of Association collection"); + declareProperty("AssociationCollection", w_AssociationCol, "Handle of Association collection"); } StatusCode DCHDigiAlg::initialize() { - /* m_geosvc = service<IGeomSvc>("GeomSvc"); if ( !m_geosvc ) throw "DCHDigiAlg :Failed to find GeomSvc ..."; dd4hep::Detector* m_dd4hep = m_geosvc->lcdd(); if ( !m_dd4hep ) throw "DCHDigiAlg :Failed to get dd4hep::Detector ..."; - m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep); - */ + dd4hep::Readout readout = m_dd4hep->readout(m_readout_name); + m_segmentation = dynamic_cast<dd4hep::DDSegmentation::GridDriftChamber*>(readout.segmentation().segmentation()); + m_decoder = m_geosvc->getDecoder(m_readout_name); + if (!m_decoder) { + error() << "Failed to get the decoder. " << endmsg; + return StatusCode::FAILURE; + } + + if(m_WriteAna){ + + NTuplePtr nt( ntupleSvc(), "MyTuples/DCH_digi_evt" ); + if ( nt ) m_tuple = nt; + else { + m_tuple = ntupleSvc()->book( "MyTuples/DCH_digi_evt", CLID_ColumnWiseTuple, "DCH_digi_evt" ); + if ( m_tuple ) { + m_tuple->addItem( "N_sim" , m_n_sim , 0, 100000 ).ignore(); + m_tuple->addItem( "N_digi", m_n_digi, 0, 100000 ).ignore(); + m_tuple->addItem( "simhit_x", m_n_sim, m_simhit_x).ignore(); + m_tuple->addItem( "simhit_y", m_n_sim, m_simhit_y).ignore(); + m_tuple->addItem( "simhit_z", m_n_sim, m_simhit_z).ignore(); + m_tuple->addItem( "chamber" , m_n_digi, m_chamber ).ignore(); + m_tuple->addItem( "layer" , m_n_digi, m_layer ).ignore(); + m_tuple->addItem( "cell" , m_n_digi, m_cell ).ignore(); + m_tuple->addItem( "cell_x" , m_n_digi, m_cell_x ).ignore(); + m_tuple->addItem( "cell_y" , m_n_digi, m_cell_y ).ignore(); + m_tuple->addItem( "hit_x" , m_n_digi,m_hit_x ).ignore(); + m_tuple->addItem( "hit_y" , m_n_digi,m_hit_y ).ignore(); + m_tuple->addItem( "hit_z" , m_n_digi,m_hit_z ).ignore(); + m_tuple->addItem( "dca" , m_n_digi,m_dca ).ignore(); + m_tuple->addItem( "hit_dE" , m_n_digi,m_hit_dE ).ignore(); + m_tuple->addItem( "hit_dE_dx", m_n_digi,m_hit_dE_dx ).ignore(); + } + else { // did not manage to book the N tuple.... + error() << " Cannot book N-tuple:" << long( m_tuple ) << endmsg; + return StatusCode::FAILURE; + } + } + } std::cout<<"DCHDigiAlg::initialized"<< std::endl; return GaudiAlgorithm::initialize(); } + StatusCode DCHDigiAlg::execute() { + info() << "Processing " << _nEvt << " events " << endmsg; std::map<unsigned long long, std::vector<edm4hep::SimTrackerHit> > id_hits_map; edm4hep::TrackerHitCollection* Vec = w_DigiDCHCol.createAndPut(); edm4hep::MCRecoTrackerAssociationCollection* AssoVec = w_AssociationCol.createAndPut(); const edm4hep::SimTrackerHitCollection* SimHitCol = r_SimDCHCol.get(); - if(SimHitCol == 0) - { - std::cout<<"not found SimCalorimeterHitCollection"<< std::endl; - return StatusCode::SUCCESS; - } std::cout<<"input sim hit size="<< SimHitCol->size() <<std::endl; for( int i = 0; i < SimHitCol->size(); i++ ) { edm4hep::SimTrackerHit SimHit = SimHitCol->at(i); unsigned long long id = SimHit.getCellID(); + float sim_hit_mom = sqrt( SimHit.getMomentum()[0]*SimHit.getMomentum()[0] + SimHit.getMomentum()[1]*SimHit.getMomentum()[1] + SimHit.getMomentum()[2]*SimHit.getMomentum()[2] );//GeV + if(sim_hit_mom < m_mom_threshold) continue; + if(SimHit.getEDep() <= 0) continue; if ( id_hits_map.find(id) != id_hits_map.end()) id_hits_map[id].push_back(SimHit); else @@ -72,13 +110,15 @@ StatusCode DCHDigiAlg::execute() } } + m_n_sim = 0; + m_n_digi = 0 ; for(std::map<unsigned long long, std::vector<edm4hep::SimTrackerHit> >::iterator iter = id_hits_map.begin(); iter != id_hits_map.end(); iter++) { + unsigned long long wcellid = iter->first; auto trkHit = Vec->create(); - trkHit.setCellID(iter->first); + trkHit.setCellID(wcellid); double tot_edep = 0 ; double tot_length = 0 ; - double tot_time = 0 ; double pos_x = 0 ; double pos_y = 0 ; double pos_z = 0 ; @@ -87,45 +127,80 @@ StatusCode DCHDigiAlg::execute() { tot_edep += iter->second.at(i).getEDep();//GeV } + int chamber = m_decoder->get(wcellid, "chamber"); + int layer = m_decoder->get(wcellid, "layer" ); + int cellID = m_decoder->get(wcellid, "cellID" ); + TVector3 Wstart(0,0,0); + TVector3 Wend (0,0,0); + m_segmentation->cellposition(wcellid, Wstart, Wend); + Wstart = 10*Wstart;// from DD4HEP cm to mm + Wend = 10*Wend ; + //std::cout<<"wcellid="<<wcellid<<",chamber="<<chamber<<",layer="<<layer<<",cellID="<<cellID<<",s_x="<<Wstart.x()<<",s_y="<<Wstart.y()<<",s_z="<<Wstart.z()<<",E_x="<<Wend.x()<<",E_y="<<Wend.y()<<",E_z="<<Wend.z()<<std::endl; + + TVector3 denominator = (Wend-Wstart) ; float min_distance = 999 ; for(unsigned int i=0; i< simhit_size; i++) { - dd4hep::rec::Vector3D west(0,0,0); - dd4hep::rec::Vector3D east(0,0,0); - dd4hep::rec::Vector3D pos(iter->second.at(i).getPosition()[0], iter->second.at(i).getPosition()[1], iter->second.at(i).getPosition()[2]); - dd4hep::rec::Vector3D numerator = (east-west).cross(west-pos) ; - dd4hep::rec::Vector3D denominator = (east-west) ; - float distance = numerator.r()/denominator.r() ; - std::cout<<"distance="<<distance<<std::endl; - if(distance < min_distance){ - min_distance = distance; + float sim_hit_mom = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] + iter->second.at(i).getMomentum()[2]*iter->second.at(i).getMomentum()[2] );//GeV + float sim_hit_pt = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] );//GeV + TVector3 pos(iter->second.at(i).getPosition()[0], iter->second.at(i).getPosition()[1], iter->second.at(i).getPosition()[2]); + TVector3 numerator = denominator.Cross(Wstart-pos) ; + float tmp_distance = numerator.Mag()/denominator.Mag() ; + //std::cout<<"tmp_distance="<<tmp_distance<<",x="<<pos.x()<<",y="<<pos.y()<<",z="<<pos.z()<<",mom="<<sim_hit_mom<<",pt="<<sim_hit_pt<<std::endl; + if(tmp_distance < min_distance){ + min_distance = tmp_distance; pos_x = pos.x(); pos_y = pos.y(); pos_z = pos.z(); } tot_length += iter->second.at(i).getPathLength();//mm - /* - tot_x += iter->second.at(i).getEDep()*iter->second.at(i).getPosition()[0]; - tot_y += iter->second.at(i).getEDep()*iter->second.at(i).getPosition()[1]; - tot_z += iter->second.at(i).getEDep()*iter->second.at(i).getPosition()[2]; - */ auto asso = AssoVec->create(); asso.setRec(trkHit); asso.setSim(iter->second.at(i)); asso.setWeight(iter->second.at(i).getEDep()/tot_edep); + + if(m_WriteAna){ + m_simhit_x[m_n_sim] = pos.x(); + m_simhit_y[m_n_sim] = pos.y(); + m_simhit_z[m_n_sim] = pos.z(); + m_n_sim ++ ; + } } trkHit.setTime(min_distance*1e3/m_velocity);//m_velocity is um/ns, drift time in ns trkHit.setEDep(tot_edep);// GeV - trkHit.setEdx (tot_edep*1000/(tot_length/10) ); // MeV/cm, need check! - //trkHit.setPosition (edm4hep::Vector3d(tot_x/tot_edep, tot_y/tot_edep, tot_z/tot_edep));//center mass + trkHit.setEdx (tot_edep/tot_length); // GeV/mm trkHit.setPosition (edm4hep::Vector3d(pos_x, pos_y, pos_z));//position of closest sim hit trkHit.setCovMatrix(std::array<float, 6>{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm + + if(m_WriteAna){ + m_chamber [m_n_digi] = chamber; + m_layer [m_n_digi] = layer ; + m_cell [m_n_digi] = cellID; + m_cell_x [m_n_digi] = Wstart.x(); + m_cell_y [m_n_digi] = Wstart.y(); + m_hit_x [m_n_digi] = pos_x; + m_hit_y [m_n_digi] = pos_y; + m_hit_z [m_n_digi] = pos_z; + m_dca [m_n_digi] = min_distance; + m_hit_dE [m_n_digi] = trkHit.getEDep(); + m_hit_dE_dx[m_n_digi] = trkHit.getEdx() ; + m_n_digi ++ ; + } } + + std::cout<<"output digi DCHhit size="<< Vec->size() <<std::endl; _nEvt ++ ; + if(m_WriteAna){ + StatusCode status = m_tuple->write(); + if ( status.isFailure() ) { + error() << " Cannot fill N-tuple:" << long( m_tuple ) << endmsg; + return StatusCode::FAILURE; + } + } return StatusCode::SUCCESS; } diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.h b/Digitisers/DCHDigi/src/DCHDigiAlg.h index 6a6bb320675aae225b5589f2a03599b643f8e120..339f85baaa48970b241ae215b3256bc55955794e 100644 --- a/Digitisers/DCHDigi/src/DCHDigiAlg.h +++ b/Digitisers/DCHDigi/src/DCHDigiAlg.h @@ -2,6 +2,7 @@ #define DCH_DIGI_ALG_H #include "FWCore/DataHandle.h" +#include "GaudiKernel/NTuple.h" #include "GaudiAlg/GaudiAlgorithm.h" #include "edm4hep/SimTrackerHitCollection.h" #include "edm4hep/TrackerHitCollection.h" @@ -10,7 +11,10 @@ #include <DDRec/DetectorData.h> #include <DDRec/CellIDPositionConverter.h> #include "DetInterface/IGeomSvc.h" +#include "DDSegmentation/Segmentation.h" +#include "DetSegmentation/GridDriftChamber.h" +#include "TVector3.h" @@ -40,13 +44,39 @@ protected: typedef std::vector<float> FloatVec; int _nEvt ; - //float m_length; + NTuple::Tuple* m_tuple = nullptr ; + NTuple::Item<long> m_n_sim; + NTuple::Item<long> m_n_digi; + NTuple::Array<int > m_chamber ; + NTuple::Array<int > m_layer ; + NTuple::Array<int > m_cell ; + NTuple::Array<float> m_cell_x ; + NTuple::Array<float> m_cell_y ; + NTuple::Array<float> m_simhit_x ; + NTuple::Array<float> m_simhit_y ; + NTuple::Array<float> m_simhit_z ; + NTuple::Array<float> m_hit_x ; + NTuple::Array<float> m_hit_y ; + NTuple::Array<float> m_hit_z ; + NTuple::Array<float> m_dca ; + NTuple::Array<float> m_hit_dE ; + NTuple::Array<float> m_hit_dE_dx ; + + + dd4hep::rec::CellIDPositionConverter* m_cellIDConverter; + dd4hep::DDSegmentation::GridDriftChamber* m_segmentation; + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + + Gaudi::Property<std::string> m_readout_name{ this, "readout", "DriftChamberHitsCollection"};//readout for getting segmentation Gaudi::Property<float> m_res_x { this, "res_x", 0.11};//mm Gaudi::Property<float> m_res_y { this, "res_y", 0.11};//mm Gaudi::Property<float> m_res_z { this, "res_z", 1 };//mm Gaudi::Property<float> m_velocity { this, "drift_velocity", 40};// um/ns + Gaudi::Property<float> m_mom_threshold { this, "mom_threshold", 0};// GeV + Gaudi::Property<bool> m_WriteAna { this, "WriteAna", false}; + // Input collections DataHandle<edm4hep::SimTrackerHitCollection> r_SimDCHCol{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this}; diff --git a/Examples/options/tut_detsim_digi_SDT.py b/Examples/options/tut_detsim_digi_SDT.py new file mode 100644 index 0000000000000000000000000000000000000000..6511cf25954101ca5322d6b727e08eea4a63bdab --- /dev/null +++ b/Examples/options/tut_detsim_digi_SDT.py @@ -0,0 +1,189 @@ +#!/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'"] +############################################################################## +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 = "digi_DCH.root" +out.outputCommands = ["keep *"] + +############################################################################## +# ApplicationMgr +############################################################################## + +from Configurables import ApplicationMgr +ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, out], + EvtSel = 'NONE', + EvtMax = 10, + ExtSvc = [rndmengine, dsvc, geosvc], + HistogramPersistency = "ROOT", + OutputLevel=INFO +)