From 49ac9a0be06b64c79c2270fdb32dce9966c234be Mon Sep 17 00:00:00 2001
From: "tanggy@ihep.ac.cn" <tanggy@ihep.ac.cn>
Date: Thu, 9 May 2024 01:44:43 +0000
Subject: [PATCH] The function of calculating absorbed dose/dose
eq/Si-1MeV-neutron fluence/High-energy-hadron fluence
---
Examples/options/tut_dosesim.py | 162 ++++++
Simulation/DetSimAna/CMakeLists.txt | 1 +
.../DetSimAna/src/ExampleAnaDoseElemTool.cpp | 462 ++++++++++++++++++
.../DetSimAna/src/ExampleAnaDoseElemTool.h | 120 +++++
Simulation/DetSimCore/src/RunAction.cpp | 1 +
5 files changed, 746 insertions(+)
create mode 100644 Examples/options/tut_dosesim.py
create mode 100644 Simulation/DetSimAna/src/ExampleAnaDoseElemTool.cpp
create mode 100644 Simulation/DetSimAna/src/ExampleAnaDoseElemTool.h
diff --git a/Examples/options/tut_dosesim.py b/Examples/options/tut_dosesim.py
new file mode 100644
index 00000000..90d0e218
--- /dev/null
+++ b/Examples/options/tut_dosesim.py
@@ -0,0 +1,162 @@
+#!/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 = "CepC_v4-onlyVXD.xml"
+geometry_option = "CepC_v4.xml"
+
+if not os.getenv("DETCEPCV4ROOT"):
+ print("Can't find the geometry. Please setup envvar DETCEPCV4ROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCEPCV4ROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+ print("Can't find the compact geometry file: %s"%geometry_path)
+ sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+##############################################################################
+# Physics Generator
+##############################################################################
+from Configurables import GenAlgo
+from Configurables import GtGunTool
+from Configurables import StdHepRdr
+from Configurables import SLCIORdr
+from Configurables import HepMCRdr
+from Configurables import GenPrinter
+
+gun = GtGunTool("GtGunTool")
+gun.Particles = ["mu-"]
+gun.EnergyMins = [120.] # GeV
+gun.EnergyMaxs = [120.] # GeV
+
+gun.ThetaMins = [0] # rad; 45deg
+gun.ThetaMaxs = [180.] # rad; 45deg
+
+gun.PhiMins = [0] # rad; 0deg
+gun.PhiMaxs = [0.] # rad; 360deg
+
+gun.PositionXs = [0.]
+gun.PositionYs = [0.]
+gun.PositionZs = [0.]
+# 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",
+]
+detsimalg.AnaElems = [
+ # example_anatool.name()
+ # "ExampleAnaElemTool"
+ "ExampleAnaDoseElemTool"
+]
+detsimalg.RootDetElem = "WorldDetElemTool"
+
+from Configurables import ExampleAnaDoseElemTool
+dosesimtool = ExampleAnaDoseElemTool("ExampleAnaDoseElemTool")
+dosesimtool.Gridnbins = [450,1,1100]
+dosesimtool.Coormin = [0,-0.5,-550]
+dosesimtool.Coormax = [450,0.5,550]
+dosesimtool.Regionhist = [20,1.e-12,10]
+dosesimtool.filename = "testdose_"
+dosesimtool.Dosecofffilename = "dosecoffe/"
+dosesimtool.HEHadroncut = 0.02
+
+from Configurables import CalorimeterSensDetTool
+from Configurables import DriftChamberSensDetTool
+
+from Configurables import TimeProjectionChamberSensDetTool
+tpc_sensdettool = TimeProjectionChamberSensDetTool("TimeProjectionChamberSensDetTool")
+tpc_sensdettool.TypeOption = 1
+
+
+##############################################################################
+# POD I/O
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "test-dosesim.root"
+out.outputCommands = ["keep *"]
+
+##############################################################################
+# ApplicationMgr
+##############################################################################
+#from Configurables import NTupleSvc
+#ntsvc = NTupleSvc("NTupleSvc")
+#ntsvc.Output = ["MyTuples DATAFILE='testout.root' OPT='NEW' TYP='ROOT'"]
+
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [genalg, detsimalg, out],
+ EvtSel = 'NONE',
+ EvtMax = 10,
+ ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc],
+ OutputLevel=ERROR
+)
diff --git a/Simulation/DetSimAna/CMakeLists.txt b/Simulation/DetSimAna/CMakeLists.txt
index b1bd935f..6e0adf5b 100644
--- a/Simulation/DetSimAna/CMakeLists.txt
+++ b/Simulation/DetSimAna/CMakeLists.txt
@@ -4,6 +4,7 @@ include(${Geant4_USE_FILE})
gaudi_add_module(DetSimAna
SOURCES src/Edm4hepWriterAnaElemTool.cpp
+ SOURCES src/ExampleAnaDoseElemTool.cpp
LINK DetSimInterface
DetSimSDLib
${DD4hep_COMPONENT_LIBRARIES}
diff --git a/Simulation/DetSimAna/src/ExampleAnaDoseElemTool.cpp b/Simulation/DetSimAna/src/ExampleAnaDoseElemTool.cpp
new file mode 100644
index 00000000..10dc7ed3
--- /dev/null
+++ b/Simulation/DetSimAna/src/ExampleAnaDoseElemTool.cpp
@@ -0,0 +1,462 @@
+#include "ExampleAnaDoseElemTool.h"
+
+#include "G4Step.hh"
+#include "G4Event.hh"
+#include "G4THitsCollection.hh"
+#include "G4EventManager.hh"
+#include "G4TrackingManager.hh"
+#include "G4SteppingManager.hh"
+#include "G4UnitsTable.hh"
+#include "G4SystemOfUnits.hh"
+
+#include "G4Run.hh"
+#include "G4RunManager.hh"
+
+#include "DD4hep/Detector.h"
+#include "DD4hep/Plugins.h"
+#include "DDG4/Geant4Converter.h"
+#include "DDG4/Geant4Mapping.h"
+#include "DDG4/Geant4HitCollection.h"
+#include "DDG4/Geant4Data.h"
+#include "DetSimSD/Geant4Hits.h"
+
+#include "GaudiKernel/DataObject.h"
+#include "GaudiKernel/IHistogramSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/SmartDataPtr.h"
+
+DECLARE_COMPONENT(ExampleAnaDoseElemTool)
+
+void ExampleAnaDoseElemTool::Preparecoeff(const G4ThreeVector inputnbins, const G4ThreeVector inputcoormin, const G4ThreeVector inputcoormax, const G4ThreeVector inputregionhist, const std::string path, const G4double inputhehadcut)
+{
+ // nbins = G4ThreeVector(20,20,20);
+ // // coormin = G4ThreeVector(-1000.,-1000.,-1000.);
+ // // coormax = G4ThreeVector(1000.,1000.,1000.);
+ // // regionhist = G4ThreeVector(26,1.e-14,0.02);
+ nbins = inputnbins;
+ regionhist = inputregionhist;
+ coormin = inputcoormin;
+ coormax = inputcoormax;
+ vecwidth =G4ThreeVector((coormax[0]-coormin[0])/nbins[0], (coormax[1]-coormin[1])/nbins[1], (coormax[2]-coormin[2])/nbins[2]);
+ hehadcut = inputhehadcut;
+ G4cout<< "in Run.cc: begin gen vector..."<<G4endl;
+ G4cout << "Nbins: " <<nbins<<coormin<<coormax<<regionhist<<vecwidth<<G4endl;
+ G4cout << "Dose coeff. files are in "<<path<<G4endl;
+ G4cout << "High energy hadron cut: "<<hehadcut<<G4endl;
+
+ Readcoffe(path, "n_GeV_icru95joint.txt", nbinscoeffdoseeqneu, coeffneu);
+ Readcoffe(path, "photon_GeV_icru95joint.txt", nbinscoeffdoseeqphoton, coeffproton);
+ Readcoffe(path, "heion_GeV_icru95joint.txt", nbinscoeffdoseeqhe, coeffele);
+ Readcoffe(path, "proton_GeV_icru95joint.txt", nbinscoeffdoseeqproton, coeffpos);
+ Readcoffe(path, "electron_GeV_icru95joint.txt", nbinscoeffdoseeqele, coeffhe);
+ Readcoffe(path, "positron_GeV_icru95joint.txt", nbinscoeffdoseeqpos, coeffphoton);
+ Readcoffe(path, "muonpositive_GeV_icru95joint.txt",nbinscoeffdoseeqmup, coeffmup);
+ Readcoffe(path, "muonnegative_GeV_icru95joint.txt",nbinscoeffdoseeqmum, coeffmum);
+ Readcoffe(path, "pionpositive_GeV_icru95joint.txt",nbinscoeffdoseeqpip, coeffpip);
+ Readcoffe(path, "pionnegative_GeV_icru95joint.txt",nbinscoeffdoseeqpim, coeffpim);
+ Readcoffe(path, "nflu_GeV_rd50.txt", nbinscoeff1mevneueqfluenceneu,coeff1mevneueqfluenceneu);
+ Readcoffe(path, "protonflu_GeV_rd50.txt", nbinscoeff1mevneueqfluenceproton,coeff1mevneueqfluenceproton);
+ Readcoffe(path, "pionflu_GeV_rd50.txt", nbinscoeff1mevneueqfluencepion,coeff1mevneueqfluencepion);
+ Readcoffe(path, "electronflu_GeV_rd50.txt", nbinscoeff1mevneueqfluenceele,coeff1mevneueqfluenceele);
+
+ fVector.clear();
+ fnneu.clear();
+ fdet1.clear();
+ fair.clear();
+ fdoseeqright.clear();
+ f1mevnfluen.clear();
+ fhehadfluen.clear();
+ for (G4int i = 0;i<nbins[0]*nbins[1]*nbins[2];i++){
+ (fVector).push_back(0.);
+ (fnneu).push_back(0.);
+ (f1mevnfluen).push_back(0.);
+ (fhehadfluen).push_back(0.);
+ (fdoseeqright).push_back(0.);
+ }
+ for (G4int i = 0;i<regionhist[0];i++){
+ (fdet1).push_back(0);
+ (fair).push_back(0);
+ }
+
+}
+
+void ExampleAnaDoseElemTool::Initialize(){
+ //*********************************
+ // initialize Run...
+ //********************************
+ // fRun = new Run();
+ info() <<"Grid nbins: "<<endmsg;
+ if (m_gridnbins.value().size()!=3){
+ error() <<"Grid nbins size != 3" <<endmsg;
+ }
+ for (auto tempnbins: m_gridnbins.value()) {
+ info() << tempnbins << endmsg;
+ }
+ info() <<"Coormin: "<<endmsg;
+ if (m_coormin.value().size()!=3){
+ error() <<"Coormin size != 3" <<endmsg;
+ }
+ for (auto tempnbins: m_coormin.value()) {
+ info() << tempnbins << endmsg;
+ }
+ info() <<"Coormax: "<<endmsg;
+ if (m_coormax.value().size()!=3){
+ error() <<"Coormax size != 3" <<endmsg;
+ }
+ for (auto tempnbins: m_coormax.value()) {
+ info() << tempnbins << endmsg;
+ }
+ info() <<"Regionhist: "<<endmsg;
+ if (m_regionhist.value().size()!=3){
+ error() <<"Regionhist size != 3" <<endmsg;
+ }
+ for (auto tempnbins: m_regionhist.value()) {
+ info() << tempnbins << endmsg;
+ }
+
+ G4ThreeVector innbins(m_gridnbins.value()[0],m_gridnbins.value()[1],m_gridnbins.value()[2]);
+ G4ThreeVector incoormin(m_coormin.value()[0],m_coormin.value()[1],m_coormin.value()[2]);
+ G4ThreeVector incoormax(m_coormax.value()[0],m_coormax.value()[1],m_coormax.value()[2]);
+ G4ThreeVector inregionhist(m_regionhist.value()[0],m_regionhist.value()[1],m_regionhist.value()[2]);
+ Preparecoeff(innbins,incoormin,incoormax,inregionhist,m_dosecoefffilename,m_hehadcut);
+ info() << "in initialize: finish gen vector ..."<<endmsg;
+}
+
+void
+ExampleAnaDoseElemTool::BeginOfRunAction(const G4Run*) {
+ auto msglevel = msgLevel();
+ if (msglevel == MSG::VERBOSE || msglevel == MSG::DEBUG) {
+ verboseOutput = true;
+ }
+
+ Initialize();
+ G4cout << "Begin Run of detector simultion..." << G4endl;
+}
+
+void
+ExampleAnaDoseElemTool::EndOfRunAction(const G4Run*) {
+ G4cout << "End Run of detector simultion... "<<m_filename.value() << G4endl;
+ Writefile(m_filename.value());
+}
+
+void
+ExampleAnaDoseElemTool::BeginOfEventAction(const G4Event* anEvent) {
+ msg() << "Event " << anEvent->GetEventID() << endmsg;
+}
+
+void
+ExampleAnaDoseElemTool::EndOfEventAction(const G4Event* anEvent) {
+ msg() << "Event " << anEvent->GetEventID() << endmsg;
+ // save all data
+}
+
+void
+ExampleAnaDoseElemTool::PreUserTrackingAction(const G4Track* track) {
+}
+
+void
+ExampleAnaDoseElemTool::PostUserTrackingAction(const G4Track* track) {
+}
+
+void
+ExampleAnaDoseElemTool::UserSteppingAction(const G4Step* aStep) {
+ auto aTrack = aStep->GetTrack();
+ // try to get user track info
+
+ Setdosevalue(aStep);
+}
+
+StatusCode
+ExampleAnaDoseElemTool::initialize() {
+ StatusCode sc;
+
+ return sc;
+}
+
+StatusCode
+ExampleAnaDoseElemTool::finalize() {
+ StatusCode sc;
+
+ return sc;
+}
+
+////....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void ExampleAnaDoseElemTool::Setdosevalue(const G4Step* step)
+{
+ auto steppoint=step->GetPreStepPoint();
+ auto prepoint=steppoint->GetPosition();
+ auto postpoint=step->GetPostStepPoint()->GetPosition();
+ auto momentum = steppoint->GetMomentum();
+ auto momentumdirection = steppoint->GetMomentumDirection();
+ auto totlength = step->GetStepLength();
+ auto totdepo = step->GetTotalEnergyDeposit();
+ auto density = step->GetTrack()->GetMaterial()->GetDensity()/(g/cm3);
+
+ std::vector< std::pair<G4ThreeVector,G4double> > vecpoint;
+ vecpoint.clear();
+
+ if ((prepoint[0]>=coormin[0] and prepoint[0]<=coormax[0]) and (prepoint[1]>=coormin[1] and prepoint[1]<=coormax[1]) and (prepoint[2]>=coormin[2] and prepoint[2]<=coormax[2])){
+ std::pair<G4ThreeVector,G4double> ori=std::make_pair(prepoint,0.);
+ vecpoint.push_back(ori);
+ }
+
+ G4double tempxyz,tempstepxyz;
+ G4ThreeVector temppoint;
+ std::pair<G4ThreeVector,G4double> temppair;
+
+ // x axis
+ if (momentumdirection[0]>0) {
+ for (G4int i =std::max(0,G4int(1+(prepoint[0]-coormin[0])/vecwidth[0]));i<=std::min(G4int(nbins[0]),G4int((postpoint[0]-coormin[0])/vecwidth[0]));i++){
+ tempxyz=coormin[0]+vecwidth[0]*i;
+ tempstepxyz=(tempxyz-prepoint[0])/momentumdirection[0];
+ temppoint= G4ThreeVector(tempxyz,prepoint[1]+tempstepxyz*momentumdirection[1],prepoint[2]+tempstepxyz*momentumdirection[2]);
+ if ((temppoint[0]>=coormin[0] and temppoint[0]<=coormax[0]) and (temppoint[1]>=coormin[1] and temppoint[1]<=coormax[1]) and (temppoint[2]>=coormin[2] and temppoint[2]<=coormax[2])){
+ temppair=std::make_pair(temppoint,(temppoint-prepoint).mag());
+ vecpoint.push_back(temppair);
+ // G4cout<<"in loop x: "<<i<<G4endl;
+ // G4cout<<"temppoint: "<<temppoint<<G4endl;
+ }
+ }
+ }else if(momentumdirection[0]<0) {
+ for (G4int i =std::min(G4int(nbins[0]),G4int((prepoint[0]-coormin[0])/vecwidth[0]));i>=std::max(0,G4int((postpoint[0]-coormin[0])/vecwidth[0]+1));i--){
+ tempxyz=coormin[0]+vecwidth[0]*i;
+ tempstepxyz=(tempxyz-prepoint[0])/momentumdirection[0];
+ temppoint= G4ThreeVector(tempxyz,prepoint[1]+tempstepxyz*momentumdirection[1],prepoint[2]+tempstepxyz*momentumdirection[2]);
+ if ((temppoint[0]>=coormin[0] and temppoint[0]<=coormax[0]) and (temppoint[1]>=coormin[1] and temppoint[1]<=coormax[1]) and (temppoint[2]>=coormin[2] and temppoint[2]<=coormax[2])){
+ temppair=std::make_pair(temppoint,(temppoint-prepoint).mag());
+ vecpoint.push_back(temppair);
+ //G4cout<<"in loop x: "<<i<<G4endl;
+ //G4cout<<"temppoint: "<<temppoint<<G4endl;
+ }
+ }
+ }
+
+ // y axis
+ if (momentumdirection[1]>0) {
+ for (G4int i =std::max(0,G4int(1+(prepoint[1]-coormin[1])/vecwidth[1]));i<=std::min(G4int(nbins[1]),G4int((postpoint[1]-coormin[1])/vecwidth[1]));i++){
+ tempxyz=coormin[1]+vecwidth[1]*i;
+ tempstepxyz=(tempxyz-prepoint[1])/momentumdirection[1];
+ temppoint= G4ThreeVector(prepoint[0]+tempstepxyz*momentumdirection[0],tempxyz,prepoint[2]+tempstepxyz*momentumdirection[2]);
+ if ((temppoint[0]>=coormin[0] and temppoint[0]<=coormax[0]) and (temppoint[1]>=coormin[1] and temppoint[1]<=coormax[1]) and (temppoint[2]>=coormin[2] and temppoint[2]<=coormax[2])){
+ temppair=std::make_pair(temppoint,(temppoint-prepoint).mag());
+ vecpoint.push_back(temppair);
+ // G4cout<<"in loop y: "<<i<<G4endl;
+ // G4cout<<"temppoint: "<<temppoint<<G4endl;
+ }
+ }
+ }else if(momentumdirection[1]<0) {
+ for (G4int i =std::min(G4int(nbins[1]),G4int((prepoint[1]-coormin[1])/vecwidth[1]));i>=std::max(0,G4int((postpoint[1]-coormin[1])/vecwidth[1]+1));i--){
+ tempxyz=coormin[1]+vecwidth[1]*i;
+ tempstepxyz=(tempxyz-prepoint[1])/momentumdirection[1];
+ temppoint= G4ThreeVector(prepoint[0]+tempstepxyz*momentumdirection[0],tempxyz,prepoint[2]+tempstepxyz*momentumdirection[2]);
+ if ((temppoint[0]>=coormin[0] and temppoint[0]<=coormax[0]) and (temppoint[1]>=coormin[1] and temppoint[1]<=coormax[1]) and (temppoint[2]>=coormin[2] and temppoint[2]<=coormax[2])){
+ temppair=std::make_pair(temppoint,(temppoint-prepoint).mag());
+ vecpoint.push_back(temppair);
+ // G4cout<<"in loop y: "<<i<<G4endl;
+ // G4cout<<"temppoint: "<<temppoint<<G4endl;
+ }
+ }
+ }
+
+ // z axis
+ if (momentumdirection[2]>0) {
+ for (G4int i =std::max(0,G4int(1+(prepoint[2]-coormin[2])/vecwidth[2]));i<=std::min(G4int(nbins[2]),G4int((postpoint[2]-coormin[2])/vecwidth[2]));i++){
+ // G4cout<<i<<" "<<std::max(0,G4int(1+(prepoint[2]-coormin[2])/vecwidth[2]))<<" "<<std::min(G4int(nbins[2]),G4int((postpoint[2]-coormin[2])/vecwidth[2]))<<G4endl;
+ tempxyz=coormin[2]+vecwidth[2]*i;
+ tempstepxyz=(tempxyz-prepoint[2])/momentumdirection[2];
+ temppoint= G4ThreeVector(prepoint[0]+tempstepxyz*momentumdirection[0],prepoint[1]+tempstepxyz*momentumdirection[1],tempxyz);
+ if ((temppoint[0]>=coormin[0] and temppoint[0]<=coormax[0]) and (temppoint[1]>=coormin[1] and temppoint[1]<=coormax[1]) and (temppoint[2]>=coormin[2] and temppoint[2]<=coormax[2])){
+ temppair=std::make_pair(temppoint,(temppoint-prepoint).mag());
+ vecpoint.push_back(temppair);
+ // G4cout<<"in loop z: "<<i<<G4endl;
+ // G4cout<<"temppoint: "<<temppoint<<G4endl;
+ }
+ }
+ }else if(momentumdirection[2]<0) {
+ for (G4int i =std::min(G4int(nbins[2]),G4int((prepoint[2]-coormin[2])/vecwidth[2]));i>=std::max(0,G4int((postpoint[2]-coormin[2])/vecwidth[2]+1));i--){
+ // G4cout<<i<<" "<<std::min(G4int(nbins[2]),G4int((prepoint[2]-coormin[2])/vecwidth[2]))<<" "<<std::max(0,G4int((postpoint[2]-coormin[2])/vecwidth[2]+1))<<G4endl;
+ tempxyz=coormin[2]+vecwidth[2]*i;
+ tempstepxyz=(tempxyz-prepoint[2])/momentumdirection[2];
+ temppoint= G4ThreeVector(prepoint[0]+tempstepxyz*momentumdirection[0],prepoint[1]+tempstepxyz*momentumdirection[1],tempxyz);
+ if ((temppoint[0]>=coormin[0] and temppoint[0]<=coormax[0]) and (temppoint[1]>=coormin[1] and temppoint[1]<=coormax[1]) and (temppoint[2]>=coormin[2] and temppoint[2]<=coormax[2])){
+ temppair=std::make_pair(temppoint,(temppoint-prepoint).mag());
+ vecpoint.push_back(temppair);
+ // G4cout<<"in loop z: "<<i<<G4endl;
+ // G4cout<<"temppoint: "<<temppoint<<G4endl;
+ }
+ }
+ }
+ if ((postpoint[0]>=coormin[0] and postpoint[0]<=coormax[0]) and (postpoint[1]>=coormin[1] and postpoint[1]<=coormax[1]) and (postpoint[2]>=coormin[2] and postpoint[2]<=coormax[2])){
+ std::pair<G4ThreeVector,G4double> terminal=std::make_pair(postpoint,step->GetStepLength());
+ vecpoint.push_back(terminal);
+
+ }
+
+ if (vecpoint.size()==0){return;}
+ if (vecpoint.size()==1){
+ // G4cout<<"only one point in vecpoint......"<<G4endl;
+ // G4cout<<prepoint[0]<<" "<<prepoint[1]<<" "<<prepoint[2]<<G4endl;
+ // G4cout<<postpoint[0]<<" "<<postpoint[1]<<" "<<postpoint[2]<<G4endl;
+ // G4cout<<vecpoint[0].first[0]<<" "<<vecpoint[0].first[1]<<" "<<vecpoint[0].first[2]<<G4endl;
+ vecpoint.push_back(vecpoint[0]);
+ }
+ sort(vecpoint.begin(),vecpoint.end(),cmppair);
+
+ // Run* run = static_cast<Run*>(
+ // G4RunManager::GetRunManager()->GetNonConstCurrentRun());
+ G4double ratio;
+
+ for(std::vector< std::pair<G4ThreeVector,G4double> >::iterator iter=vecpoint.begin(); iter!=vecpoint.end()-1;iter++){
+ if (totlength<=0){
+ ratio=1.;
+ }else{
+ ratio=abs(((*iter).second-(*(iter+1)).second))/totlength;
+ }
+ // G4cout<<(*iter).first[0]<<" "<<(*iter).first[1]<<" "<<(*iter).first[2]<<G4endl;
+ // G4cout<<(*(iter+1)).first[0]<<" "<<(*(iter+1)).first[1]<<" "<<(*(iter+1)).first[2]<<G4endl;
+ temppoint = ((*iter).first+(*(iter+1)).first)/2.;
+ G4int ijkbin[3];
+ for (G4int i=0;i<3;i++){
+ ijkbin[i] = G4int((temppoint[i] - coormin[i])/vecwidth[i]);
+ }
+ // G4cout<<"depo: "<<std::setprecision(8)<<ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]<<" "<<ratio*totdepo<<G4endl;
+ // if (ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]>nbins[0]*nbins[1]*nbins[2] or ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]<0){
+ if (ijkbin[0]<0 or ijkbin[0]>nbins[0] or ijkbin[1]<0 or ijkbin[1]>nbins[1] or ijkbin[2]<0 or ijkbin[2]>nbins[2]){
+ G4cout<<"out of grids....."<<G4endl;
+ G4cout<<"coor: "<<temppoint[0]<<" "<<temppoint[1]<<" "<<temppoint[2]<<G4endl;
+ G4cout<<" "<<ijkbin[0]<<" "<<ijkbin[1]<<" "<<ijkbin[2]<<G4endl;
+ G4cout<<"vecwidth: "<<vecwidth[0]<<" "<<vecwidth[1]<<" "<<vecwidth[2]<<G4endl;
+ G4cout<<"nbins: "<<nbins[0]<<" "<<nbins[1]<<" "<<nbins[2]<<G4endl;
+ G4cout<<"coormin: "<<coormin[0]<<" "<<coormin[1]<<" "<<coormin[2]<<G4endl;
+ G4cout<<"coormax: "<<coormax[0]<<" "<<coormax[1]<<" "<<coormax[2]<<G4endl;
+ G4cout<<"test: "<<(*iter).first<<" "<<(*iter).second<<G4endl;
+ G4cout<<"test: "<<(*(iter+1)).first<<" "<<(*iter).second<<G4endl;
+ }else{
+ G4double temptracklength = abs((*iter).second-(*(iter+1)).second)/cm;
+ G4double kinenergy = step->GetPostStepPoint()->GetKineticEnergy()/GeV;
+ fVector[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=ratio*(totdepo/GeV)/density;
+ fnneu[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength;
+ if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == 2112 ){
+ if(kinenergy>hehadcut) {fhehadfluen[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength;}
+ Calcudose(fdoseeqright, nbinscoeffdoseeqneu, coeffneu, ijkbin, temptracklength, kinenergy);
+ Calcudose(f1mevnfluen, nbinscoeff1mevneueqfluenceneu, coeff1mevneueqfluenceneu, ijkbin, temptracklength, kinenergy);
+ }else if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == 2212 ){
+ if(kinenergy>hehadcut) {fhehadfluen[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength;}
+ Calcudose(fdoseeqright, nbinscoeffdoseeqproton, coeffproton, ijkbin, temptracklength, kinenergy);
+ Calcudose(f1mevnfluen, nbinscoeff1mevneueqfluenceproton, coeff1mevneueqfluenceproton, ijkbin, temptracklength, kinenergy);
+ }else if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == 1000020030 or (step->GetTrack()->GetDefinition()->GetPDGEncoding()) == 1000020040 ){
+ // He3 He4
+ if(kinenergy>hehadcut) {fhehadfluen[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength;}
+ Calcudose(fdoseeqright, nbinscoeffdoseeqhe, coeffhe, ijkbin, temptracklength, kinenergy);
+ }else if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == 22 ){
+ Calcudose(fdoseeqright, nbinscoeffdoseeqphoton, coeffphoton, ijkbin, temptracklength, kinenergy);
+ }else if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == 11 ){
+ Calcudose(fdoseeqright, nbinscoeffdoseeqele, coeffele, ijkbin, temptracklength, kinenergy);
+ Calcudose(f1mevnfluen, nbinscoeff1mevneueqfluenceele, coeff1mevneueqfluenceele, ijkbin, temptracklength, kinenergy);
+ }else if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == -11 ){
+ Calcudose(fdoseeqright, nbinscoeffdoseeqpos, coeffpos, ijkbin, temptracklength, kinenergy);
+ }else if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == 13 ){
+ Calcudose(fdoseeqright, nbinscoeffdoseeqmum, coeffmum, ijkbin, temptracklength, kinenergy);
+ }else if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == -13 ){
+ Calcudose(fdoseeqright, nbinscoeffdoseeqmup, coeffmup, ijkbin, temptracklength, kinenergy);
+ }else if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == 211 ){
+ if(kinenergy>hehadcut) {fhehadfluen[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength;}
+ Calcudose(fdoseeqright, nbinscoeffdoseeqpip, coeffpip, ijkbin, temptracklength, kinenergy);
+ Calcudose(f1mevnfluen, nbinscoeff1mevneueqfluencepion, coeff1mevneueqfluencepion, ijkbin, temptracklength, kinenergy);
+ }else if ((step->GetTrack()->GetDefinition()->GetPDGEncoding()) == -211 ){
+ if(kinenergy>hehadcut) {fhehadfluen[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength;}
+ Calcudose(fdoseeqright, nbinscoeffdoseeqpim, coeffpim, ijkbin, temptracklength, kinenergy);
+ Calcudose(f1mevnfluen, nbinscoeff1mevneueqfluencepion, coeff1mevneueqfluencepion, ijkbin, temptracklength, kinenergy);
+ }else if ((abs(step->GetTrack()->GetDefinition()->GetPDGEncoding()) >= 310) and (abs(step->GetTrack()->GetDefinition()->GetPDGEncoding()) <= 321)){
+ if(kinenergy>hehadcut) {fhehadfluen[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength;}
+ }else if (abs(step->GetTrack()->GetDefinition()->GetPDGEncoding()) == 130){
+ if(kinenergy>hehadcut) {fhehadfluen[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength;}
+ }
+ }
+ }
+}
+////....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void ExampleAnaDoseElemTool::Writefile(std::string prefixion)
+{
+ G4cout << "write file paras. " <<nbins<<coormin<<coormax<<regionhist<< G4endl;
+ std::ofstream file(prefixion+"edep.dat");
+ for (G4int i = 0;i<nbins[0]*nbins[1]*nbins[2];i++){
+ file<<fVector[i]<<std::endl;
+ }
+ file.close();
+ std::ofstream file2(prefixion+"nneu.dat");
+ for (G4int i = 0;i<nbins[0]*nbins[1]*nbins[2];i++){
+ file2<<fnneu[i]<<std::endl;
+ }
+ file2.close();
+ // std::ofstream file3("det1.dat");
+ // for (G4int i = 0;i<regionhist[0];i++){
+ // file3<<fdet1[i]<<std::endl;
+ // // std::cout<<fdet1[i]<<std::endl;
+ // }
+ // file3.close();
+ // std::ofstream file4("air.dat");
+ // for (G4int i = 0;i<regionhist[0];i++){
+ // file4<<fair[i]<<std::endl;
+ // // std::cout<<fair[i]<<std::endl;
+ // }
+ // file4.close();
+ std::ofstream file5(prefixion+"doseeq.dat");
+ for (G4int i = 0;i<nbins[0]*nbins[1]*nbins[2];i++){
+ file5<<fdoseeqright[i]<<std::endl;
+ }
+ file5.close();
+ std::ofstream file6(prefixion+"1mevneqfluence.dat");
+ for (G4int i = 0;i<nbins[0]*nbins[1]*nbins[2];i++){
+ file6<<f1mevnfluen[i]<<std::endl;
+ }
+ file6.close();
+ std::ofstream file7(prefixion+"hehadfluence.dat");
+ for (G4int i = 0;i<nbins[0]*nbins[1]*nbins[2];i++){
+ file7<<fhehadfluen[i]<<std::endl;
+ }
+ file7.close();
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void ExampleAnaDoseElemTool::Readcoffe(std::string filepath, std::string filename, const G4int constnbins, G4double coeffexample[][2]){
+ // std::string filename="n_GeV_icru95joint.txt";
+ m_inputFile_doseeqneu.open(filepath+filename);
+ if (!m_inputFile_doseeqneu){
+ std::cout << "SteppingAction: PROBLEMS OPENING FILE " <<filename<< std::endl;
+ exit(0);
+ }
+ for (int i = 0; i < constnbins; i++){
+ m_inputFile_doseeqneu >> coeffexample[i][0]>> coeffexample[i][1];
+ if( coeffexample[i][0]<=0. or coeffexample[i][1]<=0.){
+ std::cout<<"There is Zero value: "<<coeffexample[i][0]<<" "<<coeffexample[i][1]<<std::endl;
+ }
+ }
+ std::cout<<"coff: "<<filename<<" "<<constnbins<<std::endl;
+ // for (int i = 0; i < constnbins; i++){
+ // std::cout<<coeffexample[i][0]<<" "<<coeffexample[i][1]<<std::endl;
+ // }
+ m_inputFile_doseeqneu.close();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void ExampleAnaDoseElemTool::Calcudose(std::vector<G4double>& fexample, const G4int constnbins, const G4double coeffexample[][2], const G4int ijkbin[], const G4double inputtracklength, const G4double inputkinenergy){
+ if(inputkinenergy >= coeffexample[constnbins-1][0]){
+ fexample[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=inputtracklength* coeffexample[constnbins-1][1];
+ // fdoseeqright[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength* coeffexample[constnbins-1][1];
+ // fdoseeqleft[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]] +=temptracklength* coeffexample[constnbins-1][1];
+ }else if(inputkinenergy <= coeffexample[0][0]){
+ fexample[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=inputtracklength* coeffexample[0][1];
+ // fdoseeqright[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength* coeffexample[0][1];
+ // fdoseeqleft[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]] +=temptracklength* coeffexample[0][1];
+ }else{
+ for(G4int jj=1;jj<nbinscoeffdoseeqneu;jj++){
+ if( inputkinenergy<coeffexample[jj][0]){
+ fexample[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=inputtracklength* coeffexample[jj][1];
+ // fdoseeqright[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]]+=temptracklength* coeffexample[jj][1];
+ // fdoseeqleft[ijkbin[0]*nbins[1]*nbins[2]+ijkbin[1]*nbins[2]+ijkbin[2]] +=temptracklength* coeffexample[jj-1][1];
+ break;
+ }
+ }
+ }
+}
diff --git a/Simulation/DetSimAna/src/ExampleAnaDoseElemTool.h b/Simulation/DetSimAna/src/ExampleAnaDoseElemTool.h
new file mode 100644
index 00000000..4c01af35
--- /dev/null
+++ b/Simulation/DetSimAna/src/ExampleAnaDoseElemTool.h
@@ -0,0 +1,120 @@
+#ifndef Edm4hepWriterAnaElemTool_h
+#define Edm4hepWriterAnaElemTool_h
+
+#include <map>
+
+#include "GaudiKernel/AlgTool.h"
+#include "k4FWCore/DataHandle.h"
+#include "DetSimInterface/IAnaElemTool.h"
+#include "DetSimInterface/CommonUserEventInfo.hh"
+#include "DetSimInterface/CommonUserTrackInfo.hh"
+
+#include "DetInterface/IGeomSvc.h"
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/CaloHitContributionCollection.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/NTuple.h"
+#include "G4ThreeVector.hh"
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+const G4int nbinscoeffdoseeqneu = 68;
+const G4int nbinscoeffdoseeqproton = 33;
+const G4int nbinscoeffdoseeqele = 49;
+const G4int nbinscoeffdoseeqpos = 49;
+const G4int nbinscoeffdoseeqhe = 24;
+const G4int nbinscoeffdoseeqphoton = 64;
+const G4int nbinscoeffdoseeqmup = 33;
+const G4int nbinscoeffdoseeqmum = 33;
+const G4int nbinscoeffdoseeqpip = 43;
+const G4int nbinscoeffdoseeqpim = 43;
+const G4int nbinscoeff1mevneueqfluenceneu = 1381;
+const G4int nbinscoeff1mevneueqfluenceproton = 927;
+const G4int nbinscoeff1mevneueqfluencepion = 900;
+const G4int nbinscoeff1mevneueqfluenceele = 15;
+
+
+
+class ExampleAnaDoseElemTool: public extends<AlgTool, IAnaElemTool> {
+
+public:
+
+ using extends::extends;
+
+ /// IAnaElemTool interface
+ // Run
+ virtual void BeginOfRunAction(const G4Run*) override;
+ virtual void EndOfRunAction(const G4Run*) override;
+
+ // Event
+ virtual void BeginOfEventAction(const G4Event*) override;
+ virtual void EndOfEventAction(const G4Event*) override;
+
+ // Tracking
+ virtual void PreUserTrackingAction(const G4Track*) override;
+ virtual void PostUserTrackingAction(const G4Track*) override;
+
+ // Stepping
+ virtual void UserSteppingAction(const G4Step*) override;
+
+
+ /// Overriding initialize and finalize
+ StatusCode initialize() override;
+ StatusCode finalize() override;
+
+
+private:
+ bool verboseOutput = false;
+ Gaudi::Property<std::vector<int>> m_gridnbins{this, "Gridnbins"};
+ Gaudi::Property<std::vector<double>> m_coormin{this, "Coormin"};
+ Gaudi::Property<std::vector<double>> m_coormax{this, "Coormax"};
+ Gaudi::Property<std::vector<double>> m_regionhist{this, "Regionhist"};
+ Gaudi::Property<double> m_hehadcut{this, "HEHadroncut"};
+ Gaudi::Property<std::string> m_filename{this, "filename"};
+ Gaudi::Property<std::string> m_dosecoefffilename{this, "Dosecofffilename","Simulation/DetSimAna/src/dosecoffe/"};
+
+ std::vector<G4double> fVector; //absorbed dose
+ std::vector<G4double> fnneu; //track length
+ std::vector<G4double> fdoseeqright; //dose eq
+ std::vector<G4double> f1mevnfluen; //1mev neutron equivalent fluence
+ std::vector<G4double> fhehadfluen; //high energy hadron fluence
+ std::vector<G4int> fdet1;
+ std::vector<G4int> fair;
+
+ G4ThreeVector nbins,regionhist,coormin,coormax,vecwidth;
+ G4double hehadcut;
+
+ void Writefile(std::string);
+ void Setdosevalue(const G4Step* );
+ void Calcudose(std::vector<G4double>& fexample, const G4int constnbins, const G4double coeffexample[][2], const G4int ijkbin[], const G4double inputtracklength, const G4double inputkinenergy);
+ void Readcoffe(std::string filepath, std::string filename, const G4int constnbins, G4double coeffexample[][2]);
+ void Initialize();
+ void Preparecoeff(const G4ThreeVector, const G4ThreeVector , const G4ThreeVector, const G4ThreeVector, const std::string, const G4double);
+
+ static int cmppair(const std::pair< G4ThreeVector,G4double > p1, const std::pair< G4ThreeVector,G4double > p2) {return p1.second<p2.second;};
+
+ G4double coeffneu[nbinscoeffdoseeqneu][2];
+ G4double coeffproton[nbinscoeffdoseeqproton][2];
+ G4double coeffele[nbinscoeffdoseeqele][2];
+ G4double coeffpos[nbinscoeffdoseeqpos][2];
+ G4double coeffhe[nbinscoeffdoseeqhe][2];
+ G4double coeffphoton[nbinscoeffdoseeqphoton][2];
+ G4double coeffmup[nbinscoeffdoseeqmup][2];
+ G4double coeffmum[nbinscoeffdoseeqmum][2];
+ G4double coeffpip[nbinscoeffdoseeqpip][2];
+ G4double coeffpim[nbinscoeffdoseeqpim][2];
+ G4double coeff1mevneueqfluenceneu[nbinscoeff1mevneueqfluenceneu][2];
+ G4double coeff1mevneueqfluenceproton[nbinscoeff1mevneueqfluenceproton][2];
+ G4double coeff1mevneueqfluencepion[nbinscoeff1mevneueqfluencepion][2];
+ G4double coeff1mevneueqfluenceele[nbinscoeff1mevneueqfluenceele][2];
+
+ std::ifstream m_inputFile_doseeqneu;
+};
+
+#endif
diff --git a/Simulation/DetSimCore/src/RunAction.cpp b/Simulation/DetSimCore/src/RunAction.cpp
index 4bf8c5bd..22eaa41c 100644
--- a/Simulation/DetSimCore/src/RunAction.cpp
+++ b/Simulation/DetSimCore/src/RunAction.cpp
@@ -12,6 +12,7 @@ RunAction::~RunAction() {
}
+
void
RunAction::BeginOfRunAction(const G4Run* aRun)
{
--
GitLab