Merge pull request #248 from wenxingfang/master

add TrackHeedSimTool for DC cell simulation

Detector/DetDriftChamber/compact/det.xml

@@ -19,6 +19,8 @@
 
 
   <define>
+    <constant name="tracker_region_rmax" value="1723*mm" />
+    <constant name="tracker_region_zmax" value="3050*mm" />
     <constant name="world_size" value="2226*mm"/>
     <constant name="world_x" value="world_size"/>
     <constant name="world_y" value="world_size"/>

Examples/options/tut_detsim_SDT_Heed.py

+#!/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
+    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'
+    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
+)

Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp

@@ -42,6 +42,7 @@ Edm4hepWriterAnaElemTool::BeginOfRunAction(const G4Run*) {
     } else {
         error() << "Failed to find GeomSvc." << endmsg;
     }
+
 }
 
 void
@@ -86,7 +87,7 @@ Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) {
 
     // reset
     m_track2primary.clear();
-
+ 
 }
 
 void
@@ -94,7 +95,6 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) {
 
     msg() << "mcCol size (after simulation) : " << mcCol->size() << endmsg;
     // save all data
-
     // create collections.
     auto trackercols = m_trackerCol.createAndPut();
     auto calorimetercols = m_calorimeterCol.createAndPut();

Simulation/DetSimDedx/CMakeLists.txt

 
 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}")
 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
+                      #${OnnxRuntime_LIBRARY}
+                      #/cvmfs/sft.cern.ch/lcg/views/LCG_103/x86_64-centos7-gcc11-opt/lib/libonnxruntime.so
+                      ${CLHEP_LIBRARIES}
+ 
 )
 

Simulation/DetSimDedx/src/TrackHeedSimTool.h

+#ifndef TrackHeedSimTool_h
+#define TrackHeedSimTool_h
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/MsgStream.h"
+#include "DetSimInterface/IDedxSimTool.h"
+#include <GaudiKernel/AlgTool.h>
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimPrimaryIonizationClusterCollection.h"
+#include "TVector3.h"
+#include <G4StepPoint.hh>
+
+#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 <map>
+#include <string>
+
+#include "core/session/onnxruntime_cxx_api.h"
+#include "core/session/onnxruntime_c_api.h"
+using namespace Garfield;
+
+class TrackHeedSimTool: public extends<AlgTool, IDedxSimTool> {
+    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_beginEvt = true;
+            m_isFirst = true;
+            m_previous_track_ID = 0;
+            m_previous_KE = 0;
+            m_tot_edep = 0;
+            m_tot_length = 0;
+        }
+        void endOfEvent();
+        long long getCellID(float x, float y, float z);
+        void wire_xy(float x1, float y1, float z1, float x2, float y2, float z2, float z, float &x, float &y);
+        float* NNPred(std::vector<float>& inputs);
+        float xy2phi(float x, float y);
+        void getLocal(float x1, float y1, float x2, float y2, float& dx, float& dy);
+    private:
+        //ServiceHandle<IDataProviderSvc> m_eds;
+        SmartIF<IGeomSvc> m_geosvc;
+        dd4hep::Detector* m_dd4hep; 
+        dd4hep::Readout* m_readout;
+        dd4hep::DDSegmentation::GridDriftChamber* m_segmentation;
+        Gaudi::Property<std::string> m_readout_name{ this, "readout", "DriftChamberHitsCollection"};//readout for getting segmentation
+        Gaudi::Property<std::string> m_gas_file{ this, "gas_file", "He_50_isobutane_50.gas"};//gas
+        Gaudi::Property<std::string> m_IonMobility{ this, "IonMobility_file", "IonMobility_He+_He.txt"};
+        Gaudi::Property<float> m_isob  {this, "isob", 50, ""};
+        Gaudi::Property<float> m_he    {this, "he", 50, ""};
+        Gaudi::Property<bool> m_debug{this, "debug", false};
+        Gaudi::Property<bool> m_use_max_step{this, "use_max_step", false};
+        Gaudi::Property<bool> m_update_KE{this, "update_KE", true};
+        Gaudi::Property<float> m_max_step   {this, "max_step", 1};//mm
+        Gaudi::Property<bool> m_only_primary{this, "only_primary", false};
+        Gaudi::Property<bool> m_save_mc{this, "save_mc", false};
+        Gaudi::Property<bool> m_save_cellID{this, "save_cellID", true};
+        Gaudi::Property<float> m_delta_threshold{this, "delta_threshold", 50};//keV
+        Gaudi::Property<float> m_change_threshold {this, "change_threshold", 0.05};
+        Gaudi::Property<float> m_BField   {this, "BField", -3};
+        Gaudi::Property<float> m_eps     { this, "eps"   , 1e-6  };//very small value, it is returned dedx for unsimulated step (may needed for SimTrackerHit)
+        // Output collections
+        DataHandle<edm4hep::SimPrimaryIonizationClusterCollection>    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<edm4hep::MCParticleCollection> 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<int, std::string> 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_beginEvt;
+        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<int> m_intra_op_nthreads{ this, "intraOpNumThreads", 1};
+        Gaudi::Property<int> m_inter_op_nthreads{ this, "interOpNumThreads", 1};
+        std::shared_ptr<Ort::Env> m_env;
+        std::shared_ptr<Ort::SessionOptions> m_seesion_options;
+        std::shared_ptr<Ort::Session> m_session;
+        Ort::AllocatorWithDefaultOptions m_allocator;
+        std::vector<const char*> m_input_node_names;
+        std::vector<std::vector<int64_t>> m_input_node_dims;
+        std::vector<const char*> m_output_node_names;
+        std::vector<int64_t> m_output_node_dims;
+        #if (ORT_API_VERSION >=13)
+        std::vector<Ort::AllocatedStringPtr> m_inputNodeNameAllocatedStrings;
+        std::vector<Ort::AllocatedStringPtr> m_outputNodeNameAllocatedStrings;
+        #else
+        std::vector<const char*> m_inputNodeNameAllocatedStrings;
+        std::vector<const char*> m_outputNodeNameAllocatedStrings;
+        #endif
+  
+        Gaudi::Property<bool> m_sim_pulse    { this, "sim_pulse"   , true  };
+        Gaudi::Property<std::string> m_model_file{ this, "model", "model_test.onnx"};
+        Gaudi::Property<int> m_batchsize     { this, "batchsize", 100};
+        Gaudi::Property<float> m_time_scale  { this, "time_scale", 503.0};
+        Gaudi::Property<float> m_time_shift  { this, "time_shift", 814.0};
+        Gaudi::Property<float> m_amp_scale   { this, "amp_scale" , 1.15 };
+        Gaudi::Property<float> m_amp_shift   { this, "amp_shift" , 0.86 };
+        Gaudi::Property<float> m_x_scale     { this, "x_scale"   , 9.  };// in mm
+        Gaudi::Property<float> m_y_scale     { this, "y_scale"   , 9.  };// in mm
+  
+  
+
+};
+
+#endif

Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h

@@ -28,6 +28,7 @@ public:
     virtual double dedx(const G4Step* aStep) = 0;
     virtual double dedx(const edm4hep::MCParticle& mc) = 0;
     virtual double dndx(double betagamma) = 0;
+    virtual void endOfEvent() {}
 
 };
 

Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp

@@ -22,7 +22,6 @@ StatusCode DriftChamberSensDetTool::initialize() {
         error() << "Failed to find dedx simtoo." << endmsg;
         return StatusCode::FAILURE;
     }
-
     return sc;
 }
 

Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp

@@ -71,5 +71,5 @@ DriftChamberSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) {
 
 void
 DriftChamberSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE) {
-
+    m_dedx_simtool->endOfEvent();
 }

cmake/FindOnnxRuntime.cmake

+# 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)

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