From 6a5eab0052fd3cc963a0f21e8193418adac2905a Mon Sep 17 00:00:00 2001 From: Nam Tran Date: Thu, 14 Jan 2021 16:29:21 -0600 Subject: [PATCH] initial working version: 80 scintillators in 8 layers, root output, mono-energetic muon on -Z direction --- .gitignore | 2 + CMakeLists.txt | 68 +++++++++++++ include/ActionInitialization.hh | 23 +++++ include/Analysis.hh | 6 ++ include/DetectorConstruction.hh | 30 ++++++ include/EventAction.hh | 32 ++++++ include/PrimaryGeneratorAction.hh | 37 +++++++ include/RunAction.hh | 31 ++++++ macros/batch_10k.mac | 3 + macros/batch_1M.mac | 3 + macros/init_vis.mac | 13 +++ macros/run1.mac | 23 +++++ macros/run2.mac | 22 +++++ macros/steal.sh | 17 ++++ macros/vis.mac | 101 +++++++++++++++++++ muography.cc | 117 ++++++++++++++++++++++ src/ActionInitialization.cc | 35 +++++++ src/DetectorConstruction.cc | 157 ++++++++++++++++++++++++++++++ src/EventAction.cc | 58 +++++++++++ src/PrimaryGeneratorAction.cc | 77 +++++++++++++++ src/RunAction.cc | 60 ++++++++++++ 21 files changed, 915 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 include/ActionInitialization.hh create mode 100644 include/Analysis.hh create mode 100644 include/DetectorConstruction.hh create mode 100644 include/EventAction.hh create mode 100644 include/PrimaryGeneratorAction.hh create mode 100644 include/RunAction.hh create mode 100644 macros/batch_10k.mac create mode 100644 macros/batch_1M.mac create mode 100644 macros/init_vis.mac create mode 100644 macros/run1.mac create mode 100644 macros/run2.mac create mode 100755 macros/steal.sh create mode 100644 macros/vis.mac create mode 100644 muography.cc create mode 100644 src/ActionInitialization.cc create mode 100644 src/DetectorConstruction.cc create mode 100644 src/EventAction.cc create mode 100644 src/PrimaryGeneratorAction.cc create mode 100644 src/RunAction.cc diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..9ef9604 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +build + diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..63bf051 --- /dev/null +++ b/CMakeLists.txt @@ -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) diff --git a/include/ActionInitialization.hh b/include/ActionInitialization.hh new file mode 100644 index 0000000..e17e0f7 --- /dev/null +++ b/include/ActionInitialization.hh @@ -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 + + diff --git a/include/Analysis.hh b/include/Analysis.hh new file mode 100644 index 0000000..8b47338 --- /dev/null +++ b/include/Analysis.hh @@ -0,0 +1,6 @@ +#ifndef Analysis_h +#define Analysis_h + +#include "g4root.hh" + +#endif /* end of include guard */ diff --git a/include/DetectorConstruction.hh b/include/DetectorConstruction.hh new file mode 100644 index 0000000..8195759 --- /dev/null +++ b/include/DetectorConstruction.hh @@ -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 + diff --git a/include/EventAction.hh b/include/EventAction.hh new file mode 100644 index 0000000..7d32c9a --- /dev/null +++ b/include/EventAction.hh @@ -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 + + diff --git a/include/PrimaryGeneratorAction.hh b/include/PrimaryGeneratorAction.hh new file mode 100644 index 0000000..fc18e36 --- /dev/null +++ b/include/PrimaryGeneratorAction.hh @@ -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 diff --git a/include/RunAction.hh b/include/RunAction.hh new file mode 100644 index 0000000..052a901 --- /dev/null +++ b/include/RunAction.hh @@ -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 + diff --git a/macros/batch_10k.mac b/macros/batch_10k.mac new file mode 100644 index 0000000..b4659bf --- /dev/null +++ b/macros/batch_10k.mac @@ -0,0 +1,3 @@ +/run/initialize +/tracking/verbose 0 +/run/beamOn 10000 diff --git a/macros/batch_1M.mac b/macros/batch_1M.mac new file mode 100644 index 0000000..6e41788 --- /dev/null +++ b/macros/batch_1M.mac @@ -0,0 +1,3 @@ +/run/initialize +/tracking/verbose 0 +/run/beamOn 1000000 diff --git a/macros/init_vis.mac b/macros/init_vis.mac new file mode 100644 index 0000000..7f15a59 --- /dev/null +++ b/macros/init_vis.mac @@ -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 diff --git a/macros/run1.mac b/macros/run1.mac new file mode 100644 index 0000000..028208d --- /dev/null +++ b/macros/run1.mac @@ -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 diff --git a/macros/run2.mac b/macros/run2.mac new file mode 100644 index 0000000..1d00add --- /dev/null +++ b/macros/run2.mac @@ -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 diff --git a/macros/steal.sh b/macros/steal.sh new file mode 100755 index 0000000..ff40d95 --- /dev/null +++ b/macros/steal.sh @@ -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 diff --git a/macros/vis.mac b/macros/vis.mac new file mode 100644 index 0000000..4625ee0 --- /dev/null +++ b/macros/vis.mac @@ -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 diff --git a/muography.cc b/muography.cc new file mode 100644 index 0000000..99aecb5 --- /dev/null +++ b/muography.cc @@ -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..... diff --git a/src/ActionInitialization.cc b/src/ActionInitialization.cc new file mode 100644 index 0000000..ced8f31 --- /dev/null +++ b/src/ActionInitialization.cc @@ -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...... diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc new file mode 100644 index 0000000..970b0b3 --- /dev/null +++ b/src/DetectorConstruction.cc @@ -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); +} diff --git a/src/EventAction.cc b/src/EventAction.cc new file mode 100644 index 0000000..071476a --- /dev/null +++ b/src/EventAction.cc @@ -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 *evtMap = + (G4THitsMap *)(HCE->GetHC(fCollID_detX)); + std::map::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 *)(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...... diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc new file mode 100644 index 0000000..b0b475e --- /dev/null +++ b/src/PrimaryGeneratorAction.cc @@ -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(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...... diff --git a/src/RunAction.cc b/src/RunAction.cc new file mode 100644 index 0000000..15c9e5e --- /dev/null +++ b/src/RunAction.cc @@ -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(); +}