fixed beam and macros

Armory/AnasenG4/G4History.macro

@@ -36,7 +36,7 @@
 /vis/scene/notifyHandlers
 /vis/scene/add/axes 0 0 0 20 mm
 /vis/scene/notifyHandlers
-/vis/viewer/set/background 1 1 1
+/vis/viewer/set/background 0 0 0
 /vis/viewer/set/style surface
 /vis/viewer/set/hiddenMarker true
 /vis/viewer/set/viewpointThetaPhi 120 150
@@ -46,7 +46,18 @@
 /vis/viewer/flush
 /vis/viewer/refresh viewer-0
 /vis/viewer/update viewer-0
-/vis/viewer/set/style w
+/control/execute run.mac
-/vis/viewer/refresh
+/control/verbose 1
-/vis/viewer/set/hiddenEdge 1
+/run/verbose 1
-/vis/viewer/refresh
+/event/verbose 0
+/tracking/verbose 0
+/run/initialize
+/vis/filtering/trajectories/create/particleFilter
+/vis/filtering/trajectories/particleFilter-0/add e-
+/vis/scene/notifyHandlers scene-0
+/vis/filtering/trajectories/particleFilter-0/add e+
+/vis/scene/notifyHandlers scene-0
+/vis/filtering/trajectories/particleFilter-0/invert true
+/vis/scene/notifyHandlers scene-0
+/run/beamOn 100
+/vis/scene/notifyHandlers scene-0

Armory/AnasenG4/build/CMakeFiles/AnasenG4.dir/DependInfo.cmake

@@ -9,6 +9,7 @@ set(CMAKE_DEPENDS_LANGUAGES
 # The set of dependency files which are needed:
 set(CMAKE_DEPENDS_DEPENDENCY_FILES
   "/home/jamesszalkie/anasen/Armory/AnasenG4/src/ActionInitialization.cc" "CMakeFiles/AnasenG4.dir/src/ActionInitialization.cc.o" "gcc" "CMakeFiles/AnasenG4.dir/src/ActionInitialization.cc.o.d"
+  "/home/jamesszalkie/anasen/Armory/AnasenG4/src/BeamReactionSteppingAction.cc" "CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o" "gcc" "CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o.d"
   "/home/jamesszalkie/anasen/Armory/AnasenG4/src/DetectorConstruction.cc" "CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o" "gcc" "CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o.d"
   "/home/jamesszalkie/anasen/Armory/AnasenG4/src/PrimaryGeneratorAction.cc" "CMakeFiles/AnasenG4.dir/src/PrimaryGeneratorAction.cc.o" "gcc" "CMakeFiles/AnasenG4.dir/src/PrimaryGeneratorAction.cc.o.d"
   "/home/jamesszalkie/anasen/Armory/AnasenG4/src/ReactionGenerator.cc" "CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o" "gcc" "CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o.d"

Armory/AnasenG4/build/CMakeFiles/AnasenG4.dir/build.make

@@ -111,10 +111,24 @@ CMakeFiles/AnasenG4.dir/src/PrimaryGeneratorAction.cc.s: cmake_force
 	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green "Compiling CXX source to assembly CMakeFiles/AnasenG4.dir/src/PrimaryGeneratorAction.cc.s"
 	/usr/bin/c++ $(CXX_DEFINES) $(CXX_INCLUDES) $(CXX_FLAGS) -S /home/jamesszalkie/anasen/Armory/AnasenG4/src/PrimaryGeneratorAction.cc -o CMakeFiles/AnasenG4.dir/src/PrimaryGeneratorAction.cc.s
 
+CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o: CMakeFiles/AnasenG4.dir/flags.make
+CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o: ../src/BeamReactionSteppingAction.cc
+CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o: CMakeFiles/AnasenG4.dir/compiler_depend.ts
+	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green --progress-dir=/home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles --progress-num=$(CMAKE_PROGRESS_4) "Building CXX object CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o"
+	/usr/bin/c++ $(CXX_DEFINES) $(CXX_INCLUDES) $(CXX_FLAGS) -MD -MT CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o -MF CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o.d -o CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o -c /home/jamesszalkie/anasen/Armory/AnasenG4/src/BeamReactionSteppingAction.cc
+
+CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.i: cmake_force
+	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green "Preprocessing CXX source to CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.i"
+	/usr/bin/c++ $(CXX_DEFINES) $(CXX_INCLUDES) $(CXX_FLAGS) -E /home/jamesszalkie/anasen/Armory/AnasenG4/src/BeamReactionSteppingAction.cc > CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.i
+
+CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.s: cmake_force
+	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green "Compiling CXX source to assembly CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.s"
+	/usr/bin/c++ $(CXX_DEFINES) $(CXX_INCLUDES) $(CXX_FLAGS) -S /home/jamesszalkie/anasen/Armory/AnasenG4/src/BeamReactionSteppingAction.cc -o CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.s
+
 CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o: CMakeFiles/AnasenG4.dir/flags.make
 CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o: ../src/ReactionGenerator.cc
 CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o: CMakeFiles/AnasenG4.dir/compiler_depend.ts
-	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green --progress-dir=/home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles --progress-num=$(CMAKE_PROGRESS_4) "Building CXX object CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o"
+	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green --progress-dir=/home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles --progress-num=$(CMAKE_PROGRESS_5) "Building CXX object CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o"
 	/usr/bin/c++ $(CXX_DEFINES) $(CXX_INCLUDES) $(CXX_FLAGS) -MD -MT CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o -MF CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o.d -o CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o -c /home/jamesszalkie/anasen/Armory/AnasenG4/src/ReactionGenerator.cc
 
 CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.i: cmake_force
@@ -128,7 +142,7 @@ CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.s: cmake_force
 CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o: CMakeFiles/AnasenG4.dir/flags.make
 CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o: ../src/DetectorConstruction.cc
 CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o: CMakeFiles/AnasenG4.dir/compiler_depend.ts
-	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green --progress-dir=/home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles --progress-num=$(CMAKE_PROGRESS_5) "Building CXX object CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o"
+	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green --progress-dir=/home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles --progress-num=$(CMAKE_PROGRESS_6) "Building CXX object CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o"
 	/usr/bin/c++ $(CXX_DEFINES) $(CXX_INCLUDES) $(CXX_FLAGS) -MD -MT CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o -MF CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o.d -o CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o -c /home/jamesszalkie/anasen/Armory/AnasenG4/src/DetectorConstruction.cc
 
 CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.i: cmake_force
@@ -144,6 +158,7 @@ AnasenG4_OBJECTS = \
 "CMakeFiles/AnasenG4.dir/src/main.cc.o" \
 "CMakeFiles/AnasenG4.dir/src/ActionInitialization.cc.o" \
 "CMakeFiles/AnasenG4.dir/src/PrimaryGeneratorAction.cc.o" \
+"CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o" \
 "CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o" \
 "CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o"
 
@@ -153,6 +168,7 @@ AnasenG4_EXTERNAL_OBJECTS =
 AnasenG4: CMakeFiles/AnasenG4.dir/src/main.cc.o
 AnasenG4: CMakeFiles/AnasenG4.dir/src/ActionInitialization.cc.o
 AnasenG4: CMakeFiles/AnasenG4.dir/src/PrimaryGeneratorAction.cc.o
+AnasenG4: CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o
 AnasenG4: CMakeFiles/AnasenG4.dir/src/ReactionGenerator.cc.o
 AnasenG4: CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o
 AnasenG4: CMakeFiles/AnasenG4.dir/build.make
@@ -193,7 +209,7 @@ AnasenG4: /home/jamesszalkie/geant4_installation/geant4-install/lib/libG4global.
 AnasenG4: /home/jamesszalkie/geant4_installation/geant4-install/lib/libG4clhep.so
 AnasenG4: /home/jamesszalkie/geant4_installation/geant4-install/lib/libG4ptl.so.2.3.3
 AnasenG4: CMakeFiles/AnasenG4.dir/link.txt
-	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green --bold --progress-dir=/home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles --progress-num=$(CMAKE_PROGRESS_6) "Linking CXX executable AnasenG4"
+	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --green --bold --progress-dir=/home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles --progress-num=$(CMAKE_PROGRESS_7) "Linking CXX executable AnasenG4"
 	$(CMAKE_COMMAND) -E cmake_link_script CMakeFiles/AnasenG4.dir/link.txt --verbose=$(VERBOSE)
 
 # Rule to build all files generated by this target.

Armory/AnasenG4/build/CMakeFiles/AnasenG4.dir/cmake_clean.cmake

@@ -3,6 +3,8 @@ file(REMOVE_RECURSE
   "AnasenG4.pdb"
   "CMakeFiles/AnasenG4.dir/src/ActionInitialization.cc.o"
   "CMakeFiles/AnasenG4.dir/src/ActionInitialization.cc.o.d"
+  "CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o"
+  "CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o.d"
   "CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o"
   "CMakeFiles/AnasenG4.dir/src/DetectorConstruction.cc.o.d"
   "CMakeFiles/AnasenG4.dir/src/PrimaryGeneratorAction.cc.o"

Armory/AnasenG4/build/CMakeFiles/AnasenG4.dir/progress.make

@@ -4,4 +4,5 @@ CMAKE_PROGRESS_3 = 3
 CMAKE_PROGRESS_4 = 4
 CMAKE_PROGRESS_5 = 5
 CMAKE_PROGRESS_6 = 6
+CMAKE_PROGRESS_7 = 7
 

Armory/AnasenG4/build/CMakeFiles/Makefile2

@@ -81,12 +81,12 @@ clean: CMakeFiles/AnasenG4.dir/clean
 CMakeFiles/AnasenG4.dir/all:
 	$(MAKE) $(MAKESILENT) -f CMakeFiles/AnasenG4.dir/build.make CMakeFiles/AnasenG4.dir/depend
 	$(MAKE) $(MAKESILENT) -f CMakeFiles/AnasenG4.dir/build.make CMakeFiles/AnasenG4.dir/build
-	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --progress-dir=/home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles --progress-num=1,2,3,4,5,6 "Built target AnasenG4"
+	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --progress-dir=/home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles --progress-num=1,2,3,4,5,6,7 "Built target AnasenG4"
 .PHONY : CMakeFiles/AnasenG4.dir/all
 
 # Build rule for subdir invocation for target.
 CMakeFiles/AnasenG4.dir/rule: cmake_check_build_system
-	$(CMAKE_COMMAND) -E cmake_progress_start /home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles 6
+	$(CMAKE_COMMAND) -E cmake_progress_start /home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles 7
 	$(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 CMakeFiles/AnasenG4.dir/all
 	$(CMAKE_COMMAND) -E cmake_progress_start /home/jamesszalkie/anasen/Armory/AnasenG4/build/CMakeFiles 0
 .PHONY : CMakeFiles/AnasenG4.dir/rule

Armory/AnasenG4/build/CMakeFiles/progress.marks

@@ -1 +1 @@
-6
+7

Armory/AnasenG4/build/Makefile

@@ -153,6 +153,30 @@ src/ActionInitialization.cc.s:
 	$(MAKE) $(MAKESILENT) -f CMakeFiles/AnasenG4.dir/build.make CMakeFiles/AnasenG4.dir/src/ActionInitialization.cc.s
 .PHONY : src/ActionInitialization.cc.s
 
+src/BeamReactionSteppingAction.o: src/BeamReactionSteppingAction.cc.o
+.PHONY : src/BeamReactionSteppingAction.o
+
+# target to build an object file
+src/BeamReactionSteppingAction.cc.o:
+	$(MAKE) $(MAKESILENT) -f CMakeFiles/AnasenG4.dir/build.make CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.o
+.PHONY : src/BeamReactionSteppingAction.cc.o
+
+src/BeamReactionSteppingAction.i: src/BeamReactionSteppingAction.cc.i
+.PHONY : src/BeamReactionSteppingAction.i
+
+# target to preprocess a source file
+src/BeamReactionSteppingAction.cc.i:
+	$(MAKE) $(MAKESILENT) -f CMakeFiles/AnasenG4.dir/build.make CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.i
+.PHONY : src/BeamReactionSteppingAction.cc.i
+
+src/BeamReactionSteppingAction.s: src/BeamReactionSteppingAction.cc.s
+.PHONY : src/BeamReactionSteppingAction.s
+
+# target to generate assembly for a file
+src/BeamReactionSteppingAction.cc.s:
+	$(MAKE) $(MAKESILENT) -f CMakeFiles/AnasenG4.dir/build.make CMakeFiles/AnasenG4.dir/src/BeamReactionSteppingAction.cc.s
+.PHONY : src/BeamReactionSteppingAction.cc.s
+
 src/DetectorConstruction.o: src/DetectorConstruction.cc.o
 .PHONY : src/DetectorConstruction.o
 
@@ -261,6 +285,9 @@ help:
 	@echo "... src/ActionInitialization.o"
 	@echo "... src/ActionInitialization.i"
 	@echo "... src/ActionInitialization.s"
+	@echo "... src/BeamReactionSteppingAction.o"
+	@echo "... src/BeamReactionSteppingAction.i"
+	@echo "... src/BeamReactionSteppingAction.s"
 	@echo "... src/DetectorConstruction.o"
 	@echo "... src/DetectorConstruction.i"
 	@echo "... src/DetectorConstruction.s"

Armory/AnasenG4/include/BeamReactionSteppingAction.hh

@@ -0,0 +1,25 @@
+#ifndef BeamReactionSteppingAction_h
+#define BeamReactionSteppingAction_h 1
+
+#include "G4UserSteppingAction.hh"
+#include "G4String.hh"
+
+class G4Step;
+class ReactionGenerator;
+
+class BeamReactionSteppingAction : public G4UserSteppingAction
+{
+  public:
+    BeamReactionSteppingAction();
+    virtual ~BeamReactionSteppingAction();
+
+    virtual void UserSteppingAction(const G4Step* step) override;
+
+  private:
+    ReactionGenerator* fReaction;
+    G4String fTargetVolumeName;
+    G4int fBeamZ;
+    G4int fBeamA;
+};
+
+#endif

Armory/AnasenG4/include/PrimaryGeneratorAction.hh

@@ -6,7 +6,6 @@
 
 class G4ParticleGun;
 class G4Event;
-class ReactionGenerator;
 
 class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
 {
@@ -18,9 +17,9 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
 
   private:
     G4ParticleGun* fParticleGun;
-    ReactionGenerator* fReaction;
+    G4double fBeamRadius;
-    G4double fTargetRadius;
+    G4double fBeamStartZ;
-    G4double fTargetHalfLength;
+    G4double fBeamEnergyPerU;
 };
 
 #endif

Armory/AnasenG4/include/ReactionGenerator.hh

@@ -1,7 +1,15 @@
 #ifndef ReactionGenerator_h
 #define ReactionGenerator_h 1
 
+#include "G4ThreeVector.hh"
+
+enum class ReactionChannel {
+  Proton,
+  Alpha
+};
+
 struct ReactionOutput {
+  ReactionChannel channel;
   double kineticEnergy;
   double px;
   double py;
@@ -13,11 +21,13 @@ class ReactionGenerator {
   ReactionGenerator();
   ~ReactionGenerator();
 
-  ReactionOutput SampleEvent();
+  ReactionOutput SampleEvent(double beamKineticEnergy, const G4ThreeVector& beamDirection);
+  double GetEffectiveCrossSection() const;
 
  private:
   class Impl;
   Impl* fImpl;
+  double fEffectiveSigma;
 };
 
 #endif

Armory/AnasenG4/run.mac

@@ -3,4 +3,10 @@
 /event/verbose 0
 /tracking/verbose 0
 /run/initialize
-/run/beamOn 10000
+
+/vis/filtering/trajectories/create/particleFilter
+/vis/filtering/trajectories/particleFilter-0/add e-
+/vis/filtering/trajectories/particleFilter-0/add e+
+/vis/filtering/trajectories/particleFilter-0/invert true
+
+/run/beamOn 100

Armory/AnasenG4/src/ActionInitialization.cc

@@ -1,5 +1,6 @@
 #include "ActionInitialization.hh"
 #include "PrimaryGeneratorAction.hh"
+#include "BeamReactionSteppingAction.hh"
 
 ActionInitialization::ActionInitialization()
 : G4VUserActionInitialization()
@@ -10,4 +11,5 @@ ActionInitialization::~ActionInitialization() {}
 void ActionInitialization::Build() const
 {
   SetUserAction(new PrimaryGeneratorAction());
+  SetUserAction(new BeamReactionSteppingAction());
 }

Armory/AnasenG4/src/BeamReactionSteppingAction.cc

@@ -0,0 +1,103 @@
+#include "BeamReactionSteppingAction.hh"
+
+#include "ReactionGenerator.hh"
+
+#include "G4Step.hh"
+#include "G4Track.hh"
+#include "G4DynamicParticle.hh"
+#include "G4EventManager.hh"
+#include "G4StackManager.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4Proton.hh"
+#include "G4Alpha.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4RandomTools.hh"
+
+#include <cmath>
+
+BeamReactionSteppingAction::BeamReactionSteppingAction()
+: G4UserSteppingAction(),
+  fReaction(new ReactionGenerator()),
+  fTargetVolumeName("Target"),
+  fBeamZ(13),
+  fBeamA(27)
+{}
+
+BeamReactionSteppingAction::~BeamReactionSteppingAction()
+{
+  delete fReaction;
+}
+
+void BeamReactionSteppingAction::UserSteppingAction(const G4Step* step)
+{
+  const G4Track* track = step->GetTrack();
+  if (track->GetTrackStatus() != fAlive) {
+    return;
+  }
+
+  const G4ParticleDefinition* particle = track->GetDefinition();
+  if (particle->GetParticleType() != "nucleus") {
+    return;
+  }
+
+  if (particle->GetAtomicNumber() != fBeamZ || particle->GetAtomicMass() != fBeamA) {
+    return;
+  }
+
+  const auto* prePoint = step->GetPreStepPoint();
+  if (!prePoint || !prePoint->GetPhysicalVolume()) {
+    return;
+  }
+
+  if (prePoint->GetPhysicalVolume()->GetName() != fTargetVolumeName) {
+    return;
+  }
+
+  const G4double stepLength = step->GetStepLength();
+  if (stepLength <= 0.0) {
+    return;
+  }
+
+  const auto* material = prePoint->GetMaterial();
+  if (!material) {
+    return;
+  }
+
+  const G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
+  const G4double sigma = fReaction->GetEffectiveCrossSection();
+  const G4double probability = 1.0 - std::exp(-atomDensity * sigma * stepLength);
+
+  if (G4UniformRand() > probability) {
+    return;
+  }
+
+  ReactionOutput output = fReaction->SampleEvent(track->GetKineticEnergy(), track->GetMomentumDirection());
+  if (output.kineticEnergy <= 0.0) {
+    return;
+  }
+
+  const G4ThreeVector momentum(output.px, output.py, output.pz);
+  G4ThreeVector direction = momentum;
+  if (direction.mag2() > 0.0) {
+    direction = direction.unit();
+  } else {
+    direction = G4ThreeVector(0.0, 0.0, 1.0);
+  }
+
+  G4ParticleDefinition* product = nullptr;
+  if (output.channel == ReactionChannel::Proton) {
+    product = G4Proton::ProtonDefinition();
+  } else {
+    product = G4Alpha::AlphaDefinition();
+  }
+
+  auto* dynamicParticle = new G4DynamicParticle(product, direction, output.kineticEnergy);
+  auto* secondaryTrack = new G4Track(dynamicParticle, track->GetGlobalTime(), prePoint->GetPosition());
+  secondaryTrack->SetParentID(track->GetTrackID());
+  secondaryTrack->SetTouchableHandle(track->GetTouchableHandle());
+
+  G4EventManager::GetEventManager()->GetStackManager()->PushOneTrack(secondaryTrack, nullptr);
+
+  G4Track* mutableTrack = const_cast<G4Track*>(track);
+  mutableTrack->SetTrackStatus(fStopAndKill);
+}

Armory/AnasenG4/src/DetectorConstruction.cc

@@ -9,6 +9,7 @@
 #include "G4SystemOfUnits.hh"
 #include "G4RotationMatrix.hh"
 #include "G4ThreeVector.hh"
+#include "G4VisAttributes.hh"
 
 DetectorConstruction::DetectorConstruction()
 : G4VUserDetectorConstruction(),
@@ -36,12 +37,29 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
   G4VPhysicalVolume* physWorld = new G4PVPlacement(0, G4ThreeVector(), logicWorld, "World", 0, false, 0, true);
 
   // Helium target
+
+  G4double pressure = 379.0 * 133.322368 * pascal;
+  G4double temperature = 293.15 * kelvin;
+  G4double density = 8.29e-5 * g/cm3;
+
+  G4Material* HeGas =
+      new G4Material("HeGas",
+                    density,
+                    1,
+                    kStateGas,
+                    temperature,
+                    pressure);
+
+  HeGas->AddElement(nist->FindOrBuildElement("He"), 1);
   G4Material* he = nist->FindOrBuildMaterial("G4_He");
-  const G4double targetRadius = 20.*mm;
+
+  const G4double targetRadius = 200.*mm;
   const G4double targetHalfLength = 175.*mm;
   G4Tubs* solidTarget = new G4Tubs("Target", 0, targetRadius, targetHalfLength, 0, 360*deg);
-  G4LogicalVolume* logicTarget = new G4LogicalVolume(solidTarget, he, "Target");
+  G4LogicalVolume* logicTarget =
+    new G4LogicalVolume(solidTarget, HeGas, "Target");
   new G4PVPlacement(0, G4ThreeVector(), logicTarget, "Target", logicWorld, false, 0, true);
+  logicTarget->SetVisAttributes(G4VisAttributes::GetInvisible());
   fScoringVolume = logicTarget;
 
   // Axes (optional, for visualization)
@@ -63,7 +81,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
 
   // ANASEN geometry
   const int nWire = 24;
-  const int wireShift = 3;
+  const int wireShift = 4;
   const G4double zLen = 350.*mm;
   const G4double radiusA = 38.*mm;
   const G4double radiusC = 43.*mm;
@@ -73,7 +91,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
   G4double wireALength = sqrt(zLen*zLen + pow(2 * radiusA * sin(dAngle/2.), 2));
   G4double wireATheta = atan2(2 * radiusA * sin(dAngle / 2.), zLen);
 
-  G4Tubs* solidPC_A = new G4Tubs("PC_A", 0, 1*mm, wireALength/2., 0, 360*deg);
+  G4Tubs* solidPC_A = new G4Tubs("PC_A", 0, .1*mm, wireALength/2., 0, 360*deg);
   G4LogicalVolume* logicPC_A = new G4LogicalVolume(solidPC_A, al, "PC_A");
 
   for(int i = 0; i < nWire; i++){
@@ -91,7 +109,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
   G4double wireCLength = sqrt(zLen*zLen + pow(2 * radiusC * sin(dAngle/2.), 2));
   G4double wireCTheta = atan2(2 * radiusC * sin(dAngle / 2.), zLen);
 
-  G4Tubs* solidPC_C = new G4Tubs("PC_C", 0, 1*mm, wireCLength/2., 0, 360*deg);
+  G4Tubs* solidPC_C = new G4Tubs("PC_C", 0, .1*mm, wireCLength/2., 0, 360*deg);
   G4LogicalVolume* logicPC_C = new G4LogicalVolume(solidPC_C, al, "PC_C");
 
   for(int i = 0; i < nWire; i++){

Armory/AnasenG4/src/PrimaryGeneratorAction.cc

@@ -1,6 +1,7 @@
 #include "PrimaryGeneratorAction.hh"
 #include "G4ParticleGun.hh"
 #include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
 #include "G4Proton.hh"
 #include "G4SystemOfUnits.hh"
 #include "G4Event.hh"
@@ -8,54 +9,42 @@
 #include "G4RandomTools.hh"
 #include "G4UnitsTable.hh"
 
-#include "ReactionGenerator.hh"
-
 #include <cmath>
 
 PrimaryGeneratorAction::PrimaryGeneratorAction()
 : G4VUserPrimaryGeneratorAction(),
   fParticleGun(nullptr),
-  fReaction(nullptr),
+  fBeamRadius(2.0 * mm),
-  fTargetRadius(20.0 * mm),
+  fBeamStartZ(-190.0 * mm),
-  fTargetHalfLength(175.0 * mm)
+  fBeamEnergyPerU(10.0 * MeV)
 {
   G4int nParticle = 1;
   fParticleGun = new G4ParticleGun(nParticle);
   fParticleGun->SetParticleDefinition(G4Proton::ProtonDefinition());
   fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
-  fParticleGun->SetParticleEnergy(1.0 * MeV);
+  fParticleGun->SetParticleEnergy(27.0 * fBeamEnergyPerU);
-
-  fReaction = new ReactionGenerator();
 }
 
 PrimaryGeneratorAction::~PrimaryGeneratorAction()
 {
   delete fParticleGun;
-  delete fReaction;
 }
 
 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
 {
-  ReactionOutput reaction = fReaction->SampleEvent();
+  auto* ion = G4IonTable::GetIonTable()->GetIon(13, 27, 0.0);
-
+  if (ion) {
-  G4double kineticEnergy = reaction.kineticEnergy;
+    fParticleGun->SetParticleDefinition(ion);
-  if (kineticEnergy < 0) kineticEnergy = 0;
-
-  G4ThreeVector direction(reaction.px, reaction.py, reaction.pz);
-  if (direction.mag() > 0) {
-    direction = direction.unit();
-  } else {
-    direction = G4ThreeVector(0., 0., 1.);
   }
 
-  G4double r = fTargetRadius * std::sqrt(G4UniformRand());
+  G4double r = fBeamRadius * std::sqrt(G4UniformRand());
   G4double phi = twopi * G4UniformRand();
   G4double x = r * std::cos(phi);
   G4double y = r * std::sin(phi);
-  G4double z = (2.0 * G4UniformRand() - 1.0) * fTargetHalfLength;
+  G4double z = fBeamStartZ;
 
   fParticleGun->SetParticlePosition(G4ThreeVector(x, y, z));
-  fParticleGun->SetParticleMomentumDirection(direction);
+  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
-  fParticleGun->SetParticleEnergy(kineticEnergy);
+  fParticleGun->SetParticleEnergy(27.0 * fBeamEnergyPerU);
   fParticleGun->GeneratePrimaryVertex(event);
 }

Armory/AnasenG4/src/ReactionGenerator.cc

@@ -7,21 +7,29 @@
 #include <cmath>
 
 struct ReactionGenerator::Impl {
-  TransferReaction transfer;
+  TransferReaction transferP;
+  TransferReaction transferAlpha;
 };
 
 ReactionGenerator::ReactionGenerator()
-: fImpl(nullptr)
+: fImpl(nullptr),
+  fEffectiveSigma(5e3 * barn) //originally 5e6
 {
   // Ensure the isotope mass table is found when running from the AnasenG4 folder.
   massData = "/home/jamesszalkie/anasen/Armory/mass20.txt";
   fImpl = new Impl();
-  fImpl->transfer.SetA(18, 10, 0);
+
-  fImpl->transfer.Seta(4, 2);
+  // p-channel: 27Al + 4He -> p + 30Si
-  fImpl->transfer.Setb(1, 1);
+  fImpl->transferP.SetA(27, 13, 0);
-  fImpl->transfer.SetB(21, 11, 0);
+  fImpl->transferP.Seta(4, 2);
-  fImpl->transfer.SetIncidentEnergyAngle(10.0, 0.0, 0.0);
+  fImpl->transferP.Setb(1, 1);
-  fImpl->transfer.CalReactionConstant();
+  fImpl->transferP.SetB(30, 14, 0);
+
+  // alpha-channel: 27Al + 4He -> alpha + 27Al
+  fImpl->transferAlpha.SetA(27, 13, 0);
+  fImpl->transferAlpha.Seta(4, 2);
+  fImpl->transferAlpha.Setb(4, 2);
+  fImpl->transferAlpha.SetB(27, 13, 0);
 }
 
 ReactionGenerator::~ReactionGenerator()
@@ -29,20 +37,46 @@ ReactionGenerator::~ReactionGenerator()
   delete fImpl;
 }
 
-ReactionOutput ReactionGenerator::SampleEvent()
+ReactionOutput ReactionGenerator::SampleEvent(double beamKineticEnergy, const G4ThreeVector& beamDirection)
 {
-  const double cosTheta = 2.0 * gRandom->Rndm() - 1.0;
+  // Forward-peaked sampling gives the expected inverse E(theta)-like trend in lab observables.
-  const double thetaCM = std::acos(cosTheta);
+  double thetaCM = 0.0;
+  while (true) {
+    const double trial = gRandom->Rndm() * M_PI;
+    const double weight = 1.0 / (1.0 + 20.0 * trial * trial);
+    if (gRandom->Rndm() < weight) {
+      thetaCM = trial;
+      break;
+    }
+  }
   const double phiCM = gRandom->Rndm() * 2.0 * M_PI;
 
-  TLorentzVector* reactionOutput = fImpl->transfer.Event(thetaCM, phiCM);
+  const double beamEnergyPerU = beamKineticEnergy / (27.0 * MeV);
-  TLorentzVector proton = reactionOutput[2];
+  const double beamThetaDeg = beamDirection.theta() / deg;
+  const double beamPhiDeg = beamDirection.phi() / deg;
+
+  // Use a mild energy-dependent branching to mix p and alpha events.
+  const double protonFraction = (beamEnergyPerU > 8.0) ? 0.65 : 0.45;
+  const bool isProtonChannel = (gRandom->Rndm() < protonFraction);
+
+  TransferReaction* transfer = isProtonChannel ? &fImpl->transferP : &fImpl->transferAlpha;
+  transfer->SetIncidentEnergyAngle(beamEnergyPerU, beamThetaDeg, beamPhiDeg);
+  transfer->CalReactionConstant();
+
+  TLorentzVector* reactionOutput = transfer->Event(thetaCM, phiCM);
+  TLorentzVector light = reactionOutput[2];
   delete[] reactionOutput;
 
   ReactionOutput output;
-  output.kineticEnergy = proton.E() - proton.M();
+  output.channel = isProtonChannel ? ReactionChannel::Proton : ReactionChannel::Alpha;
-  output.px = proton.Px();
+  output.kineticEnergy = light.E() - light.M();
-  output.py = proton.Py();
+  output.px = light.Px();
-  output.pz = proton.Pz();
+  output.py = light.Py();
+  output.pz = light.Pz();
   return output;
 }
+
+double ReactionGenerator::GetEffectiveCrossSection() const
+{
+  return fEffectiveSigma;
+}

Armory/AnasenG4/src/main.cc

@@ -14,7 +14,7 @@ int main(int argc, char** argv)
     ui = new G4UIExecutive(argc, argv);
   }
 
-  auto* runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default);
+  auto* runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::SerialOnly);
 
   runManager->SetUserInitialization(new DetectorConstruction());
 

Armory/AnasenG4/vis.mac

@@ -1,4 +1,3 @@
-# Interactive visualization setup modeled on Geant4 basic example B1
 
 /vis/open
 /vis/viewer/set/autoRefresh false
@@ -22,7 +21,7 @@
 /vis/scene/endOfEventAction accumulate 20
 /vis/geometry/set/visibility World 0 false
 /vis/scene/add/axes 0 0 0 20 mm
-/vis/viewer/set/background 1 1 1
+/vis/viewer/set/background 0 0 0
 /vis/viewer/set/style surface
 /vis/viewer/set/hiddenMarker true
 /vis/viewer/set/viewpointThetaPhi 120 150