Skip to content
Snippets Groups Projects
Commit cb54836b authored by lintao@ihep.ac.cn's avatar lintao@ihep.ac.cn
Browse files

Merge branch 'master' into 'master' 

The function of calculating absorbed dose/dose eq/Si-1MeV-neutron fluence/High-energy-hadron fluence

See merge request !8
parents b7f4f2ef 49ac9a0b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Examples/options/tut_dosesim.py 0 → 100644
+ 162
− 0
View file @ cb54836b
#!/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
)
Simulation/DetSimAna/CMakeLists.txt
+ 1
− 0
View file @ cb54836b
......@@ -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}
......
Simulation/DetSimAna/src/ExampleAnaDoseElemTool.cpp 0 → 100644
+ 462
− 0
View file @ cb54836b
#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;
}
}
}
}
Simulation/DetSimAna/src/ExampleAnaDoseElemTool.h 0 → 100644
+ 120
− 0
View file @ cb54836b
#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
Simulation/DetSimCore/src/RunAction.cpp
+ 1
− 0
View file @ cb54836b
......@@ -12,6 +12,7 @@ RunAction::~RunAction() {
}
void
RunAction::BeginOfRunAction(const G4Run* aRun)
{
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment