diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/sim_Heed.py b/Detector/DetCRD/scripts/TDR_o1_v01/sim_Heed.py
new file mode 100644
index 0000000000000000000000000000000000000000..aefaa8643254dd9c1c311f34adb941543abf892b
--- /dev/null
+++ b/Detector/DetCRD/scripts/TDR_o1_v01/sim_Heed.py
@@ -0,0 +1,116 @@
+#!/usr/bin/env python
+import os
+from Gaudi.Configuration import *
+
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc")
+
+from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
+seed = [12340]
+# 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()
+
+# option for standalone tracker study
+geometry_option = "TDR_o1_v01/TDR_o1_v01.xml"
+
+if not os.getenv("DETCRDROOT"):
+ print("Can't find the geometry. Please setup envvar DETCRDROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCRDROOT"), "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.PositionXs = [0]
+gun.PositionYs = [0]
+gun.PositionZs = [0]
+gun.Particles = ["mu-"]
+gun.EnergyMins = [1]
+gun.EnergyMaxs = [100]
+gun.ThetaMins = [8]
+gun.ThetaMaxs = [172]
+gun.PhiMins = [0]
+gun.PhiMaxs = [360]
+
+genprinter = GenPrinter("GenPrinter")
+
+genalg = GenAlgo("GenAlgo")
+genalg.GenTools = ["GtGunTool"]
+
+##############################################################################
+# Detector Simulation
+##############################################################################
+from Configurables import DetSimSvc
+detsimsvc = DetSimSvc("DetSimSvc")
+
+from Configurables import DetSimAlg
+detsimalg = DetSimAlg("DetSimAlg")
+detsimalg.RandomSeeds = seed
+# detsimalg.VisMacs = ["vis.mac"]
+detsimalg.RunCmds = [
+# "/tracking/verbose 1",
+]
+detsimalg.AnaElems = [
+ # example_anatool.name()
+ # "ExampleAnaElemTool"
+ "Edm4hepWriterAnaElemTool"
+]
+detsimalg.RootDetElem = "WorldDetElemTool"
+
+from Configurables import TimeProjectionChamberSensDetTool
+tpc_sensdettool = TimeProjectionChamberSensDetTool("TimeProjectionChamberSensDetTool")
+tpc_sensdettool.TypeOption = 1
+tpc_sensdettool.DoHeedSim = True
+dedxoption = "TrackHeedSimTool"
+tpc_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import TrackHeedSimTool
+dedx_simtool = TrackHeedSimTool("TrackHeedSimTool")
+dedx_simtool.detector = "TPC"
+dedx_simtool.only_primary = False#True
+dedx_simtool.use_max_step = False#True
+dedx_simtool.max_step = 1#mm
+dedx_simtool.save_mc = True
+
+
+
+from Configurables import MarlinEvtSeeder
+evtseeder = MarlinEvtSeeder("EventSeeder")
+
+# output
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "simHeed00.root"
+out.outputCommands = ["keep *"]
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+mgr = ApplicationMgr(
+ TopAlg = [genalg, detsimalg, out],
+ EvtSel = 'NONE',
+ EvtMax = 10,
+ ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc],
+ HistogramPersistency = 'ROOT',
+ OutputLevel = ERROR
+)
diff --git a/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp b/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp
index f66c582b2a90d19ad86ede717d7e538ca51d9de0..a4fe56c70f84ea3c7b1d27dc97f2d0645ee3db2f 100644
--- a/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp
+++ b/Simulation/DetSimDedx/src/TrackHeedSimTool.cpp
@@ -142,19 +142,17 @@ double TrackHeedSimTool::dedx(const G4Step* Step)
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);
+ m_track->TransportDeltaElectron(m_init_x, m_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)};
+ double epos[3] = { cm_to_mm*( (xe - m_init_x)+position_x/CLHEP::cm) , cm_to_mm*((ye - m_init_y)+position_y/CLHEP::cm), cm_to_mm*(ze + position_z/CLHEP::cm)};
ehit.setPosition(edm4hep::Vector3d(epos));
ehit.setType(11);
}
@@ -186,7 +184,13 @@ double TrackHeedSimTool::dedx(const G4Step* Step)
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
+ m_tv3.SetX(track_dx);
+ m_tv3.SetY(track_dy);
+ m_tv3.SetZ(track_dz);
+ if(m_det=="TPC"){
+ doXRotation(false, position_z, m_tv3);
+ }
+ m_track->NewTrack(m_init_x, m_init_y, 0, track_time/CLHEP::ns, m_tv3.X(), m_tv3.Y(), m_tv3.Z());//cm
double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
int nc = 0;
int ic = 0;
@@ -197,12 +201,18 @@ double TrackHeedSimTool::dedx(const G4Step* Step)
//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)};
+ m_tv3.SetX(xc - m_init_x);
+ m_tv3.SetY(yc - m_init_y);
+ m_tv3.SetZ(zc);
+ if(m_det=="TPC"){
+ doXRotation(true, position_z, m_tv3);
+ }
+ double cpos[3] = { cm_to_mm*( m_tv3.X()+position_x/CLHEP::cm) , cm_to_mm*(m_tv3.Y()+position_y/CLHEP::cm), cm_to_mm*(m_tv3.Z() + 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_cellID && m_det=="DC") 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="<<mcCol->at(track_ID-1).getObjectID().index<<std::endl;
@@ -217,8 +227,14 @@ double TrackHeedSimTool::dedx(const G4Step* Step)
//auto ehit = SimIonizationCol->create();
//ehit.setPrimaryIonization(chit);
chit.addToElectronTime(te);
+ m_tv3.SetX(xe - m_init_x);
+ m_tv3.SetY(ye - m_init_y);
+ m_tv3.SetZ(ze);
+ if(m_det=="TPC"){
+ doXRotation(true, position_z, m_tv3);
+ }
//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)};
+ double epos[3] = { cm_to_mm*( m_tv3.X()+position_x/CLHEP::cm) , cm_to_mm*(m_tv3.Y()+position_y/CLHEP::cm), cm_to_mm*(m_tv3.Z() + position_z/CLHEP::cm)};
//ehit.setPosition(edm4hep::Vector3d(epos));
//ehit.setPosition(edm4hep::Vector3d(epos));
chit.addToElectronPosition(edm4hep::Vector3d(epos));
@@ -233,7 +249,7 @@ double TrackHeedSimTool::dedx(const G4Step* Step)
ehit.setMomentum(edm4hep::Vector3f(emom));
*/
//if(m_save_cellID) ehit.setCellID( getCellID(epos[0], epos[1], epos[2]) );
- if(m_save_cellID) chit.addToElectronCellID( getCellID(epos[0], epos[1], epos[2]) );
+ if(m_save_cellID && m_det=="DC") chit.addToElectronCellID( getCellID(epos[0], epos[1], epos[2]) );
//ehit.setQuality(2);
//ehit.setType(0);//default
}
@@ -300,15 +316,16 @@ void TrackHeedSimTool::getMom(float ee, float dx, float dy,float dz, float mom[3
StatusCode TrackHeedSimTool::initialize()
{
-
- m_geosvc = service<IGeomSvc>("GeomSvc");
- if ( !m_geosvc ) throw "TrackHeedSimTool :Failed to find GeomSvc ...";
- m_dd4hep = m_geosvc->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<dd4hep::DDSegmentation::GridDriftChamber*>(m_readout->segmentation().segmentation());
- if ( !m_segmentation ) throw "TrackHeedSimTool :Failed to get segmentation ...";
+ if(m_det=="DC"){
+ m_geosvc = service<IGeomSvc>("GeomSvc");
+ if ( !m_geosvc ) throw "TrackHeedSimTool :Failed to find GeomSvc ...";
+ m_dd4hep = m_geosvc->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<dd4hep::DDSegmentation::GridDriftChamber*>(m_readout->segmentation().segmentation());
+ if ( !m_segmentation ) throw "TrackHeedSimTool :Failed to get segmentation ...";
+ }
m_particle_map[ 11] = "e-";
m_particle_map[-11] = "e+";
@@ -323,62 +340,85 @@ StatusCode TrackHeedSimTool::initialize()
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");
- */
+
+ if(m_det=="DC"){
+ 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");
+ */
+ }
+ else if(m_det=="TPC"){
+ m_gas.SetTemperature(293.15);
+ m_gas.SetPressure(750.);
+ m_gas.SetComposition("Ar", 95., "isobutane", 2., "CF4", 3.);
+ }
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");
-
+
+ if(m_det=="DC"){
+ m_init_x = 10;//cm
+ m_init_y = -10;//cm
+ // 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");
+ }
+ else if(m_det=="TPC"){
+ //鍦℅arfield涓病鏈堿ddPlaneZ杩欎竴閫夐」锛屽洜姝ゅ綋鍓嶇殑Y鍧愭爣涓哄叏灞€涓殑Z鍧愭爣
+ //鍗充负鏉熸祦鏂瑰悜(鍚庨潰浼氳繘琛屽潗鏍囧彉鎹�)
+ cmp.AddPlaneY(0, -60000., "pad_plane");
+ cmp.AddPlaneY(300, 0., "HV");
+ cmp.SetMagneticField(0, m_BField, 0);
+ m_init_x = 0;//cm
+ m_init_y = 150;//cm
+ }
///
/// 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();
+ if(m_det=="DC"){
+ 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();
+ }
@@ -411,84 +451,105 @@ StatusCode TrackHeedSimTool::initialize()
m_pre_t = 0;
- // for NN pulse simulation//
- m_env = std::make_shared<Ort::Env>(ORT_LOGGING_LEVEL_WARNING, "ENV");
- m_seesion_options = std::make_shared<Ort::SessionOptions>();
- 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<Ort::Session>(*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) {
- #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();
- 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 << "]"
- #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;
- }
- // 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;
-
- }
+ if(m_det=="DC" && m_sim_pulse){
+ // for NN pulse simulation//
+ m_env = std::make_shared<Ort::Env>(ORT_LOGGING_LEVEL_WARNING, "ENV");
+ m_seesion_options = std::make_shared<Ort::SessionOptions>();
+ 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<Ort::Session>(*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) {
+ #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();
+ 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 << "]"
+ #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;
+ }
+ // 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;
+
+ }
+ }
return StatusCode::SUCCESS;
}
+void TrackHeedSimTool::doXRotation(bool transfer_back, float z_real, TVector3& v3){
+ if(transfer_back){
+ if(z_real>0){
+ v3.RotateX(-TMath::Pi()/2);
+ }
+ else{
+ v3.RotateX(TMath::Pi()/2);
+ }
+ }
+ else{
+ if(z_real>0){
+ v3.RotateX(TMath::Pi()/2);
+ }
+ else{
+ v3.RotateX(-TMath::Pi()/2);
+ }
+ }
+}
+
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)
@@ -556,7 +617,7 @@ float* TrackHeedSimTool::NNPred(std::vector<float>& inputs)
void TrackHeedSimTool::endOfEvent() {
- if(m_sim_pulse){
+ if(m_sim_pulse && m_det=="DC"){
if(m_debug) G4cout<<"SimPrimaryIonizationCol size="<<m_SimPrimaryIonizationCol->size()<<G4endl;
clock_t t01 = clock();
std::vector<float> inputs;
diff --git a/Simulation/DetSimDedx/src/TrackHeedSimTool.h b/Simulation/DetSimDedx/src/TrackHeedSimTool.h
index 4568c4f9016842d9aa5d8ab017a88e7cdfc90447..eb66b4ed09234dff8a4601165c3e9c01bcf5d818 100644
--- a/Simulation/DetSimDedx/src/TrackHeedSimTool.h
+++ b/Simulation/DetSimDedx/src/TrackHeedSimTool.h
@@ -51,7 +51,6 @@ class TrackHeedSimTool: public extends<AlgTool, IDedxSimTool> {
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;
@@ -61,17 +60,20 @@ class TrackHeedSimTool: public extends<AlgTool, IDedxSimTool> {
m_tot_length = 0;
}
void endOfEvent();
+ private:
+ void doXRotation(bool transfer_back, float z_real, TVector3& v3);
+ void getMom(float ee, float dx, float dy,float dz, float mom[3] );
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:
+ float xy2phi(float x, float y);
+ float* NNPred(std::vector<float>& inputs);
//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_det {this, "detector", "DC"};//choice ["DC","TPC"]
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"};
@@ -126,6 +128,9 @@ class TrackHeedSimTool: public extends<AlgTool, IDedxSimTool> {
G4double m_pre_dy ;
G4double m_pre_dz ;
G4double m_pre_t ;
+ float m_init_x;//cm
+ float m_init_y;//cm
+ TVector3 m_tv3;
//// sim pulse from NN ///
Gaudi::Property<int> m_intra_op_nthreads{ this, "intraOpNumThreads", 1};
diff --git a/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.cpp b/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.cpp
index 94c1e223aa124d6f739acfa1eb2d829400148115..eeab06e7c0b02aa15823e8cf7fc33d9d640d7d9c 100644
--- a/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.cpp
+++ b/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.cpp
@@ -20,6 +20,14 @@ StatusCode TimeProjectionChamberSensDetTool::initialize() {
error() << "Failed to find GeomSvc." << endmsg;
return StatusCode::FAILURE;
}
+
+ if(m_do_heed_sim){
+ m_dedx_simtool = ToolHandle<IDedxSimTool>(m_dedx_sim_option.value());
+ if (!m_dedx_simtool) {
+ error() << "Failed to find dedx simtoo." << endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
return AlgTool::initialize();
}
@@ -45,6 +53,8 @@ G4VSensitiveDetector* TimeProjectionChamberSensDetTool::createSD(const std::stri
tpcsd->setWriteMCTruthForLowPtHits(m_writeMCTruthForLowPtHits);
tpcsd->setLowPtCut(m_lowPtCut/dd4hep::MeV*CLHEP::MeV);
tpcsd->setLowPtMaxHitSeparation(m_lowPtMaxHitSeparation/dd4hep::mm*CLHEP::mm);
+ tpcsd->setDedxSimTool(m_dedx_simtool);
+ tpcsd->setDoHeedSim(m_do_heed_sim);
sd = tpcsd;
info() << "TPC will use TimeProjectionChamberSensitiveDetector" << endmsg;
}
diff --git a/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.h b/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.h
index 89a8593ec99f2a54a059075e9096af3794ed669a..9ee751f832c4014a9207f1f71214cf40bee304f4 100644
--- a/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.h
+++ b/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.h
@@ -8,6 +8,7 @@
#include "GaudiKernel/AlgTool.h"
#include "GaudiKernel/ToolHandle.h"
#include "DetSimInterface/ISensDetTool.h"
+#include "DetSimInterface/IDedxSimTool.h"
#include "DetInterface/IGeomSvc.h"
#include "DD4hep/DD4hepUnits.h"
@@ -35,6 +36,9 @@ private:
Gaudi::Property<double> m_lowPtMaxHitSeparation{this, "LowPtMaxHitSeparation", 5*dd4hep::mm};
Gaudi::Property<bool> m_sameStepLimit{this, "SameStepLimit", true};
Gaudi::Property<bool> m_writeMCTruthForLowPtHits{this, "WriteMCTruthForLowPtHits", false};
+ ToolHandle<IDedxSimTool> m_dedx_simtool;
+ Gaudi::Property<std::string> m_dedx_sim_option{this, "DedxSimTool",""};
+ Gaudi::Property<bool> m_do_heed_sim{this, "DoHeedSim", false};
};
#endif
diff --git a/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.cpp b/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.cpp
index fdb1ff09ce8b5fff32a2929b90d5255db8858c15..4455f3876a6545e47e812c7cc130b6490cdc3721 100644
--- a/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.cpp
+++ b/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.cpp
@@ -21,6 +21,12 @@ TimeProjectionChamberSensitiveDetector::TimeProjectionChamberSensitiveDetector(c
collectionName.insert(CollName3);
}
+bool TimeProjectionChamberSensitiveDetector::setDedxSimTool(ToolHandle<IDedxSimTool> simtool) {
+ m_dedx_simtool = simtool;
+
+ return true;
+}
+
void TimeProjectionChamberSensitiveDetector::Initialize(G4HCofThisEvent* HCE){
m_hc = new HitCollection(GetName(), collectionName[0]);
int HCID = G4SDManager::GetSDMpointer()->GetCollectionID(m_hc);
@@ -82,6 +88,10 @@ G4bool TimeProjectionChamberSensitiveDetector::ProcessHits(G4Step* step, G4Touch
// (ii) a particle that crosses the boundry between two pad-ring halves will have the hit
// placed on this surface at the last crossing point, and will be assinged the total energy
// deposited in the whole pad-ring. This is a possible source of bias for the hit
+ double dedx = 0.0;
+ if(m_doHeedSim){
+ dedx = m_dedx_simtool->dedx(step);
+ }
G4TouchableHandle touchPost = step->GetPostStepPoint()->GetTouchableHandle();
G4TouchableHandle touchPre = step->GetPreStepPoint()->GetTouchableHandle();
@@ -278,6 +288,9 @@ void TimeProjectionChamberSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE){
if (CumulativeEnergyDeposit > m_thresholdEnergyDeposit) {
DepositLowPtHit();
}
+ if(m_doHeedSim) {
+ m_dedx_simtool->endOfEvent();
+ }
}
void TimeProjectionChamberSensitiveDetector::DepositLowPtHit(){
diff --git a/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.h b/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.h
index bcb91b435f2506d5ee6cf1eda8c61a12f26ff2bf..c8c4cdb633db81372e119453345bfdfe1c0454e1 100644
--- a/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.h
+++ b/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.h
@@ -9,6 +9,8 @@
#include "DetSimSD/DDG4SensitiveDetector.h"
#include "DDG4/Defs.h"
+#include "DetSimInterface/IDedxSimTool.h"
+#include "GaudiKernel/ToolHandle.h"
class TimeProjectionChamberSensitiveDetector: public DDG4SensitiveDetector {
public:
@@ -33,6 +35,8 @@ class TimeProjectionChamberSensitiveDetector: public DDG4SensitiveDetector {
/// helper function to avoid code duplication,
/// adds energy, track length and momentum of a low pt step to the cumulative variables
void CumulateLowPtStep(dd4hep::sim::Geant4StepHandler& h);
+ bool setDedxSimTool(ToolHandle<IDedxSimTool>);
+ void setDoHeedSim(bool doHeedSim){m_doHeedSim=doHeedSim;};
protected:
@@ -63,6 +67,8 @@ class TimeProjectionChamberSensitiveDetector: public DDG4SensitiveDetector {
int CurrentTrackID; //< the TrackID of the particle causing the cumulative energy deposit
double CurrentGlobalTime; ///< the global time of the track causing the cumulative energy deposit
int CurrentCopyNumber; ///< copy number of the preStepPoint's TouchableHandle for the cumulative energy deposit
+ ToolHandle<IDedxSimTool> m_dedx_simtool;
+ bool m_doHeedSim = false;
};
#endif