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..e20b6d10f7d8bb5035367959d8d668eac2e6dec7 100644 --- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp +++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp @@ -29,39 +29,83 @@ 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; + } + + std::string s_output=m_Output; + if(m_WriteAna){ + m_fout = new TFile(s_output.c_str(),"RECREATE"); + m_tree = new TTree("evt","tree"); + m_tree->Branch("chamber" , &m_chamber ); + m_tree->Branch("layer" , &m_layer ); + m_tree->Branch("cell" , &m_cell ); + m_tree->Branch("cell_x" , &m_cell_x ); + m_tree->Branch("cell_y" , &m_cell_y ); + m_tree->Branch("simhit_x" , &m_simhit_x ); + m_tree->Branch("simhit_y" , &m_simhit_y ); + m_tree->Branch("simhit_z" , &m_simhit_z ); + m_tree->Branch("hit_x" , &m_hit_x ); + m_tree->Branch("hit_y" , &m_hit_y ); + m_tree->Branch("hit_z" , &m_hit_z ); + m_tree->Branch("dca" , &m_dca ); + m_tree->Branch("hit_dE" , &m_hit_dE ); + m_tree->Branch("hit_dE_dx", &m_hit_dE_dx); + } std::cout<<"DCHDigiAlg::initialized"<< std::endl; return GaudiAlgorithm::initialize(); } +void DCHDigiAlg::Reset() +{ + + std::vector<int >().swap( m_chamber ); + std::vector<int >().swap( m_layer ); + std::vector<int >().swap( m_cell ); + std::vector<float>().swap( m_cell_x ); + std::vector<float>().swap( m_cell_y ); + std::vector<float>().swap( m_simhit_x ); + std::vector<float>().swap( m_simhit_y ); + std::vector<float>().swap( m_simhit_z ); + std::vector<float>().swap( m_hit_x ); + std::vector<float>().swap( m_hit_y ); + std::vector<float>().swap( m_hit_z ); + std::vector<float>().swap( m_dca ); + std::vector<float>().swap( m_hit_dE ); + std::vector<float>().swap( m_hit_dE_dx ); + +} + StatusCode DCHDigiAlg::execute() { + info() << "Processing " << _nEvt << " events " << endmsg; + Reset(); 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 @@ -74,11 +118,11 @@ StatusCode DCHDigiAlg::execute() 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,50 +131,81 @@ 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.push_back(pos.x()); + m_simhit_y.push_back(pos.y()); + m_simhit_z.push_back(pos.z()); + } } 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.push_back(chamber); + m_layer .push_back(layer ); + m_cell .push_back(cellID); + m_cell_x .push_back(Wstart.x()); + m_cell_y .push_back(Wstart.y()); + m_hit_x .push_back(pos_x); + m_hit_y .push_back(pos_y); + m_hit_z .push_back(pos_z); + m_dca .push_back(min_distance); + m_hit_dE .push_back(trkHit.getEDep()); + m_hit_dE_dx .push_back(trkHit.getEdx() ); + } } std::cout<<"output digi DCHhit size="<< Vec->size() <<std::endl; _nEvt ++ ; + if(m_WriteAna) m_tree->Fill(); + return StatusCode::SUCCESS; } StatusCode DCHDigiAlg::finalize() { info() << "Processed " << _nEvt << " events " << endmsg; + if(m_WriteAna){ + m_fout->cd(); + m_tree->Write(); + m_fout->Close(); + } return GaudiAlgorithm::finalize(); } diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.h b/Digitisers/DCHDigi/src/DCHDigiAlg.h index 6a6bb320675aae225b5589f2a03599b643f8e120..7e38bf259fe756f9cf280ee99bbbd84cd881864e 100644 --- a/Digitisers/DCHDigi/src/DCHDigiAlg.h +++ b/Digitisers/DCHDigi/src/DCHDigiAlg.h @@ -10,7 +10,13 @@ #include <DDRec/DetectorData.h> #include <DDRec/CellIDPositionConverter.h> #include "DetInterface/IGeomSvc.h" +#include "DDSegmentation/Segmentation.h" +#include "DetSegmentation/GridDriftChamber.h" +#include "TVector3.h" +#include "TROOT.h" +#include "TTree.h" +#include "TFile.h" @@ -33,6 +39,7 @@ public: /** Called after data processing for clean up. */ virtual StatusCode finalize() ; + void Reset(); protected: @@ -40,13 +47,39 @@ protected: typedef std::vector<float> FloatVec; int _nEvt ; - //float m_length; + TFile* m_fout; + TTree* m_tree; + std::vector<int > m_chamber ; + std::vector<int > m_layer ; + std::vector<int > m_cell ; + std::vector<float> m_cell_x ; + std::vector<float> m_cell_y ; + std::vector<float> m_simhit_x ; + std::vector<float> m_simhit_y ; + std::vector<float> m_simhit_z ; + std::vector<float> m_hit_x ; + std::vector<float> m_hit_y ; + std::vector<float> m_hit_z ; + std::vector<float> m_dca ; + std::vector<float> m_hit_dE ; + std::vector<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}; + + Gaudi::Property<std::string> m_Output { this, "output", "ana_DCH_digi.root"}; // 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..a234b179f8591a92f4315f2896bf55115412f80c --- /dev/null +++ b/Examples/options/tut_detsim_digi_SDT.py @@ -0,0 +1,185 @@ +#!/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 = [0] # 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 +dCHDigiAlg.output = "ana_DCH_digi_v1.root" + +############################################################################## +# POD I/O +############################################################################## +from Configurables import PodioOutput +out = PodioOutput("outputalg") +out.filename = "detsim_digi_DCH.root" +out.outputCommands = ["keep *"] + +############################################################################## +# ApplicationMgr +############################################################################## + +from Configurables import ApplicationMgr +#ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, out], +ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg], + EvtSel = 'NONE', + EvtMax = 10, + ExtSvc = [rndmengine, dsvc, geosvc], +)