nam/muon-tomograpy
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

initial working version: 80 scintillators in 8 layers, root output, mono-energetic muon on -Z direction

Browse Source
This commit is contained in:
Nam Tran
2021-01-14 16:29:21 -06:00
commit 6a5eab0052
21 changed files with 915 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored Normal file
View File

@@ -0,0 +1,2 @@
build

68
CMakeLists.txt Normal file
View File

@@ -0,0 +1,68 @@
#----------------------------------------------------------------------------
# Setup the project
cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
project(muograpy)
#----------------------------------------------------------------------------
# Find Geant4 package, activating all available UI and Vis drivers by default
# You can set WITH_GEANT4_UIVIS to OFF via the command line or ccmake/cmake-gui
# to build a batch mode only executable
#
option(WITH_GEANT4_UIVIS "Build example with Geant4 UI and Vis drivers" ON)
if(WITH_GEANT4_UIVIS)
find_package(Geant4 REQUIRED ui_all vis_all)
else()
find_package(Geant4 REQUIRED)
endif()
#----------------------------------------------------------------------------
# Setup Geant4 include directories and compile definitions
# Setup include directory for this project
#
include(${Geant4_USE_FILE})
include_directories(${PROJECT_SOURCE_DIR}/include)
#----------------------------------------------------------------------------
# Locate sources and headers for this project
# NB: headers are included so they will show up in IDEs
#
file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cc)
file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh)
#----------------------------------------------------------------------------
# Add the executable, and link it to the Geant4 libraries
#
add_executable(muography muography.cc ${sources} ${headers})
target_link_libraries(muography ${Geant4_LIBRARIES})
#----------------------------------------------------------------------------
# Copy all scripts to the build directory, i.e. the directory in which we
# build B1. This is so that we can run the executable directly because it
# relies on these scripts being in the current working directory.
#
FILE(GLOB RUN_MACROS macros/*.mac)
foreach(_script ${RUN_MACROS})
# message( ${_script} )
get_filename_component(shortname ${_script} NAME)
# message(${shortname})
configure_file(
macros/${shortname}
${PROJECT_BINARY_DIR}/${shortname}
COPYONLY
)
endforeach()
#----------------------------------------------------------------------------
# For internal Geant4 use - but has no effect if you build this
# example standalone
#
# add_custom_target(muography DEPENDS muography)
#----------------------------------------------------------------------------
# Install the executable to 'bin' directory under CMAKE_INSTALL_PREFIX
#
install(TARGETS muography DESTINATION bin)

23
include/ActionInitialization.hh Normal file
View File

@@ -0,0 +1,23 @@
#ifndef ActionInitialization_h
#define ActionInitialization_h 1
#include "G4VUserActionInitialization.hh"
/// Action initialization class.
class ActionInitialization : public G4VUserActionInitialization
{
public:
ActionInitialization();
virtual ~ActionInitialization();
virtual void BuildForMaster() const;
virtual void Build() const;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

6
include/Analysis.hh Normal file
View File

@@ -0,0 +1,6 @@
#ifndef Analysis_h
#define Analysis_h
#include "g4root.hh"
#endif /* end of include guard */

30
include/DetectorConstruction.hh Normal file
View File

@@ -0,0 +1,30 @@
#ifndef DetectorConstruction_h
#define DetectorConstruction_h 1
#include "G4VUserDetectorConstruction.hh"
#include "globals.hh"
class G4VPhysicalVolume;
class G4LogicalVolume;
/// Detector construction class to define materials and geometry.
class DetectorConstruction : public G4VUserDetectorConstruction
{
public:
DetectorConstruction();
virtual ~DetectorConstruction();
virtual G4VPhysicalVolume* Construct();
virtual void ConstructSDandField();
protected:
// void DefineMaterials();
G4bool fCheckOverlaps;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

32
include/EventAction.hh Normal file
View File

@@ -0,0 +1,32 @@
#ifndef EventAction_h
#define EventAction_h 1
#include "G4UserEventAction.hh"
#include "globals.hh"
class RunAction;
/// Event action class
///
class EventAction : public G4UserEventAction
{
public:
EventAction(RunAction* runAction);
virtual ~EventAction();
virtual void BeginOfEventAction(const G4Event* event);
virtual void EndOfEventAction(const G4Event* event);
private:
// RunAction *fRunAction;
G4int fCollID_detX;
G4int fCollID_detY;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

37
include/PrimaryGeneratorAction.hh Normal file
View File

@@ -0,0 +1,37 @@
#ifndef PrimaryGeneratorAction_h
#define PrimaryGeneratorAction_h 1
#include "G4VUserPrimaryGeneratorAction.hh"
#include "G4ParticleGun.hh"
#include "globals.hh"
class G4ParticleGun;
class G4Event;
class G4Box;
/// The primary generator action class with particle gun.
///
/// The default kinematic is a 6 MeV gamma, randomly distribued
/// in front of the phantom across 80% of the (X,Y) phantom size.
class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
{
public:
PrimaryGeneratorAction();
virtual ~PrimaryGeneratorAction();
// method from the base class
virtual void GeneratePrimaries(G4Event*);
// method to access particle gun
const G4ParticleGun* GetParticleGun() const { return fParticleGun; }
private:
G4ParticleGun* fParticleGun; // pointer a to G4 gun class
G4Box* fEnvelopeBox;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

31
include/RunAction.hh Normal file
View File

@@ -0,0 +1,31 @@
#ifndef RunAction_h
#define RunAction_h 1
#include "G4UserRunAction.hh"
#include "G4Accumulable.hh"
#include "globals.hh"
class G4Run;
/// Run action class
///
/// In EndOfRunAction(), it calculates the dose in the selected volume
/// from the energy deposit accumulated via stepping and event actions.
/// The computed dose is then printed on the screen.
class RunAction : public G4UserRunAction
{
public:
RunAction();
virtual ~RunAction();
// virtual G4Run* GenerateRun();
virtual void BeginOfRunAction(const G4Run*);
virtual void EndOfRunAction(const G4Run*);
private:
};
#endif

3
macros/batch_10k.mac Normal file
View File

@@ -0,0 +1,3 @@
/run/initialize
/tracking/verbose 0
/run/beamOn 10000

3
macros/batch_1M.mac Normal file
View File

@@ -0,0 +1,3 @@
/run/initialize
/tracking/verbose 0
/run/beamOn 1000000

13
macros/init_vis.mac Normal file
View File

@@ -0,0 +1,13 @@
# Set some default verbose
/control/verbose 2
/control/saveHistory
/run/verbose 2
#
# Change the default number of threads (in multi-threaded mode)
#/run/numberOfThreads 4
#
# Initialize kernel
/run/initialize
#
# Visualization setting
/control/execute vis.mac

23
macros/run1.mac Normal file
View File

@@ -0,0 +1,23 @@
#
# Initialize kernel
/run/initialize
#
/control/verbose 2
/run/verbose 2
/event/verbose 0
/tracking/verbose 1
#
# gamma 6 MeV to the direction (0.,0.,1.)
#
/gun/particle gamma
/gun/energy 6 MeV
#
/run/beamOn 5
#
# proton 210 MeV to the direction (0.,0.,1.)
#
/gun/particle proton
/gun/energy 210 MeV
/tracking/verbose 2
#
/run/beamOn 1

22
macros/run2.mac Normal file
View File

@@ -0,0 +1,22 @@
#/run/numberOfWorkers 4
/run/initialize
#
/control/verbose 2
/run/verbose 2
#
# gamma 6 MeV to the direction (0.,0.,1.)
# 10000 events
#
/gun/particle gamma
/gun/energy 6 MeV
#
/run/printProgress 100
/run/beamOn 1000
#
# proton 210 MeV to the direction (0.,0.,1.)
# 1000 events
#
/gun/particle proton
/gun/energy 210 MeV
#
/run/beamOn 1000

17
macros/steal.sh Executable file
View File

@@ -0,0 +1,17 @@
#!/bin/bash
prefix=B1
for source in $(ls include/*.hh); do
newSource=${source/$prefix/} # replace all occurence of prefix with zero character
echo $source "-->" $newSource
mv $source $newSource
sed -i '' "s/$prefix//g" $newSource
sed -i '' -e '1,28d' $newSource
done
for source in $(ls src/*.cc); do
newSource=${source/$prefix/} # replace all occurence of prefix with zero character
echo $source "-->" $newSource
mv $source $newSource
sed -i '' "s/$prefix//g" $newSource
sed -i '' -e '1,28d' $newSource
done

101
macros/vis.mac Normal file
View File

@@ -0,0 +1,101 @@
# Use these open statements to open selected visualization
#
# Use this open statement to create an OpenGL view:
/vis/open OGL 600x600-0+0
#
# Use this open statement to create an OpenInventor view:
#/vis/open OI
#
# Use this open statement to create a .prim file suitable for
# viewing in DAWN:
#/vis/open DAWNFILE
#
# Use this open statement to create a .heprep file suitable for
# viewing in HepRApp:
#/vis/open HepRepFile
#
# Use this open statement to create a .wrl file suitable for
# viewing in a VRML viewer:
#/vis/open VRML2FILE
#
# Disable auto refresh and quieten vis messages whilst scene and
# trajectories are established:
/vis/viewer/set/autoRefresh false
/vis/verbose errors
#
# Draw geometry:
/vis/drawVolume
#
# Specify view angle:
/vis/viewer/set/viewpointVector -1 0 0
/vis/viewer/set/lightsVector -1 0 0
#
# Specify style (surface, wireframe, auxiliary edges,...)
/vis/viewer/set/style wireframe
/vis/viewer/set/auxiliaryEdge true
/vis/viewer/set/lineSegmentsPerCircle 100
#
# Draw smooth trajectories at end of event, showing trajectory points
# as markers 2 pixels wide:
/vis/scene/add/trajectories smooth
/vis/modeling/trajectories/create/drawByCharge
/vis/modeling/trajectories/drawByCharge-0/default/setDrawStepPts true
/vis/modeling/trajectories/drawByCharge-0/default/setStepPtsSize 2
# (if too many tracks cause core dump => /tracking/storeTrajectory 0)
#
# Draw hits at end of event:
#/vis/scene/add/hits
#
# To draw only gammas:
#/vis/filtering/trajectories/create/particleFilter
#/vis/filtering/trajectories/particleFilter-0/add gamma
#
# To invert the above, drawing all particles except gammas,
# keep the above two lines but also add:
#/vis/filtering/trajectories/particleFilter-0/invert true
#
# Many other options are available with /vis/modeling and /vis/filtering.
# For example, to select colour by particle ID:
#/vis/modeling/trajectories/create/drawByParticleID
#/vis/modeling/trajectories/drawByParticleID-0/default/setDrawStepPts true
# To select or override default colours (note: e+ is blue by default):
#/vis/modeling/trajectories/list
#/vis/modeling/trajectories/drawByParticleID-0/set e+ yellow
#
# To superimpose all of the events from a given run:
/vis/scene/endOfEventAction accumulate
#
# Decorations
# Name
/vis/set/textColour green
/vis/set/textLayout right
/vis/set/textLayout # Revert to normal (left adjusted) layout
/vis/set/textColour # Revert to default text colour (blue)
#
# Axes, scale, etc.
# /vis/scene/add/scale # Simple scale line
/vis/scene/add/axes # Simple axes: x=red, y=green, z=blue.
/vis/scene/add/eventID # Drawn at end of event
#
# Frame
/vis/set/colour red
/vis/set/lineWidth 2
/vis/scene/add/frame # Simple frame around the view
/vis/set/colour # Revert to default colour (white)
/vis/set/lineWidth # Revert to default line width (1.)
#
# To get nice view
# Make the "World" box invisible
/vis/geometry/set/visibility World 0 false
# "Envelope" is transparent blue to represent water
# /vis/geometry/set/colour Envelope 0 0 0 1 .3
/vis/viewer/set/style surface
/vis/viewer/set/hiddenMarker true
/vis/viewer/set/viewpointThetaPhi 120 150
#
# Re-establish auto refreshing and verbosity:
/vis/viewer/set/autoRefresh true
/vis/verbose warnings
#
# For file-based drivers, use this to create an empty detector view:
#/vis/viewer/flush

117
muography.cc Normal file
View File

@@ -0,0 +1,117 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
/// \file example.cc
/// \brief Main program of the example
#include "DetectorConstruction.hh"
#include "ActionInitialization.hh"
#ifdef G4MULTITHREADED
#include "G4MTRunManager.hh"
#else
#include "G4RunManager.hh"
#endif
#include "G4UImanager.hh"
#include "QBBC.hh"
#include "G4VisExecutive.hh"
#include "G4UIExecutive.hh"
#include "Randomize.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
int main(int argc,char** argv)
{
// Detect interactive mode (if no arguments) and define UI session
//
G4UIExecutive* ui = 0;
if ( argc == 1 ) {
ui = new G4UIExecutive(argc, argv);
}
// Choose the Random engine
G4Random::setTheEngine(new CLHEP::RanecuEngine);
// Construct the default run manager
//
#ifdef G4MULTITHREADED
G4MTRunManager* runManager = new G4MTRunManager;
runManager->SetNumberOfThreads(4);
#else
G4RunManager* runManager = new G4RunManager;
#endif
// Set mandatory initialization classes
//
// Detector construction
runManager->SetUserInitialization(new DetectorConstruction());
// Physics list
G4VModularPhysicsList* physicsList = new QBBC;
physicsList->SetVerboseLevel(0);
runManager->SetUserInitialization(physicsList);
// User action initialization
runManager->SetUserInitialization(new ActionInitialization());
// Initialize visualization
//
G4VisManager* visManager = new G4VisExecutive;
// G4VisExecutive can take a verbosity argument - see /vis/verbose guidance.
// G4VisManager* visManager = new G4VisExecutive("Quiet");
visManager->Initialize();
// Get the pointer to the User Interface manager
G4UImanager* UImanager = G4UImanager::GetUIpointer();
// Process macro or start UI session
//
if ( ! ui ) {
// batch mode
G4String command = "/control/execute ";
G4String fileName = argv[1];
UImanager->ApplyCommand(command+fileName);
}
else {
// interactive mode
UImanager->ApplyCommand("/control/execute init_vis.mac");
ui->SessionStart();
delete ui;
}
// Job termination
// Free the store: user actions, physics_list and detector_description are
// owned and deleted by the run manager, so they should not be deleted
// in the main() program !
delete visManager;
delete runManager;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

35
src/ActionInitialization.cc Normal file
View File

@@ -0,0 +1,35 @@
#include "ActionInitialization.hh"
#include "EventAction.hh"
#include "PrimaryGeneratorAction.hh"
#include "RunAction.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ActionInitialization::ActionInitialization() : G4VUserActionInitialization() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ActionInitialization::~ActionInitialization() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ActionInitialization::BuildForMaster() const {
RunAction *runAction = new RunAction;
SetUserAction(runAction);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ActionInitialization::Build() const {
SetUserAction(new PrimaryGeneratorAction);
RunAction *runAction = new RunAction;
SetUserAction(runAction);
EventAction *eventAction = new EventAction(runAction);
SetUserAction(eventAction);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

157
src/DetectorConstruction.cc Normal file
View File

@@ -0,0 +1,157 @@
#include "DetectorConstruction.hh"
#include "G4Box.hh"
#include "G4LogicalVolume.hh"
#include "G4MultiFunctionalDetector.hh"
#include "G4NistManager.hh"
#include "G4PSDoseDeposit.hh"
#include "G4PSEnergyDeposit.hh"
#include "G4PVPlacement.hh"
#include "G4PhysicalConstants.hh"
#include "G4RotationMatrix.hh"
#include "G4RunManager.hh"
#include "G4SDManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4Transform3D.hh"
#include "G4Tubs.hh"
#include "G4VPrimitiveScorer.hh"
#include "G4VisAttributes.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
DetectorConstruction::DetectorConstruction()
: G4VUserDetectorConstruction(), fCheckOverlaps(true) {
// DefineMaterials();
}
DetectorConstruction::~DetectorConstruction() {}
G4VPhysicalVolume *DetectorConstruction::Construct() {
G4NistManager *nist = G4NistManager::Instance();
G4Material *default_mat = nist->FindOrBuildMaterial("G4_AIR");
// G4Material *target_mat = nist->FindOrBuildMaterial("G4_Cu");
G4Material *scMaterial =
nist->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
// World and crystals
G4double scSizeX = 5.0 * cm;
G4double scSizeY = 50.0 * cm;
G4double scSizeZ = 1. * cm;
G4double scGapX = 0.2 * cm;
G4double scGapZ = 0.2 * cm;
G4int nScBar = 10;
G4int nLayerTop = 2;
G4int nLayerBottom = 2;
G4double layerGapZ = 1 * cm;
G4double topBotGap = 50 * cm;
G4double layerSizeX = nScBar * scSizeX + nScBar * scGapX;
G4double layerSizeY = nScBar * scSizeX + nScBar * scGapX;
G4double layerSizeZ = scSizeZ + scGapZ;
G4double worldSizeX = 1.2 * layerSizeX;
G4double worldSizeY = 1.2 * layerSizeY;
G4double worldSizeZ =
1.2 * ((nLayerBottom + nLayerTop) * (layerSizeZ + layerGapZ) + topBotGap);
// world volume
G4Box *solidWorld =
new G4Box("World", 0.5 * worldSizeX, 0.5 * worldSizeY, 0.5 * worldSizeZ);
G4LogicalVolume *logicWorld =
new G4LogicalVolume(solidWorld, default_mat, "World");
G4VPhysicalVolume *physWorld =
new G4PVPlacement(0, // no rotation
G4ThreeVector(), // at (0,0,0)
logicWorld, // its logical volume
"World", // its name
0, // its mother volume
false, // no boolean operation
0, // copy number
fCheckOverlaps); // overlaps checking
// scintillators
G4Box *solidScX =
new G4Box("solidScX", 0.5 * scSizeX, 0.5 * scSizeY, 0.5 * scSizeZ);
G4LogicalVolume *logicScX =
new G4LogicalVolume(solidScX, scMaterial, "logicScX");
G4Box *solidScY =
new G4Box("solidScY", 0.5 * scSizeY, 0.5 * scSizeX, 0.5 * scSizeZ);
G4LogicalVolume *logicScY =
new G4LogicalVolume(solidScY, scMaterial, "logicScY");
// G4Box *solidLayer =
// new G4Box("Layer", 0.5 * layerSizeX, 0.5 * layerSizeY, 0.5 * layerSizeZ);
// G4LogicalVolume *logicLayer =
// new G4LogicalVolume(solidLayer, default_mat, "logicLayer");
G4double scOffset = 0. * cm;
G4double copyNo = 0;
G4double layerZPosition[4] = {
-topBotGap / 2 - 2 * layerSizeZ - 2 * layerGapZ, // layer 0
-topBotGap / 2, // layer 1
topBotGap / 2, // layer 2
topBotGap / 2 + 2 * layerSizeZ + 2 * layerGapZ}; // layer 3
G4RotationMatrix *rotZ90 = new G4RotationMatrix();
rotZ90->rotateZ(90 * deg);
// place solid Scintillator into X and Y layers
G4double Zpos = 0.;
for (int iLayer = 0; iLayer < nLayerTop + nLayerBottom; iLayer++) {
for (int iBar = 0; iBar < nScBar; iBar++) {
scOffset = iBar * (scGapX + scSizeX) - layerSizeX / 2 + scSizeX / 2;
copyNo = iBar + nScBar * iLayer;
Zpos = layerZPosition[iLayer];
new G4PVPlacement(0, // rotation
G4ThreeVector(scOffset, 0, Zpos), // offset
logicScX, // logical volume
"ScX", // name
logicWorld, // mother volume
false, // no boolean operation
copyNo, // copy number
fCheckOverlaps); // overlaps checking
Zpos = layerZPosition[iLayer] + layerSizeZ + layerGapZ;
new G4PVPlacement(0, // rotation
G4ThreeVector(0, scOffset, Zpos), // offset
logicScY, // logical volume
"ScY", // name
logicWorld, // mother volume
false, // no boolean operation
copyNo, // copy number
fCheckOverlaps); // overlaps checking
}
}
//------------------------------------------------
// Visualization attributes
//------------------------------------------------
logicWorld->SetVisAttributes(G4VisAttributes::Invisible);
G4VisAttributes *visAttScX = new G4VisAttributes(G4Colour(.0, 1.0, 1., 1.0));
logicScX->SetVisAttributes(visAttScX);
G4VisAttributes *visAttScY = new G4VisAttributes(G4Colour(0.0, 0.0, 1., 0.4));
logicScY->SetVisAttributes(visAttScY);
// always return the physical World
//
return physWorld;
}
void DetectorConstruction::ConstructSDandField() {
G4SDManager::GetSDMpointer()->SetVerboseLevel(0);
G4MultiFunctionalDetector *mfDetX = new G4MultiFunctionalDetector("detX");
G4SDManager::GetSDMpointer()->AddNewDetector(mfDetX);
G4VPrimitiveScorer *primitiv1 = new G4PSEnergyDeposit("edep");
mfDetX->RegisterPrimitive(primitiv1);
SetSensitiveDetector("logicScX", mfDetX);
G4MultiFunctionalDetector *mfDetY = new G4MultiFunctionalDetector("detY");
G4SDManager::GetSDMpointer()->AddNewDetector(mfDetY);
G4VPrimitiveScorer *primitiv2 = new G4PSEnergyDeposit("edep");
mfDetY->RegisterPrimitive(primitiv2);
SetSensitiveDetector("logicScY", mfDetY);
}

58
src/EventAction.cc Normal file
View File

@@ -0,0 +1,58 @@
#include "EventAction.hh"
#include "Analysis.hh"
#include "G4Event.hh"
#include "G4HCofThisEvent.hh"
#include "G4RunManager.hh"
#include "G4SDManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4THitsMap.hh"
#include "G4UnitsTable.hh"
#include "RunAction.hh"
EventAction::EventAction(RunAction *)
: G4UserEventAction(), fCollID_detX(-1), fCollID_detY(-1) {}
EventAction::~EventAction() {}
void EventAction::BeginOfEventAction(const G4Event *) {}
void EventAction::EndOfEventAction(const G4Event *evt) {
G4HCofThisEvent *HCE = evt->GetHCofThisEvent();
if (!HCE) return;
if (fCollID_detX < 0) {
G4SDManager *SDMan = G4SDManager::GetSDMpointer();
fCollID_detX = SDMan->GetCollectionID("detX/edep");
}
if (fCollID_detY < 0) {
G4SDManager *SDMan = G4SDManager::GetSDMpointer();
fCollID_detY = SDMan->GetCollectionID("detY/edep");
}
auto anaMan = G4AnalysisManager::Instance();
G4THitsMap<G4double> *evtMap =
(G4THitsMap<G4double> *)(HCE->GetHC(fCollID_detX));
std::map<G4int, G4double *>::iterator itr;
for (itr = evtMap->GetMap()->begin(); itr != evtMap->GetMap()->end(); itr++) {
G4int copyNb = (itr->first);
G4double edep = *(itr->second);
// G4cout << "\n detX" << copyNb << ": " << edep / keV << " keV ";
anaMan->FillNtupleDColumn(copyNb, edep / MeV);
}
evtMap = (G4THitsMap<G4double> *)(HCE->GetHC(fCollID_detY));
for (itr = evtMap->GetMap()->begin(); itr != evtMap->GetMap()->end(); itr++) {
G4int copyNb = (itr->first);
G4double edep = *(itr->second);
// G4cout << "\n detX" << copyNb << ": " << edep / keV << " keV ";
anaMan->FillNtupleDColumn(copyNb + 40, edep / MeV);
}
anaMan->AddNtupleRow();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

77
src/PrimaryGeneratorAction.cc Normal file
View File

@@ -0,0 +1,77 @@
#include "PrimaryGeneratorAction.hh"
#include "G4Box.hh"
#include "G4LogicalVolume.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4RunManager.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorAction::PrimaryGeneratorAction()
: G4VUserPrimaryGeneratorAction(), fParticleGun(0), fEnvelopeBox(0) {
G4int n_particle = 1;
fParticleGun = new G4ParticleGun(n_particle);
// default particle kinematic
G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
G4String particleName;
G4ParticleDefinition *particle =
particleTable->FindParticle(particleName = "mu+");
fParticleGun->SetParticleDefinition(particle);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., -1.));
fParticleGun->SetParticleEnergy(6. * GeV);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fParticleGun; }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::GeneratePrimaries(G4Event *anEvent) {
// this function is called at the begining of ecah event
//
// In order to avoid dependence of PrimaryGeneratorAction
// on DetectorConstruction class we get Envelope volume
// from G4LogicalVolumeStore.
G4double worldSizeX = 0;
G4double worldSizeY = 0;
G4double worldSizeZ = 0;
if (!fEnvelopeBox) {
G4LogicalVolume *envLV =
G4LogicalVolumeStore::GetInstance()->GetVolume("World");
if (envLV) fEnvelopeBox = dynamic_cast<G4Box *>(envLV->GetSolid());
}
if (fEnvelopeBox) {
worldSizeX = fEnvelopeBox->GetXHalfLength() * 2.;
worldSizeY = fEnvelopeBox->GetYHalfLength() * 2.;
worldSizeZ = fEnvelopeBox->GetZHalfLength() * 2.;
} else {
G4ExceptionDescription msg;
msg << "Envelope volume of box shape not found.\n";
msg << "Perhaps you have changed geometry.\n";
msg << "The gun will be place at the center.";
G4Exception("PrimaryGeneratorAction::GeneratePrimaries()", "MyCode0002",
JustWarning, msg);
}
G4double x0 = worldSizeX * (G4UniformRand() - 0.5);
G4double y0 = worldSizeY * (G4UniformRand() - 0.5);
G4double z0 = 0.5 * worldSizeZ;
fParticleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
fParticleGun->GeneratePrimaryVertex(anEvent);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

60
src/RunAction.cc Normal file
View File

@@ -0,0 +1,60 @@
#include "RunAction.hh"
#include "Analysis.hh"
#include "DetectorConstruction.hh"
#include "PrimaryGeneratorAction.hh"
// #include "Run.hh"
#include "G4AccumulableManager.hh"
#include "G4LogicalVolume.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4Run.hh"
#include "G4RunManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4UnitsTable.hh"
RunAction::RunAction() : G4UserRunAction() {
auto anaMan = G4AnalysisManager::Instance();
G4cout << "Using " << anaMan->GetType() << G4endl;
// Create ntuple, first 40 columns are x hits, last 40 columns are y hits
anaMan->SetVerboseLevel(0);
anaMan->SetNtupleMerging(true);
anaMan->CreateNtuple("edep", "Energy deposit [MeV]");
char columnName[20];
for (int i = 0; i < 40; i++) {
snprintf(columnName, 20, "x%02d", i);
anaMan->CreateNtupleDColumn(G4String(columnName));
}
for (int i = 0; i < 40; i++) {
snprintf(columnName, 20, "y%02d", i);
anaMan->CreateNtupleDColumn(G4String(columnName));
}
anaMan->FinishNtuple();
}
RunAction::~RunAction() { delete G4AnalysisManager::Instance(); }
void RunAction::BeginOfRunAction(const G4Run *) {
// inform the runManager to save random number seed
G4RunManager::GetRunManager()->SetRandomNumberStore(false);
// Get analysis manager
auto anaMan = G4AnalysisManager::Instance();
// Open an output file
//
G4String fileName = "edep";
anaMan->OpenFile(fileName);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void RunAction::EndOfRunAction(const G4Run *) {
auto anaMan = G4AnalysisManager::Instance();
anaMan->Write();
anaMan->CloseFile();
}