Skip to content

Commit 0d40e2a

Browse files
kjplowskjplows
committed
Merge branch 'feature/pgreen_g4rw' of https://github.com/pgreen135/sbncode into rc-nuint
2 parents a96cdba + a2e7fbd commit 0d40e2a

19 files changed

+1001
-9
lines changed

‎CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,7 @@ find_package( CLHEP REQUIRED )
7777
find_package( ROOT REQUIRED )
7878
find_package( Geant4 REQUIRED )
7979
find_package( Boost COMPONENTS system filesystem REQUIRED )
80+
find_package( geant4reweight REQUIRED )
8081

8182
# macros for dictionary and simple_plugin
8283
include(ArtDictionary)

‎fcl/caf/cafmaker_common_defs.fcl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
#include "eventweight_geant4_sbn.fcl"
12
#include "eventweight_genie_sbn.fcl"
23
#include "eventweight_flux_sbn.fcl"
34

@@ -15,6 +16,7 @@ cafmaker_common_producers: {
1516
rns: { module_type: "RandomNumberSaver" }
1617
genieweight: @local::sbn_eventweight_genie
1718
fluxweight: @local::sbn_eventweight_flux
19+
geant4weight: @local::sbn_eventweight_geant4
1820
}
1921

2022
END_PROLOG

‎sbncode/SBNEventWeight/App/CMakeLists.txt‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,11 @@
11
cet_build_plugin(SBNEventWeight art::module
22
LIBRARIES
3-
sbnobj::Common_SBNEventWeight sbncode_SBNEventWeight_Base
3+
sbnobj::Common_SBNEventWeight
4+
sbncode_SBNEventWeight_Base
45
sbncode_SBNEventWeight_Calculators_CrossSection
5-
sbncode_SBNEventWeight_Calculators_BNBFlux nusimdata::SimulationBase
6+
sbncode_SBNEventWeight_Calculators_BNBFlux
7+
sbncode_SBNEventWeight_Calculators_Geant4
8+
nusimdata::SimulationBase
69
ROOT::Geom)
710

811
cet_build_plugin(SystToolsEventWeight art::module

‎sbncode/SBNEventWeight/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,6 @@
1+
include_directories ( $ENV{GEANT4REWEIGHT_INC} )
2+
link_directories ( $ENV{GEANT4REWEIGHT_LIB} )
3+
14
add_subdirectory(Base)
25
add_subdirectory(Calculators)
36
add_subdirectory(App)

‎sbncode/SBNEventWeight/Calculators/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
11
add_subdirectory(CrossSections)
22
add_subdirectory(BNBFlux)
3-
#add_subdirectory(Geant4)
3+
add_subdirectory(Geant4)
44

‎sbncode/SBNEventWeight/Calculators/Geant4/BetheBlochForG4ReweightValid.h‎

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
2+
3+
double BetheBloch(double energy, double mass){
4+
5+
//Need to make this configurable? Or delete...
6+
double K = .307075;
7+
double rho = 1.390;
8+
double Z = 18;
9+
double A = 40;
10+
double I = 188E-6;
11+
double me = .511;
12+
//Need to make sure this is total energy, not KE
13+
double gamma = energy/mass;
14+
double beta = sqrt( 1. - (1. / (gamma*gamma)) );
15+
double Tmax = 2 * me * beta*beta * gamma*gamma;
16+
17+
double first = K * (Z/A) * rho / (beta*beta);
18+
double second = .5 * log(Tmax*Tmax/(I*I)) - beta*beta;
19+
20+
double dEdX = first*second;
21+
return dEdX;
22+
}
23+
24+
std::vector< std::pair<double, int> > ThinSliceBetheBloch(G4ReweightTraj * theTraj, double res, double mass, bool isElastic){
25+
26+
std::vector< std::pair<double, int> > result;
27+
28+
//First slice position
29+
// double sliceEdge = res;
30+
// double lastPos = 0.;
31+
// double nextPos = 0.;
32+
// double px,py,pz;
33+
int interactInSlice = 0;
34+
35+
//Get total distance traveled in z
36+
// double totalDeltaZ = 0.;
37+
double disp = 0.;
38+
// double oldDisp = 0.;
39+
// int crossedSlices = 0;
40+
41+
int currentSlice = 0;
42+
int oldSlice = 0;
43+
44+
double sliceEnergy = theTraj->GetEnergy();
45+
size_t nSteps = theTraj->GetNSteps();
46+
for(size_t is = 0; is < nSteps; ++is){
47+
auto theStep = theTraj->GetStep(is);
48+
49+
disp += theStep->GetStepLength();
50+
currentSlice = floor(disp/res);
51+
52+
std::string theProc = theStep->GetStepChosenProc();
53+
54+
//Check to see if in a new slice and it's not the end
55+
if( oldSlice != currentSlice && is < nSteps - 1){
56+
57+
58+
//Save Interaction info of the prev slice
59+
//and reset
60+
result.push_back( std::make_pair(sliceEnergy, interactInSlice) );
61+
interactInSlice = 0;
62+
63+
//Update the energy
64+
sliceEnergy = sliceEnergy - res*BetheBloch(sliceEnergy, mass);
65+
if( sliceEnergy - mass < 0.){
66+
//std::cout << "Warning! Negative energy " << sliceEnergy - mass << std::endl;
67+
//std::cout << "Crossed " << oldSlice - currentSlice << std::endl;
68+
sliceEnergy = 0.0001;
69+
}
70+
//If it's more than 1 slice, add in non-interacting slices
71+
for(int ic = 1; ic < abs( oldSlice - currentSlice ); ++ic){
72+
//std::cout << ic << std::endl;
73+
74+
result.push_back( std::make_pair(sliceEnergy, 0) );
75+
76+
//Update the energy again
77+
sliceEnergy = sliceEnergy - res*BetheBloch(sliceEnergy, mass);
78+
if( sliceEnergy - mass < 0.){
79+
//std::cout << "Warning! Negative energy " << sliceEnergy - mass << std::endl;
80+
//std::cout << "Crossed " << oldSlice - currentSlice << std::endl;
81+
sliceEnergy = 0.0001;
82+
}
83+
}
84+
85+
if ((!isElastic && theProc.find(std::string("Inelastic")) != std::string::npos) || (isElastic && theProc.find(std::string("hadElastic")) != std::string::npos)) {
86+
// std::cout << "found! " << theProc << '\n';
87+
interactInSlice = 1;
88+
}
89+
}
90+
//It's crossed a slice and it's the last step. Save both info
91+
else if( oldSlice != currentSlice && is == nSteps - 1 ){
92+
result.push_back( std::make_pair(sliceEnergy, interactInSlice) );
93+
interactInSlice = 0;
94+
95+
//Update the energy
96+
sliceEnergy = sliceEnergy - res*BetheBloch(sliceEnergy, mass);
97+
if( sliceEnergy - mass < 0.){
98+
//std::cout << "Warning! Negative energy " << sliceEnergy - mass << std::endl;
99+
//std::cout << "Crossed " << oldSlice - currentSlice << std::endl;
100+
sliceEnergy = 0.0001;
101+
}
102+
//If it's more than 1 slice, add in non-interacting slices
103+
for(int ic = 1; ic < abs( oldSlice - currentSlice ); ++ic){
104+
//std::cout << ic << std::endl;
105+
106+
result.push_back( std::make_pair(sliceEnergy, 0) );
107+
108+
//Update the energy again
109+
sliceEnergy = sliceEnergy - res*BetheBloch(sliceEnergy, mass);
110+
if( sliceEnergy - mass < 0.){
111+
//std::cout << "Warning! Negative energy " << sliceEnergy - mass << std::endl;
112+
//std::cout << "Crossed " << oldSlice - currentSlice << std::endl;
113+
sliceEnergy = 0.0001;
114+
}
115+
}
116+
117+
//Save the last slice
118+
if ((!isElastic && theProc.find(std::string("Inelastic")) != std::string::npos) || (isElastic && theProc.find(std::string("hadElastic")) != std::string::npos)) {
119+
// std::cout << "found! " << theProc << '\n';
120+
interactInSlice = 1;
121+
}
122+
result.push_back( std::make_pair(sliceEnergy, interactInSlice) );
123+
}
124+
//It's the end, so just save this last info
125+
else if( oldSlice == currentSlice && is == nSteps - 1 ){
126+
if ((!isElastic && theProc.find(std::string("Inelastic")) != std::string::npos) || (isElastic && theProc.find(std::string("hadElastic")) != std::string::npos)) {
127+
// std::cout << "found! " << theProc << '\n';
128+
interactInSlice = 1;
129+
}
130+
result.push_back( std::make_pair(sliceEnergy, interactInSlice) );
131+
}
132+
//Same slice, not the end. Check for interactions
133+
else{
134+
if ((!isElastic && theProc.find(std::string("Inelastic")) != std::string::npos) || (isElastic && theProc.find(std::string("hadElastic")) != std::string::npos)) {
135+
// std::cout << "found! " << theProc << '\n';
136+
interactInSlice = 1;
137+
}
138+
}
139+
140+
//Update oldslice
141+
oldSlice = currentSlice;
142+
}
143+
144+
return result;
145+
}

‎sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
include_directories ( $ENV{GEANT4_FQ_DIR}/include )
2+
include_directories ( $ENV{GEANT4REWEIGHT_INC} )
3+
link_directories ( $ENV{GEANT4REWEIGHT_LIB} )
4+
5+
art_make_library(
6+
LIBRARY_NAME sbncode_SBNEventWeight_Calculators_Geant4
7+
LIBRARIES
8+
sbncode_SBNEventWeight_Base
9+
nusimdata::SimulationBase
10+
nurandom::RandomUtils_NuRandomService_service
11+
art::Framework_Principal
12+
art::Framework_Services_Registry
13+
art_root_io::TFileService_service
14+
larcore::Geometry_Geometry_service
15+
cetlib_except::cetlib_except
16+
ReweightBaseLib
17+
PropBaseLib
18+
ROOT::Tree
19+
)
20+
21+
install_headers()
22+
install_fhicl()
23+
install_source()
24+

0 commit comments

Comments
 (0)