From 553b3aae7c0052c847493ef8c3dea9515012fddc Mon Sep 17 00:00:00 2001 From: wenxingfang <1473717798@qq.com> Date: Fri, 16 Jun 2023 16:32:22 +0800 Subject: [PATCH 1/6] add TrackHeedSimTool for DC cell simulation --- Examples/options/tut_detsim_SDT_Heed.py | 207 ++++++ .../src/Edm4hepWriterAnaElemTool.cpp | 19 + .../DetSimAna/src/Edm4hepWriterAnaElemTool.h | 3 + Simulation/DetSimDedx/CMakeLists.txt | 13 + .../DetSimDedx/src/TrackHeedSimTool.cpp | 646 ++++++++++++++++++ Simulation/DetSimDedx/src/TrackHeedSimTool.h | 156 +++++ .../include/DetSimInterface/IDedxSimTool.h | 2 + cmake/FindOnnxRuntime.cmake | 35 + 8 files changed, 1081 insertions(+) create mode 100644 Examples/options/tut_detsim_SDT_Heed.py create mode 100644 Simulation/DetSimDedx/src/TrackHeedSimTool.cpp create mode 100644 Simulation/DetSimDedx/src/TrackHeedSimTool.h create mode 100644 cmake/FindOnnxRuntime.cmake diff --git a/Examples/options/tut_detsim_SDT_Heed.py b/Examples/options/tut_detsim_SDT_Heed.py new file mode 100644 index 000000000..049c3c226 --- /dev/null +++ b/Examples/options/tut_detsim_SDT_Heed.py @@ -0,0 +1,207 @@ +#!/usr/bin/env python + +import os +import sys +# sys.exit(0) + +from Gaudi.Configuration import * + +############################################################################## +# Random Number Svc +############################################################################## +from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ + +seed = [42] + +# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi +rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine") # The default engine in Geant4 +rndmengine.SetSingleton = True +rndmengine.Seeds = seed + +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" +#geometry_option = "CepC_v4.xml" +det_root = "DETDRIFTCHAMBERROOT" +#det_root = "DETCEPCV4ROOT"#"DETDRIFTCHAMBERROOT" +if not os.getenv(det_root): + print("Can't find the geometry. Please setup envvar %s."%det_root ) + sys.exit(-1) + +geometry_path = os.path.join(os.getenv(det_root), "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") +print('geometry_path=',geometry_path) +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 = [80] # rad; 45deg +gun.ThetaMaxs = [90.] # rad; 45deg + +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") +detsimalg.RandomSeeds = seed + +if int(os.environ.get("VIS", 0)): + detsimalg.VisMacs = ["vis.mac"] + +detsimalg.RunCmds = [ +# "/tracking/verbose 1", +] + +from Configurables import DummyFastSimG4Tool +dummy_fastsim_tool = DummyFastSimG4Tool("DummyFastSimG4Tool") + +detsimalg.FastSimG4Tools = [ +# "DummyFastSimG4Tool" +] + +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" +dedxoption = "TrackHeedSimTool" + +driftchamber_sensdettool.DedxSimTool = dedxoption + +from Configurables import DummyDedxSimTool +from Configurables import BetheBlochEquationDedxSimTool +from Configurables import TrackHeedSimTool + +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 +elif dedxoption == "TrackHeedSimTool": + dedx_simtool = TrackHeedSimTool("TrackHeedSimTool") + dedx_simtool.only_primary = False#True + dedx_simtool.use_max_step = True#True + dedx_simtool.max_step = 1#mm + #dedx_simtool.he = 50 + #dedx_simtool.isob = 50 + #dedx_simtool.gas_file ="/junofs/users/wxfang/MyGit/tmp/check_G4FastSim_20210121/CEPCSW/Digitisers/DigiGarfield/He_50_isobutane_50.gas" + dedx_simtool.he = 90 + dedx_simtool.isob = 10 + #dedx_simtool.gas_file ="/junofs/users/wxfang/MyGit/tmp/check_G4FastSim_20210121/CEPCSW/Digitisers/DigiGarfield/he_90_isobutane_10.gas" + #dedx_simtool.IonMobility_file ="/junofs/users/wxfang/MyGit/tmp/check_G4FastSim_20210121/CEPCSW/Digitisers/DigiGarfield/IonMobility_He+_He.txt" + dedx_simtool.gas_file ="he_90_isobutane_10.gas" + dedx_simtool.IonMobility_file ="IonMobility_He+_He.txt" + dedx_simtool.save_mc = True##IF this is False then ... + dedx_simtool.debug = False + dedx_simtool.sim_pulse = True + #dedx_simtool.model='/junofs/users/wxfang/MyGit/tmp/fork_cepcsw_20220418/CEPCSW/Digitisers/SimCurrentONNX/src/model_test.onnx' + #dedx_simtool.model='/junofs/users/wxfang/MyGit/tmp/fork_cepcsw_20220418/CEPCSW/Digitisers/SimCurrentONNX/src/model_90He10C4H10_18mm.onnx' + dedx_simtool.model='model_90He10C4H10_18mm.onnx' + dedx_simtool.batchsize = 100 + +############################################################################## +# POD I/O +############################################################################## +from Configurables import PodioOutput +out = PodioOutput("outputalg") +out.filename = "detsim_heed.root" +out.outputCommands = ["keep *"] + +############################################################################## +# ApplicationMgr +############################################################################## + +from Configurables import ApplicationMgr +ApplicationMgr( TopAlg = [genalg, detsimalg, out], + EvtSel = 'NONE', + EvtMax = 20, + ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc], + OutputLevel=INFO +) diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp index cfbf31dfc..8e121bac3 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp @@ -14,6 +14,7 @@ #include "DDG4/Geant4HitCollection.h" #include "DDG4/Geant4Data.h" #include "DetSimSD/Geant4Hits.h" +#include DECLARE_COMPONENT(Edm4hepWriterAnaElemTool) @@ -42,6 +43,7 @@ Edm4hepWriterAnaElemTool::BeginOfRunAction(const G4Run*) { } else { error() << "Failed to find GeomSvc." << endmsg; } + } void @@ -87,6 +89,14 @@ Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) { // reset m_track2primary.clear(); + auto SimPIonCol = m_SimPrimaryIonizationCol.createAndPut(); + ToolHandleArray tmp_m_dedx_tools; + tmp_m_dedx_tools.push_back("TrackHeedSimTool"); + for (auto dedxtool: tmp_m_dedx_tools) { + debug() << "reset dedx_tool:" <reset(); + } + } void @@ -94,6 +104,15 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) { msg() << "mcCol size (after simulation) : " << mcCol->size() << endmsg; // save all data + auto SimPrimaryIonizationCol = m_SimPrimaryIonizationCol.get(); + msg() << "SimPrimaryIonizationCol size ="<size()< tmp_m_dedx_tools; + tmp_m_dedx_tools.push_back("TrackHeedSimTool"); + for (auto dedxtool: tmp_m_dedx_tools) { + debug() << "call endOfEvent() for dedx_tool:" << dedxtool << endmsg; + dedxtool->endOfEvent(); + } + // create collections. auto trackercols = m_trackerCol.createAndPut(); diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h index e4415b178..60cc930a4 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h @@ -15,6 +15,7 @@ #include "edm4hep/SimTrackerHitCollection.h" #include "edm4hep/SimCalorimeterHitCollection.h" #include "edm4hep/CaloHitContributionCollection.h" +#include "edm4hep/SimPrimaryIonizationClusterCollection.h" class Edm4hepWriterAnaElemTool: public extends { @@ -129,6 +130,8 @@ class Edm4hepWriterAnaElemTool: public extends { "DriftChamberHitsCollection", Gaudi::DataHandle::Writer, this}; + // for ionized electron + DataHandle m_SimPrimaryIonizationCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Writer, this}; private: // in order to associate the hit contribution with the primary track, diff --git a/Simulation/DetSimDedx/CMakeLists.txt b/Simulation/DetSimDedx/CMakeLists.txt index fb8e4d3a6..03bf40d9d 100644 --- a/Simulation/DetSimDedx/CMakeLists.txt +++ b/Simulation/DetSimDedx/CMakeLists.txt @@ -1,15 +1,28 @@ find_package(Geant4 REQUIRED ui_all vis_all) include(${Geant4_USE_FILE}) +find_package(Garfield REQUIRED) +message("libonnxruntime ${OnnxRuntime_LIBRARY}") +find_package(OnnxRuntime REQUIRED) +message("libonnxruntime ${OnnxRuntime_LIBRARY}") gaudi_add_module(DetSimDedx SOURCES src/DummyDedxSimTool.cpp src/BetheBlochEquationDedxSimTool.cpp src/GFDndxSimTool.cpp + src/TrackHeedSimTool.cpp LINK DetSimInterface + DetInterface + DetSegmentation ${DD4hep_COMPONENT_LIBRARIES} Gaudi::GaudiKernel EDM4HEP::edm4hep EDM4HEP::edm4hepDict + k4FWCore::k4FWCore + Garfield::Garfield + ${OnnxRuntime_LIBRARY} + #/cvmfs/sft.cern.ch/lcg/views/LCG_103/x86_64-centos7-gcc11-opt/lib/libonnxruntime.so + ${CLHEP_LIBRARIES} + ) diff --git a/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp b/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp new file mode 100644 index 000000000..fb03cefec --- /dev/null +++ b/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp @@ -0,0 +1,646 @@ +#include "TrackHeedSimTool.h" +#include "G4Step.hh" +#include "G4SystemOfUnits.hh" +#include +#include "G4TransportationManager.hh" +#include "G4Navigator.hh" +#include "G4VPhysicalVolume.hh" +#include "G4LogicalVolume.hh" +#include +#include "DDG4/Geant4Converter.h" +#include "DetSegmentation/GridDriftChamber.h" + + +#include +#include +#include +#include +#include "CLHEP/Random/RandGauss.h" + + +DECLARE_COMPONENT(TrackHeedSimTool) + +double TrackHeedSimTool::dedx(const edm4hep::MCParticle& mc) {return 0;} +double TrackHeedSimTool::dndx(double betagamma) {return 0;} + +double TrackHeedSimTool::dedx(const G4Step* Step) +{ + + clock_t t0 = clock(); + double de = 0; + float cm_to_mm = 10; + G4Step* aStep = const_cast(Step); + G4Track* g4Track = aStep->GetTrack(); + int pdg_code = g4Track->GetParticleDefinition()->GetPDGEncoding(); + G4double pdg_mass = g4Track->GetParticleDefinition()->GetPDGMass(); + G4double pdg_charge = g4Track->GetParticleDefinition()->GetPDGCharge(); + const G4VProcess* creatorProcess = g4Track->GetCreatorProcess(); + const G4String tmp_str_pro = (creatorProcess !=0) ? creatorProcess->GetProcessName() : "normal"; + G4double gammabeta=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma(); + if(g4Track->GetParticleDefinition()->GetPDGCharge() ==0) return 0;//skip neutral particle + if(gammabeta<0.01)return m_eps;//too low momentum + if(m_only_primary.value() && g4Track->GetParentID() != 0) return m_eps; + if(g4Track->GetKineticEnergy() <=0) return 0; + if(pdg_code == 11 && (tmp_str_pro=="phot" || tmp_str_pro=="hIoni" || tmp_str_pro=="eIoni" || tmp_str_pro=="muIoni" || tmp_str_pro=="ionIoni" ) ) return m_eps;//skip the electron produced by Ioni, because it is already simulated by TrackHeed + if(m_particle_map.find(pdg_code) == m_particle_map.end() ) return m_eps; + edm4hep::SimPrimaryIonizationClusterCollection* SimPrimaryIonizationCol = nullptr; + edm4hep::MCParticleCollection* mcCol = nullptr; + try{ + SimPrimaryIonizationCol = const_cast(m_SimPrimaryIonizationCol.get()); + mcCol = const_cast(m_mc_handle.get()); + } + catch(...){ + G4cout<<"Error! Can't find collection in event, please check it have been createAndPut() in Begin of event"<GetPreStepPoint()->GetKineticEnergy(); + G4double track_time = aStep->GetPreStepPoint()->GetGlobalTime(); + G4double track_dx = aStep->GetPreStepPoint()->GetMomentumDirection ().x(); + G4double track_dy = aStep->GetPreStepPoint()->GetMomentumDirection ().y(); + G4double track_dz = aStep->GetPreStepPoint()->GetMomentumDirection ().z(); + G4double track_length = aStep->GetStepLength(); + G4double position_x = aStep->GetPreStepPoint()->GetPosition().x(); + G4double position_y = aStep->GetPreStepPoint()->GetPosition().y(); + G4double position_z = aStep->GetPreStepPoint()->GetPosition().z(); + int track_ID = g4Track->GetTrackID(); + int Parent_ID = g4Track->GetParentID(); + bool update_ke = true; + if(m_use_max_step.value()){ + bool do_sim = false; + if(m_isFirst){ + m_pre_x = aStep->GetPreStepPoint()->GetPosition().x(); + m_pre_y = aStep->GetPreStepPoint()->GetPosition().y(); + m_pre_z = aStep->GetPreStepPoint()->GetPosition().z(); + m_pre_dx = aStep->GetPreStepPoint()->GetMomentumDirection().x(); + m_pre_dy = aStep->GetPreStepPoint()->GetMomentumDirection().y(); + m_pre_dz = aStep->GetPreStepPoint()->GetMomentumDirection().z(); + m_pre_t = aStep->GetPreStepPoint()->GetGlobalTime(); + m_post_point = aStep->GetPostStepPoint(); + m_total_range += track_length; + m_current_track_ID = g4Track->GetTrackID(); + m_current_Parent_ID = g4Track->GetParentID(); + m_pdg_code = g4Track->GetParticleDefinition()->GetPDGEncoding(); + m_isFirst = false; + m_pa_KE = aStep->GetPreStepPoint()->GetKineticEnergy(); + } + else{ + + if(g4Track->GetTrackID() != m_current_track_ID){ + do_sim = true; + m_change_track = true; + update_ke = false; + } + else{ + m_post_point = aStep->GetPostStepPoint(); + m_total_range += track_length; + } + } + if(m_total_range/CLHEP::mm >= m_max_step.value()){ + do_sim = true; + } + if(do_sim){ + track_KE = m_pa_KE; + pdg_code = m_pdg_code; + track_length = m_total_range; + track_time = m_pre_t; + track_dx = m_pre_dx; + track_dy = m_pre_dy; + track_dz = m_pre_dz; + position_x = m_pre_x; + position_y = m_pre_y; + position_z = m_pre_z; + track_ID = m_current_track_ID; + Parent_ID = m_current_Parent_ID; + if(m_change_track){ + m_pre_x = aStep->GetPreStepPoint()->GetPosition().x(); + m_pre_y = aStep->GetPreStepPoint()->GetPosition().y(); + m_pre_z = aStep->GetPreStepPoint()->GetPosition().z(); + m_pre_dx = aStep->GetPreStepPoint()->GetMomentumDirection().x(); + m_pre_dy = aStep->GetPreStepPoint()->GetMomentumDirection().y(); + m_pre_dz = aStep->GetPreStepPoint()->GetMomentumDirection().z(); + m_pre_t = aStep->GetPreStepPoint()->GetGlobalTime(); + m_post_point = aStep->GetPostStepPoint(); + m_total_range = aStep->GetStepLength(); + m_current_track_ID = g4Track->GetTrackID(); + m_current_Parent_ID = g4Track->GetParentID(); + m_pdg_code = g4Track->GetParticleDefinition()->GetPDGEncoding(); + m_change_track = false; + } + else{ + m_pre_x = aStep->GetPostStepPoint()->GetPosition().x(); + m_pre_y = aStep->GetPostStepPoint()->GetPosition().y(); + m_pre_z = aStep->GetPostStepPoint()->GetPosition().z(); + m_pre_dx = aStep->GetPostStepPoint()->GetMomentumDirection().x(); + m_pre_dy = aStep->GetPostStepPoint()->GetMomentumDirection().y(); + m_pre_dz = aStep->GetPostStepPoint()->GetMomentumDirection().z(); + m_pre_t = aStep->GetPostStepPoint()->GetGlobalTime(); + m_total_range = 0; + } + } + else return m_eps; + } + + float init_x = 10;//cm + float init_y = -10;//cm + /* + if(pdg_code == 11 && track_KE/CLHEP::keV < m_delta_threshold.value()){ + int nc = 0, ni=0; + m_track->TransportDeltaElectron(init_x, init_y, 0, track_time/CLHEP::ns, track_KE/CLHEP::eV, track_dx, track_dy, track_dz, nc, ni); + for (int j = 0; j < nc; ++j) { + double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.; + double dx = 0., dy = 0., dz = 0.; + m_track->GetElectron(j, xe, ye, ze, te, ee, dx, dy, dz); + auto ehit = SimIonizationCol->create(); + ehit.setTime(te); + double epos[3] = { cm_to_mm*( (xe - init_x)+position_x/CLHEP::cm) , cm_to_mm*((ye - init_y)+position_y/CLHEP::cm), cm_to_mm*(ze + position_z/CLHEP::cm)}; + ehit.setPosition(edm4hep::Vector3d(epos)); + ehit.setType(11); + } + g4Track->SetTrackStatus(fStopAndKill); + return 0; + } + */ + clock_t t01 = clock(); + //cmp.SetMagneticField(0., 0., -3.); + m_track->SetParticle(m_particle_map[pdg_code]); + + bool change_KE = false; + if( abs(m_previous_KE -(track_KE/CLHEP::eV) )/(track_KE/CLHEP::eV) > m_change_threshold.value()) { + change_KE = true; + m_previous_KE = track_KE/CLHEP::eV; + } + + bool change_ID = false; + if(m_previous_track_ID != track_ID ){ + change_ID = true; + m_previous_track_ID = track_ID; + m_previous_KE = track_KE/CLHEP::eV; + } + if(change_ID || change_KE ){ + m_track->SetKineticEnergy(track_KE/CLHEP::eV); + } + + + m_track->EnableOneStepFly(true); + m_track->SetSteppingLimits( track_length/CLHEP::cm, 1000, 0.1, 0.2); + clock_t t012 = clock(); + m_track->NewTrack(init_x, init_y, 0, track_time/CLHEP::ns, track_dx, track_dy, track_dz);//cm + double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.; + int nc = 0; + int ic = 0; + int first=true; + clock_t t02 = clock(); + while (m_track->GetCluster(xc, yc, zc, tc, nc, ec, extra)) { + //auto chit = SimHitCol->create(); + auto chit = SimPrimaryIonizationCol->create(); + chit.setTime(tc); + double cpos[3] = { cm_to_mm*( (xc - init_x)+position_x/CLHEP::cm) , cm_to_mm*((yc - init_y)+position_y/CLHEP::cm), cm_to_mm*(zc + position_z/CLHEP::cm)}; + chit.setPosition(edm4hep::Vector3d(cpos)); + //float cmom[3] = {0,0,0}; + //getMom(ec, 1, 0, 0, cmom);//FIXME direction is not important? + chit.setType(0);//default + if(m_save_cellID) chit.setCellID( getCellID(cpos[0], cpos[1], cpos[2]) ); + if(m_save_mc && Parent_ID == 0 && track_ID <= mcCol->size() && mcCol ){ + chit.setMCParticle( mcCol->at(track_ID-1) ); + //std::cout<<"mc obj index="<at(track_ID-1).getObjectID().index<lcdd(); + if ( !m_dd4hep ) throw "TrackHeedSimTool :Failed to get dd4hep::Detector ..."; + m_readout = new dd4hep::Readout( m_dd4hep->readout(m_readout_name) ); + if ( !m_readout ) throw "TrackHeedSimTool :Failed to get readout ..."; + m_segmentation = dynamic_cast(m_readout->segmentation().segmentation()); + if ( !m_segmentation ) throw "TrackHeedSimTool :Failed to get segmentation ..."; + + m_particle_map[ 11] = "e-"; + m_particle_map[-11] = "e+"; + m_particle_map[ 13] = "mu-"; + m_particle_map[-13] = "mu+"; + m_particle_map[ 211] = "pi+"; + m_particle_map[-211] = "pi-"; + m_particle_map[ 321] = "K+"; + m_particle_map[-321] = "K-"; + m_particle_map[2212] = "p"; + m_particle_map[-2212] = "pbar"; + m_particle_map[700201] = "d"; + m_particle_map[700202] = "alpha"; + + m_gas.SetComposition("he", m_he,"isobutane", m_isob); + m_gas.SetTemperature(293.15); + m_gas.SetPressure(760.0); + m_gas.SetMaxElectronEnergy(200.); + m_gas.EnablePenningTransfer(0.44, 0.0, "He"); + m_gas.LoadGasFile(m_gas_file.value()); + m_gas.LoadIonMobility(m_IonMobility.value()); + //std::this_thread::sleep_for(std::chrono::milliseconds(m_delay_time)); + //m_gas.LoadGasFile("/junofs/users/wxfang/MyGit/tmp/check_G4FastSim_20210121/CEPCSW/Digitisers/DigiGarfield/He_50_isobutane_50.gas"); + //m_gas.LoadIonMobility("/junofs/users/wxfang/MyGit/tmp/check_G4FastSim_20210121/CEPCSW/Digitisers/DigiGarfield/IonMobility_He+_He.txt"); + /* + m_gas.SetComposition("he", 90.,"isobutane", 10.); // cepc gas + m_gas.SetPressure(760.0); + m_gas.SetTemperature(293.15); + m_gas.SetFieldGrid(100., 100000., 20, true); + m_gas.GenerateGasTable(10); + m_gas.WriteGasFile("he_90_isobutane_10.gas"); + */ + + cmp.SetMedium(&m_gas); + // Field Wire radius [cm] + const double rFWire = 110.e-4; + // Signa Wire radius [cm] + const double rSWire = 25.e-4; + // Cell radius [cm] + float rCell = 50;//As the ionization process is almost not effected by cell geometry and wire voltage. Here the radius is to make sure the ionization process is completed. + // Voltages + const double vSWire = 2000.; + const double vFWire = 0.; + // Add the signal wire in the centre. + cmp.AddWire(0, 0, 2 * rSWire, vSWire, "s"); + // Add the field wire around the signal wire. + cmp.AddWire(-rCell, -rCell, 2 * rFWire, vFWire, "f"); + cmp.AddWire( 0., -rCell, 2 * rFWire, vFWire, "f"); + cmp.AddWire( rCell, -rCell, 2 * rFWire, vFWire, "f"); + cmp.AddWire(-rCell, 0., 2 * rFWire, vFWire, "f"); + cmp.AddWire( rCell, 0., 2 * rFWire, vFWire, "f"); + cmp.AddWire(-rCell, rCell, 2 * rFWire, vFWire, "f"); + cmp.AddWire( 0., rCell, 2 * rFWire, vFWire, "f"); + cmp.AddWire( rCell, rCell, 2 * rFWire, vFWire, "f"); + if(m_BField !=0 ) cmp.SetMagneticField(0., 0., m_BField); + cmp.AddReadout("s"); + + + /// + /// Make a sensor. + /// + m_sensor = new Sensor(); + m_sensor->AddComponent(&cmp); + m_sensor->AddElectrode(&cmp, "s"); + // Set the signal time window. [ns] + const double tstep = 0.5; + const double tmin = -0.5 * 0.5; + const unsigned int nbins = 1000; + m_sensor->SetTimeWindow(tmin, tstep, nbins); + m_sensor->ClearSignal(); + + + + m_track = new Garfield::TrackHeed(); + //track->EnableDebugging(); + m_track->SetSensor(m_sensor); + m_track->EnableDeltaElectronTransport(); + //track->DisableDeltaElectronTransport(); + m_track->EnableMagneticField(); + m_track->EnableElectricField();//almost no effect here + + m_current_track_ID = 0; + m_previous_track_ID =0; + m_previous_KE = 0; + m_current_Parent_ID = -1; + m_change_track = false; + m_total_range = 0; + m_isFirst = false; + m_tot_edep = 0; + m_tot_length = 0; + m_pa_KE =0; + m_pdg_code = 0; + m_pre_x = 0; + m_pre_y = 0; + m_pre_z = 0; + m_pre_dx = 0; + m_pre_dy = 0; + m_pre_dz = 0; + m_pre_t = 0; + + + // for NN pulse simulation// + m_env = std::make_shared(ORT_LOGGING_LEVEL_WARNING, "ENV"); + m_seesion_options = std::make_shared(); + m_seesion_options->SetIntraOpNumThreads(m_intra_op_nthreads); + m_seesion_options->SetInterOpNumThreads(m_inter_op_nthreads); + if(m_debug) std::cout << "before load model " << m_model_file.value() << std::endl; + m_session = std::make_shared(*m_env, m_model_file.value().c_str(), *m_seesion_options); + if(m_debug) std::cout << "after load model " << m_model_file.value() << std::endl; + // lambda function to print the dims. + auto dims_str = [&](const auto& dims) { + return std::accumulate(dims.begin(), dims.end(), std::to_string(dims[0]), + [](const std::string& a, int64_t b){ + return a + "x" + std::to_string(b); + }); + }; + // prepare the input + auto num_input_nodes = m_session->GetInputCount(); + if(m_debug) std::cout << "num_input_nodes: " << num_input_nodes << std::endl; + for (size_t i = 0; i < num_input_nodes; ++i) { + auto name = m_session->GetInputNameAllocated(i, m_allocator); + m_inputNodeNameAllocatedStrings.push_back(std::move(name)); + m_input_node_names.push_back(m_inputNodeNameAllocatedStrings.back().get()); + + Ort::TypeInfo type_info = m_session->GetInputTypeInfo(i); + auto tensor_info = type_info.GetTensorTypeAndShapeInfo(); + auto dims = tensor_info.GetShape(); + //dims[0] = 1; //wxfang, FIXME, if it is -1 (dynamic axis), need overwrite it manually + dims[0] = 10; //wxfang, FIXME, if it is -1 (dynamic axis), need overwrite it manually + m_input_node_dims.push_back(dims); + + + if(m_debug) std::cout<< "[" << i << "]" + << " input_name: " << m_inputNodeNameAllocatedStrings.back().get() + << " ndims: " << dims.size() + << " dims: " << dims_str(dims) + << std::endl; + } + // prepare the output + size_t num_output_nodes = m_session->GetOutputCount(); + for(std::size_t i = 0; i < num_output_nodes; i++) { + auto output_name = m_session->GetOutputNameAllocated(i, m_allocator); + m_outputNodeNameAllocatedStrings.push_back(std::move(output_name)); + m_output_node_names.push_back(m_outputNodeNameAllocatedStrings.back().get()); + Ort::TypeInfo type_info = m_session->GetOutputTypeInfo(i); + auto tensor_info = type_info.GetTensorTypeAndShapeInfo(); + ONNXTensorElementDataType type = tensor_info.GetElementType(); + m_output_node_dims = tensor_info.GetShape(); + if(m_debug) std::cout << "[" << i << "]" + << " output_name: " << m_outputNodeNameAllocatedStrings.back().get() + << " ndims: " << m_output_node_dims.size() + << " dims: " << dims_str(m_output_node_dims) + << std::endl; + + } + + return StatusCode::SUCCESS; +} + +void TrackHeedSimTool::wire_xy(float x1, float y1, float z1, float x2, float y2, float z2, float z, float &x, float &y){ + //linear function: + //(x-x1)/(x2-x1)=(y-y1)/(y2-y1)=(z-z1)/(z2-z1) + x = x1+(x2-x1)*(z-z1)/(z2-z1); + y = y1+(y2-y1)*(z-z1)/(z2-z1); +} + +float TrackHeedSimTool::xy2phi(float x, float y){ + float phi = acos(x/sqrt(x*x+y*y)); + if(y < 0) phi = 2*M_PI-phi; + return phi; +} +void TrackHeedSimTool::getLocal(float x1, float y1, float x2, float y2, float& dx, float& dy){ + /* . + . *(x2,y2) + . . + . . + *(x1, y1) + . + . + . + o + */ + float mo1 = sqrt(x1*x1+y1*y1); + float mo2 = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1) ); + float costheta = (x1*(x2-x1)+y1*(y2-y1))/(mo1*mo2); + dy = mo2*costheta; + dx = xy2phi(x2,y2)>xy2phi(x1,y1) ? mo2*sqrt(1-costheta*costheta) : -mo2*sqrt(1-costheta*costheta) ; +} + +float* TrackHeedSimTool::NNPred(std::vector& inputs) +{ + + + std::vector input_tensors; + auto& dims = m_input_node_dims[0]; + //std::cout << "inputs.size()="<Run(Ort::RunOptions{ nullptr }, m_input_node_names.data(), input_tensors.data(), input_tensors.size(), m_output_node_names.data(), m_output_node_names.size()); + //const auto& output_tensor = output_tensors[0]; + auto& output_tensor = output_tensors[0]; + int num_elements = output_tensor.GetTensorTypeAndShapeInfo().GetElementCount(); + //std::cout << "output_tensor num_elements=" << num_elements<< std::endl; + + float* vec2 = new float[num_elements]; + std::memcpy(vec2, output_tensor.GetTensorMutableData(), num_elements * sizeof(float)); + /* + for (int k=0;k +#include "edm4hep/MCParticle.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimPrimaryIonizationClusterCollection.h" +#include "TVector3.h" +#include + +#include "DD4hep/Segmentations.h" +#include "DD4hep/Printout.h" +#include "DD4hep/Detector.h" +#include "DetInterface/IGeomSvc.h" +#include "DetSegmentation/GridDriftChamber.h" + + +#include "Garfield/ViewCell.hh" +#include "Garfield/ViewDrift.hh" +#include "Garfield/ViewSignal.hh" +#include "Garfield/ViewMedium.hh" +#include "Garfield/ComponentAnalyticField.hh" +#include "Garfield/MediumMagboltz.hh" +#include "Garfield/Sensor.hh" +#include "Garfield/DriftLineRKF.hh" +#include "Garfield/AvalancheMicroscopic.hh" +#include "Garfield/AvalancheMC.hh" +#include "Garfield/TrackHeed.hh" +#include "Garfield/ComponentNeBem3d.hh" +#include "Garfield/SolidWire.hh" +#include "Garfield/GeometrySimple.hh" +#include "Garfield/MediumConductor.hh" +#include "Garfield/ViewField.hh" + +#include +#include + +#include "core/session/onnxruntime_cxx_api.h" +using namespace Garfield; + +class TrackHeedSimTool: public extends { + public: + using extends::extends; + + StatusCode initialize() override; + StatusCode finalize() override; + double dedx(const G4Step* Step) override; + double dedx(const edm4hep::MCParticle& mc) override; + double dndx(double betagamma) override; + void getMom(float ee, float dx, float dy,float dz, float mom[3] ); + void reset(){ + m_isFirst = true; + m_previous_track_ID = 0; + m_previous_KE = 0; + m_tot_edep = 0; + //std::cout<<"m_tot_length="<& inputs); + float xy2phi(float x, float y); + void getLocal(float x1, float y1, float x2, float y2, float& dx, float& dy); + private: + //ServiceHandle m_eds; + SmartIF m_geosvc; + dd4hep::Detector* m_dd4hep; + dd4hep::Readout* m_readout; + dd4hep::DDSegmentation::GridDriftChamber* m_segmentation; + Gaudi::Property m_readout_name{ this, "readout", "DriftChamberHitsCollection"};//readout for getting segmentation + Gaudi::Property m_gas_file{ this, "gas_file", "He_50_isobutane_50.gas"};//gas + Gaudi::Property m_IonMobility{ this, "IonMobility_file", "IonMobility_He+_He.txt"}; + Gaudi::Property m_isob {this, "isob", 50, ""}; + Gaudi::Property m_he {this, "he", 50, ""}; + Gaudi::Property m_debug{this, "debug", false}; + Gaudi::Property m_use_max_step{this, "use_max_step", false}; + Gaudi::Property m_update_KE{this, "update_KE", true}; + Gaudi::Property m_max_step {this, "max_step", 1};//mm + Gaudi::Property m_only_primary{this, "only_primary", false}; + Gaudi::Property m_save_mc{this, "save_mc", false}; + Gaudi::Property m_save_cellID{this, "save_cellID", true}; + Gaudi::Property m_delta_threshold{this, "delta_threshold", 50};//keV + Gaudi::Property m_change_threshold {this, "change_threshold", 0.05}; + Gaudi::Property m_BField {this, "BField", -3}; + Gaudi::Property m_eps { this, "eps" , 1e-6 };//very small value, it is returned dedx for unsimulated step (may needed for SimTrackerHit) + // Output collections + DataHandle m_SimPrimaryIonizationCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Writer, this}; + // In order to associate MCParticle with contribution, we need to access MC Particle. + DataHandle m_mc_handle{"MCParticle", Gaudi::DataHandle::Writer, this}; + + TrackHeed* m_track; + ComponentNeBem3d m_nebem; + ComponentAnalyticField cmp; + GeometrySimple m_geo; + MediumConductor m_metal; + MediumMagboltz m_gas; + Sensor* m_sensor; + std::map m_particle_map; + + int m_previous_track_ID; + float m_previous_KE; + int m_current_track_ID; + int m_current_Parent_ID; + int m_pdg_code; + G4StepPoint* m_pre_point; + G4StepPoint* m_post_point; + G4double m_total_range; + bool m_isFirst; + bool m_change_track; + edm4hep::MCParticle m_mc_paricle; + float m_tot_edep; + float m_tot_length; + float m_pa_KE; + + G4double m_pre_x ; + G4double m_pre_y ; + G4double m_pre_z ; + G4double m_pre_dx ; + G4double m_pre_dy ; + G4double m_pre_dz ; + G4double m_pre_t ; + + //// sim pulse from NN /// + Gaudi::Property m_intra_op_nthreads{ this, "intraOpNumThreads", 1}; + Gaudi::Property m_inter_op_nthreads{ this, "interOpNumThreads", 1}; + std::shared_ptr m_env; + std::shared_ptr m_seesion_options; + std::shared_ptr m_session; + Ort::AllocatorWithDefaultOptions m_allocator; + std::vector m_input_node_names; + std::vector> m_input_node_dims; + std::vector m_output_node_names; + std::vector m_output_node_dims; + std::vector m_inputNodeNameAllocatedStrings; + std::vector m_outputNodeNameAllocatedStrings; + + + Gaudi::Property m_sim_pulse { this, "sim_pulse" , true }; + Gaudi::Property m_model_file{ this, "model", "model_test.onnx"}; + Gaudi::Property m_batchsize { this, "batchsize", 100}; + Gaudi::Property m_time_scale { this, "time_scale", 99.0}; + Gaudi::Property m_time_shift { this, "time_shift", 166.4}; + Gaudi::Property m_amp_scale { this, "amp_scale" , 1e-2 }; + Gaudi::Property m_amp_shift { this, "amp_shift" , 0 }; + Gaudi::Property m_x_scale { this, "x_scale" , 5. };// in mm + Gaudi::Property m_y_scale { this, "y_scale" , 5. };// in mm + + + +}; + +#endif diff --git a/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h b/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h index 8c804500a..d4e8e9e4a 100644 --- a/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h +++ b/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h @@ -28,6 +28,8 @@ class IDedxSimTool: virtual public IAlgTool { virtual double dedx(const G4Step* aStep) = 0; virtual double dedx(const edm4hep::MCParticle& mc) = 0; virtual double dndx(double betagamma) = 0; + virtual void reset() {} + virtual void endOfEvent() {} }; diff --git a/cmake/FindOnnxRuntime.cmake b/cmake/FindOnnxRuntime.cmake new file mode 100644 index 000000000..5a9db01ce --- /dev/null +++ b/cmake/FindOnnxRuntime.cmake @@ -0,0 +1,35 @@ +# Find the ONNX Runtime include directory and library. +# +# This module defines the `onnxruntime` imported target that encodes all +# necessary information in its target properties. + +find_library( + OnnxRuntime_LIBRARY + NAMES onnxruntime + PATH_SUFFIXES lib lib32 lib64 + DOC "The ONNXRuntime library") + +if(NOT OnnxRuntime_LIBRARY) + message(FATAL_ERROR "onnxruntime library not found") +endif() + +find_path( + OnnxRuntime_INCLUDE_DIR + NAMES core/session/onnxruntime_cxx_api.h + PATH_SUFFIXES include include/onnxruntime + DOC "The ONNXRuntime include directory") + +if(NOT OnnxRuntime_INCLUDE_DIR) + message(FATAL_ERROR "onnxruntime includes not found") +endif() + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args( + OnnxRuntime + REQUIRED_VARS OnnxRuntime_LIBRARY OnnxRuntime_INCLUDE_DIR) + +add_library(OnnxRuntime SHARED IMPORTED) +set_property(TARGET OnnxRuntime PROPERTY IMPORTED_LOCATION ${OnnxRuntime_LIBRARY}) +set_property(TARGET OnnxRuntime PROPERTY INTERFACE_INCLUDE_DIRECTORIES ${OnnxRuntime_INCLUDE_DIR}) + +mark_as_advanced(OnnxRuntime_FOUND OnnxRuntime_INCLUDE_DIR OnnxRuntime_LIBRARY) From f50cd1af49cd3069580d4e076f0a5badaee23875 Mon Sep 17 00:00:00 2001 From: wenxingfang <1473717798@qq.com> Date: Wed, 8 Nov 2023 10:10:30 +0800 Subject: [PATCH 2/6] update TrackHeedSimTool --- Detector/DetDriftChamber/compact/det.xml | 2 ++ Simulation/DetSimDedx/CMakeLists.txt | 5 ++++- .../DetSimDedx/src/TrackHeedSimTool.cpp | 20 +++++++++++++++++++ Simulation/DetSimDedx/src/TrackHeedSimTool.h | 7 ++++++- 4 files changed, 32 insertions(+), 2 deletions(-) diff --git a/Detector/DetDriftChamber/compact/det.xml b/Detector/DetDriftChamber/compact/det.xml index b62dd532d..fb4bffdee 100644 --- a/Detector/DetDriftChamber/compact/det.xml +++ b/Detector/DetDriftChamber/compact/det.xml @@ -19,6 +19,8 @@ + + diff --git a/Simulation/DetSimDedx/CMakeLists.txt b/Simulation/DetSimDedx/CMakeLists.txt index 03bf40d9d..867ff20cd 100644 --- a/Simulation/DetSimDedx/CMakeLists.txt +++ b/Simulation/DetSimDedx/CMakeLists.txt @@ -2,7 +2,9 @@ find_package(Geant4 REQUIRED ui_all vis_all) include(${Geant4_USE_FILE}) find_package(Garfield REQUIRED) +message(Garfield::Garfield) message("libonnxruntime ${OnnxRuntime_LIBRARY}") +message("libonnxruntime include ${OnnxRuntime_INCLUDE_DIR}") find_package(OnnxRuntime REQUIRED) message("libonnxruntime ${OnnxRuntime_LIBRARY}") @@ -20,7 +22,8 @@ gaudi_add_module(DetSimDedx EDM4HEP::edm4hep EDM4HEP::edm4hepDict k4FWCore::k4FWCore Garfield::Garfield - ${OnnxRuntime_LIBRARY} + OnnxRuntime + #${OnnxRuntime_LIBRARY} #/cvmfs/sft.cern.ch/lcg/views/LCG_103/x86_64-centos7-gcc11-opt/lib/libonnxruntime.so ${CLHEP_LIBRARIES} diff --git a/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp b/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp index fb03cefec..3bafda261 100644 --- a/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp +++ b/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp @@ -428,9 +428,15 @@ StatusCode TrackHeedSimTool::initialize() auto num_input_nodes = m_session->GetInputCount(); if(m_debug) std::cout << "num_input_nodes: " << num_input_nodes << std::endl; for (size_t i = 0; i < num_input_nodes; ++i) { + #if (ORT_API_VERSION >=13) auto name = m_session->GetInputNameAllocated(i, m_allocator); m_inputNodeNameAllocatedStrings.push_back(std::move(name)); m_input_node_names.push_back(m_inputNodeNameAllocatedStrings.back().get()); + #else + auto name = m_session->GetInputName(i, m_allocator); + m_inputNodeNameAllocatedStrings.push_back(name); + m_input_node_names.push_back(m_inputNodeNameAllocatedStrings.back()); + #endif Ort::TypeInfo type_info = m_session->GetInputTypeInfo(i); auto tensor_info = type_info.GetTensorTypeAndShapeInfo(); @@ -441,7 +447,11 @@ StatusCode TrackHeedSimTool::initialize() if(m_debug) std::cout<< "[" << i << "]" + #if (ORT_API_VERSION >=13) << " input_name: " << m_inputNodeNameAllocatedStrings.back().get() + #else + << " input_name: " << m_inputNodeNameAllocatedStrings.back() + #endif << " ndims: " << dims.size() << " dims: " << dims_str(dims) << std::endl; @@ -449,15 +459,25 @@ StatusCode TrackHeedSimTool::initialize() // prepare the output size_t num_output_nodes = m_session->GetOutputCount(); for(std::size_t i = 0; i < num_output_nodes; i++) { + #if (ORT_API_VERSION >=13) auto output_name = m_session->GetOutputNameAllocated(i, m_allocator); m_outputNodeNameAllocatedStrings.push_back(std::move(output_name)); m_output_node_names.push_back(m_outputNodeNameAllocatedStrings.back().get()); + #else + auto output_name = m_session->GetOutputName(i, m_allocator); + m_outputNodeNameAllocatedStrings.push_back(output_name); + m_output_node_names.push_back(m_outputNodeNameAllocatedStrings.back()); + #endif Ort::TypeInfo type_info = m_session->GetOutputTypeInfo(i); auto tensor_info = type_info.GetTensorTypeAndShapeInfo(); ONNXTensorElementDataType type = tensor_info.GetElementType(); m_output_node_dims = tensor_info.GetShape(); if(m_debug) std::cout << "[" << i << "]" + #if (ORT_API_VERSION >=13) << " output_name: " << m_outputNodeNameAllocatedStrings.back().get() + #else + << " output_name: " << m_outputNodeNameAllocatedStrings.back() + #endif << " ndims: " << m_output_node_dims.size() << " dims: " << dims_str(m_output_node_dims) << std::endl; diff --git a/Simulation/DetSimDedx/src/TrackHeedSimTool.h b/Simulation/DetSimDedx/src/TrackHeedSimTool.h index d152fb972..46a9bd89b 100644 --- a/Simulation/DetSimDedx/src/TrackHeedSimTool.h +++ b/Simulation/DetSimDedx/src/TrackHeedSimTool.h @@ -39,6 +39,7 @@ #include #include "core/session/onnxruntime_cxx_api.h" +#include "core/session/onnxruntime_c_api.h" using namespace Garfield; class TrackHeedSimTool: public extends { @@ -135,9 +136,13 @@ class TrackHeedSimTool: public extends { std::vector> m_input_node_dims; std::vector m_output_node_names; std::vector m_output_node_dims; + #if (ORT_API_VERSION >=13) std::vector m_inputNodeNameAllocatedStrings; std::vector m_outputNodeNameAllocatedStrings; - + #else + std::vector m_inputNodeNameAllocatedStrings; + std::vector m_outputNodeNameAllocatedStrings; + #endif Gaudi::Property m_sim_pulse { this, "sim_pulse" , true }; Gaudi::Property m_model_file{ this, "model", "model_test.onnx"}; From 1b760601ef06655a07767f22b58ea609b73d349b Mon Sep 17 00:00:00 2001 From: wenxingfang <1473717798@qq.com> Date: Wed, 8 Nov 2023 11:08:10 +0800 Subject: [PATCH 3/6] update run.sh to 103.0.2 --- run.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/run.sh b/run.sh index 95f14814a..634faa1d6 100755 --- a/run.sh +++ b/run.sh @@ -67,7 +67,7 @@ function run-job() { # The current default platform lcg_platform=x86_64-centos7-gcc11-opt -lcg_version=103.0.0 +lcg_version=103.0.2 bldtool=${CEPCSW_BLDTOOL} # make, ninja # set in env var From 3aac14564740202f462b3e2778e9893da483ffc7 Mon Sep 17 00:00:00 2001 From: wenxingfang <1473717798@qq.com> Date: Thu, 9 Nov 2023 11:39:54 +0800 Subject: [PATCH 4/6] check if TrackHeedSimTool exist in Edm4hepWriterAnaElemTool --- Examples/options/tut_detsim_SDT_Heed.py | 1 - .../src/Edm4hepWriterAnaElemTool.cpp | 25 +++++++++---------- .../DetSimAna/src/Edm4hepWriterAnaElemTool.h | 3 +++ 3 files changed, 15 insertions(+), 14 deletions(-) diff --git a/Examples/options/tut_detsim_SDT_Heed.py b/Examples/options/tut_detsim_SDT_Heed.py index 049c3c226..e3ee257b7 100644 --- a/Examples/options/tut_detsim_SDT_Heed.py +++ b/Examples/options/tut_detsim_SDT_Heed.py @@ -181,7 +181,6 @@ dedx_simtool.save_mc = True##IF this is False then ... dedx_simtool.debug = False dedx_simtool.sim_pulse = True - #dedx_simtool.model='/junofs/users/wxfang/MyGit/tmp/fork_cepcsw_20220418/CEPCSW/Digitisers/SimCurrentONNX/src/model_test.onnx' #dedx_simtool.model='/junofs/users/wxfang/MyGit/tmp/fork_cepcsw_20220418/CEPCSW/Digitisers/SimCurrentONNX/src/model_90He10C4H10_18mm.onnx' dedx_simtool.model='model_90He10C4H10_18mm.onnx' dedx_simtool.batchsize = 100 diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp index 8e121bac3..4d7c3406c 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp @@ -14,7 +14,6 @@ #include "DDG4/Geant4HitCollection.h" #include "DDG4/Geant4Data.h" #include "DetSimSD/Geant4Hits.h" -#include DECLARE_COMPONENT(Edm4hepWriterAnaElemTool) @@ -90,11 +89,9 @@ Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) { m_track2primary.clear(); auto SimPIonCol = m_SimPrimaryIonizationCol.createAndPut(); - ToolHandleArray tmp_m_dedx_tools; - tmp_m_dedx_tools.push_back("TrackHeedSimTool"); - for (auto dedxtool: tmp_m_dedx_tools) { - debug() << "reset dedx_tool:" <reset(); + if(hasTrackHeedSimTool){ + debug() << "reset TrackHeedSimTool" << endmsg; + m_TrackHeedSimTool->reset(); } } @@ -106,14 +103,10 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) { // save all data auto SimPrimaryIonizationCol = m_SimPrimaryIonizationCol.get(); msg() << "SimPrimaryIonizationCol size ="<size()< tmp_m_dedx_tools; - tmp_m_dedx_tools.push_back("TrackHeedSimTool"); - for (auto dedxtool: tmp_m_dedx_tools) { - debug() << "call endOfEvent() for dedx_tool:" << dedxtool << endmsg; - dedxtool->endOfEvent(); + if(hasTrackHeedSimTool){ + debug() << "call endOfEvent() for TrackHeedSimTool" << endmsg; + m_TrackHeedSimTool->endOfEvent(); } - - // create collections. auto trackercols = m_trackerCol.createAndPut(); auto calorimetercols = m_calorimeterCol.createAndPut(); @@ -580,6 +573,12 @@ StatusCode Edm4hepWriterAnaElemTool::initialize() { StatusCode sc; + m_TrackHeedSimTool = ToolHandle("TrackHeedSimTool",nullptr,false); + if(m_TrackHeedSimTool){ + msg() << "find TrackHeedSimTool" << endmsg; + hasTrackHeedSimTool = true; + } + return sc; } diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h index 60cc930a4..6aedbdc60 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h @@ -16,6 +16,7 @@ #include "edm4hep/SimCalorimeterHitCollection.h" #include "edm4hep/CaloHitContributionCollection.h" #include "edm4hep/SimPrimaryIonizationClusterCollection.h" +#include class Edm4hepWriterAnaElemTool: public extends { @@ -155,6 +156,8 @@ class Edm4hepWriterAnaElemTool: public extends { double Z = 0; bool verboseOutput = false; + ToolHandle m_TrackHeedSimTool; + bool hasTrackHeedSimTool = false; }; #endif From 0cbb0430b1e358b4614f9b1da4f6e6eb7e30cf97 Mon Sep 17 00:00:00 2001 From: wenxingfang <1473717798@qq.com> Date: Thu, 9 Nov 2023 14:42:22 +0800 Subject: [PATCH 5/6] save the data in DriftChamberSensitiveDetector --- .../DetSimAna/src/Edm4hepWriterAnaElemTool.cpp | 16 +--------------- .../DetSimAna/src/Edm4hepWriterAnaElemTool.h | 3 --- Simulation/DetSimDedx/src/TrackHeedSimTool.h | 10 +++++----- .../DetSimSD/src/DriftChamberSensDetTool.cpp | 1 - .../src/DriftChamberSensitiveDetector.cpp | 3 ++- 5 files changed, 8 insertions(+), 25 deletions(-) diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp index 4d7c3406c..b6ab38e67 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp @@ -89,10 +89,6 @@ Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) { m_track2primary.clear(); auto SimPIonCol = m_SimPrimaryIonizationCol.createAndPut(); - if(hasTrackHeedSimTool){ - debug() << "reset TrackHeedSimTool" << endmsg; - m_TrackHeedSimTool->reset(); - } } @@ -102,11 +98,7 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) { msg() << "mcCol size (after simulation) : " << mcCol->size() << endmsg; // save all data auto SimPrimaryIonizationCol = m_SimPrimaryIonizationCol.get(); - msg() << "SimPrimaryIonizationCol size ="<size()<endOfEvent(); - } + //msg() << "SimPrimaryIonizationCol size ="<size()<("TrackHeedSimTool",nullptr,false); - if(m_TrackHeedSimTool){ - msg() << "find TrackHeedSimTool" << endmsg; - hasTrackHeedSimTool = true; - } - return sc; } diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h index 6aedbdc60..60cc930a4 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h @@ -16,7 +16,6 @@ #include "edm4hep/SimCalorimeterHitCollection.h" #include "edm4hep/CaloHitContributionCollection.h" #include "edm4hep/SimPrimaryIonizationClusterCollection.h" -#include class Edm4hepWriterAnaElemTool: public extends { @@ -156,8 +155,6 @@ class Edm4hepWriterAnaElemTool: public extends { double Z = 0; bool verboseOutput = false; - ToolHandle m_TrackHeedSimTool; - bool hasTrackHeedSimTool = false; }; #endif diff --git a/Simulation/DetSimDedx/src/TrackHeedSimTool.h b/Simulation/DetSimDedx/src/TrackHeedSimTool.h index 46a9bd89b..8696d4f66 100644 --- a/Simulation/DetSimDedx/src/TrackHeedSimTool.h +++ b/Simulation/DetSimDedx/src/TrackHeedSimTool.h @@ -102,19 +102,19 @@ class TrackHeedSimTool: public extends { Sensor* m_sensor; std::map m_particle_map; - int m_previous_track_ID; - float m_previous_KE; + int m_previous_track_ID=0; + float m_previous_KE=0; int m_current_track_ID; int m_current_Parent_ID; int m_pdg_code; G4StepPoint* m_pre_point; G4StepPoint* m_post_point; G4double m_total_range; - bool m_isFirst; + bool m_isFirst=true; bool m_change_track; edm4hep::MCParticle m_mc_paricle; - float m_tot_edep; - float m_tot_length; + float m_tot_edep=0; + float m_tot_length=0; float m_pa_KE; G4double m_pre_x ; diff --git a/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp b/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp index c911d4707..58c49d6c6 100644 --- a/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp +++ b/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp @@ -22,7 +22,6 @@ StatusCode DriftChamberSensDetTool::initialize() { error() << "Failed to find dedx simtoo." << endmsg; return StatusCode::FAILURE; } - return sc; } diff --git a/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp index c6ff3db08..c1d768021 100644 --- a/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp +++ b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp @@ -71,5 +71,6 @@ DriftChamberSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) { void DriftChamberSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE) { - + m_dedx_simtool->endOfEvent(); + m_dedx_simtool->reset(); } From 2fe6f6ed8c17ffe09c49740b4bfd49b4135d3d14 Mon Sep 17 00:00:00 2001 From: wenxingfang <1473717798@qq.com> Date: Fri, 10 Nov 2023 10:20:25 +0800 Subject: [PATCH 6/6] move SimPrimaryIonizationCol createAndPut to TrackHeedSimTool --- Examples/options/tut_detsim_SDT_Heed.py | 4 +- .../src/Edm4hepWriterAnaElemTool.cpp | 4 -- .../DetSimAna/src/Edm4hepWriterAnaElemTool.h | 3 -- .../DetSimDedx/src/TrackHeedSimTool.cpp | 51 ++++++++----------- Simulation/DetSimDedx/src/TrackHeedSimTool.h | 28 +++++----- .../include/DetSimInterface/IDedxSimTool.h | 1 - .../src/DriftChamberSensitiveDetector.cpp | 1 - 7 files changed, 39 insertions(+), 53 deletions(-) diff --git a/Examples/options/tut_detsim_SDT_Heed.py b/Examples/options/tut_detsim_SDT_Heed.py index e3ee257b7..d994ae614 100644 --- a/Examples/options/tut_detsim_SDT_Heed.py +++ b/Examples/options/tut_detsim_SDT_Heed.py @@ -178,8 +178,8 @@ #dedx_simtool.IonMobility_file ="/junofs/users/wxfang/MyGit/tmp/check_G4FastSim_20210121/CEPCSW/Digitisers/DigiGarfield/IonMobility_He+_He.txt" dedx_simtool.gas_file ="he_90_isobutane_10.gas" dedx_simtool.IonMobility_file ="IonMobility_He+_He.txt" - dedx_simtool.save_mc = True##IF this is False then ... - dedx_simtool.debug = False + dedx_simtool.save_mc = True + dedx_simtool.debug = False dedx_simtool.sim_pulse = True #dedx_simtool.model='/junofs/users/wxfang/MyGit/tmp/fork_cepcsw_20220418/CEPCSW/Digitisers/SimCurrentONNX/src/model_90He10C4H10_18mm.onnx' dedx_simtool.model='model_90He10C4H10_18mm.onnx' diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp index b6ab38e67..0c687f827 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp @@ -87,8 +87,6 @@ Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) { // reset m_track2primary.clear(); - - auto SimPIonCol = m_SimPrimaryIonizationCol.createAndPut(); } @@ -97,8 +95,6 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) { msg() << "mcCol size (after simulation) : " << mcCol->size() << endmsg; // save all data - auto SimPrimaryIonizationCol = m_SimPrimaryIonizationCol.get(); - //msg() << "SimPrimaryIonizationCol size ="<size()< { @@ -130,8 +129,6 @@ class Edm4hepWriterAnaElemTool: public extends { "DriftChamberHitsCollection", Gaudi::DataHandle::Writer, this}; - // for ionized electron - DataHandle m_SimPrimaryIonizationCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Writer, this}; private: // in order to associate the hit contribution with the primary track, diff --git a/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp b/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp index 3bafda261..f66c582b2 100644 --- a/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp +++ b/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp @@ -25,7 +25,10 @@ double TrackHeedSimTool::dndx(double betagamma) {return 0;} double TrackHeedSimTool::dedx(const G4Step* Step) { - + if(m_beginEvt){ + m_SimPrimaryIonizationCol = m_SimPrimaryIonizationColWriter.createAndPut(); + m_beginEvt = false; + } clock_t t0 = clock(); double de = 0; float cm_to_mm = 10; @@ -43,15 +46,12 @@ double TrackHeedSimTool::dedx(const G4Step* Step) if(g4Track->GetKineticEnergy() <=0) return 0; if(pdg_code == 11 && (tmp_str_pro=="phot" || tmp_str_pro=="hIoni" || tmp_str_pro=="eIoni" || tmp_str_pro=="muIoni" || tmp_str_pro=="ionIoni" ) ) return m_eps;//skip the electron produced by Ioni, because it is already simulated by TrackHeed if(m_particle_map.find(pdg_code) == m_particle_map.end() ) return m_eps; - edm4hep::SimPrimaryIonizationClusterCollection* SimPrimaryIonizationCol = nullptr; edm4hep::MCParticleCollection* mcCol = nullptr; try{ - SimPrimaryIonizationCol = const_cast(m_SimPrimaryIonizationCol.get()); mcCol = const_cast(m_mc_handle.get()); } catch(...){ - G4cout<<"Error! Can't find collection in event, please check it have been createAndPut() in Begin of event"<GetCluster(xc, yc, zc, tc, nc, ec, extra)) { //auto chit = SimHitCol->create(); - auto chit = SimPrimaryIonizationCol->create(); + //auto chit = SimPrimaryIonizationCol->create(); + auto chit = m_SimPrimaryIonizationCol->create(); chit.setTime(tc); double cpos[3] = { cm_to_mm*( (xc - init_x)+position_x/CLHEP::cm) , cm_to_mm*((yc - init_y)+position_y/CLHEP::cm), cm_to_mm*(zc + position_z/CLHEP::cm)}; chit.setPosition(edm4hep::Vector3d(cpos)); @@ -395,7 +396,8 @@ StatusCode TrackHeedSimTool::initialize() m_current_Parent_ID = -1; m_change_track = false; m_total_range = 0; - m_isFirst = false; + m_isFirst = true; + m_beginEvt = true; m_tot_edep = 0; m_tot_length = 0; m_pa_KE =0; @@ -555,26 +557,16 @@ float* TrackHeedSimTool::NNPred(std::vector& inputs) void TrackHeedSimTool::endOfEvent() { if(m_sim_pulse){ - edm4hep::SimPrimaryIonizationClusterCollection* SimPrimaryIonizationCol = nullptr; - try{ - SimPrimaryIonizationCol = const_cast(m_SimPrimaryIonizationCol.get()); - } - catch(...){ - G4cout<<"Error! Can't find collection in event, please check it have been createAndPut() in Begin of event"< { double dndx(double betagamma) override; void getMom(float ee, float dx, float dy,float dz, float mom[3] ); void reset(){ + m_beginEvt = true; m_isFirst = true; m_previous_track_ID = 0; m_previous_KE = 0; m_tot_edep = 0; - //std::cout<<"m_tot_length="< { Gaudi::Property m_BField {this, "BField", -3}; Gaudi::Property m_eps { this, "eps" , 1e-6 };//very small value, it is returned dedx for unsimulated step (may needed for SimTrackerHit) // Output collections - DataHandle m_SimPrimaryIonizationCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Writer, this}; + DataHandle m_SimPrimaryIonizationColWriter{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Writer, this}; + edm4hep::SimPrimaryIonizationClusterCollection* m_SimPrimaryIonizationCol; // In order to associate MCParticle with contribution, we need to access MC Particle. DataHandle m_mc_handle{"MCParticle", Gaudi::DataHandle::Writer, this}; @@ -102,19 +103,20 @@ class TrackHeedSimTool: public extends { Sensor* m_sensor; std::map m_particle_map; - int m_previous_track_ID=0; - float m_previous_KE=0; + int m_previous_track_ID; + float m_previous_KE; int m_current_track_ID; int m_current_Parent_ID; int m_pdg_code; G4StepPoint* m_pre_point; G4StepPoint* m_post_point; G4double m_total_range; - bool m_isFirst=true; + bool m_isFirst; + bool m_beginEvt; bool m_change_track; edm4hep::MCParticle m_mc_paricle; - float m_tot_edep=0; - float m_tot_length=0; + float m_tot_edep; + float m_tot_length; float m_pa_KE; G4double m_pre_x ; @@ -147,12 +149,12 @@ class TrackHeedSimTool: public extends { Gaudi::Property m_sim_pulse { this, "sim_pulse" , true }; Gaudi::Property m_model_file{ this, "model", "model_test.onnx"}; Gaudi::Property m_batchsize { this, "batchsize", 100}; - Gaudi::Property m_time_scale { this, "time_scale", 99.0}; - Gaudi::Property m_time_shift { this, "time_shift", 166.4}; - Gaudi::Property m_amp_scale { this, "amp_scale" , 1e-2 }; - Gaudi::Property m_amp_shift { this, "amp_shift" , 0 }; - Gaudi::Property m_x_scale { this, "x_scale" , 5. };// in mm - Gaudi::Property m_y_scale { this, "y_scale" , 5. };// in mm + Gaudi::Property m_time_scale { this, "time_scale", 503.0}; + Gaudi::Property m_time_shift { this, "time_shift", 814.0}; + Gaudi::Property m_amp_scale { this, "amp_scale" , 1.15 }; + Gaudi::Property m_amp_shift { this, "amp_shift" , 0.86 }; + Gaudi::Property m_x_scale { this, "x_scale" , 9. };// in mm + Gaudi::Property m_y_scale { this, "y_scale" , 9. };// in mm diff --git a/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h b/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h index d4e8e9e4a..4da01d532 100644 --- a/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h +++ b/Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h @@ -28,7 +28,6 @@ class IDedxSimTool: virtual public IAlgTool { virtual double dedx(const G4Step* aStep) = 0; virtual double dedx(const edm4hep::MCParticle& mc) = 0; virtual double dndx(double betagamma) = 0; - virtual void reset() {} virtual void endOfEvent() {} }; diff --git a/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp index c1d768021..44bf01d10 100644 --- a/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp +++ b/Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp @@ -72,5 +72,4 @@ DriftChamberSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) { void DriftChamberSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE) { m_dedx_simtool->endOfEvent(); - m_dedx_simtool->reset(); }