target size can also be customized

Readme.md

@@ -0,0 +1 @@
+Muon tomography simulation

include/DetectorConstruction.hh

@@ -7,6 +7,7 @@
 
 class G4VPhysicalVolume;
 class G4LogicalVolume;
+class DetectorConstructionMessenger;
 
 /// Detector construction class to define materials and geometry.
 
@@ -19,9 +20,44 @@ class DetectorConstruction : public G4VUserDetectorConstruction
     virtual G4VPhysicalVolume* Construct();
     virtual void ConstructSDandField();
 
+    void ComputeGeometry();
+
+    // public fields to avoid all getters and setters
+    G4double scSizeX;
+    G4double scSizeY;
+    G4double scSizeZ;
+    G4double scGapX;
+    G4double scGapZ;
+    G4int nScBar;
+    G4int nLayerTop;
+    G4int nLayerBottom;
+
+    G4double layerGapZ;
+    G4double topBotGap;
+
+    G4double layerSizeX;
+    G4double layerSizeY;
+    G4double layerSizeZ;
+
+    G4double concreteX;
+    G4double concreteY;
+    G4double concreteZ;
+    G4double targetX;
+    G4double targetY;
+    G4double targetZ;
+    G4double targetLocationX;
+    G4double targetLocationY;
+    G4double targetLocationZ;
+    G4String targetMaterialStr;
+
+    G4double worldSizeX;
+    G4double worldSizeY;
+    G4double worldSizeZ;
+
   protected:
     // void DefineMaterials();
-    G4bool  fCheckOverlaps;
+    G4bool fCheckOverlaps;
+    DetectorConstructionMessenger *fDetConMessenger;
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

include/DetectorConstructionMessenger.hh

@@ -0,0 +1,36 @@
+#ifndef DetectorConstructionMessenger_h
+#define DetectorConstructionMessenger_h 1
+
+class DetectorConstruction;
+class G4UIdirectory;
+class G4UIcmdWithADoubleAndUnit;
+class G4UIcmdWithAIntAndUnit;
+class G4UIcmdWith3VectorAndUnit;
+class G4UIcmdWith3Vector;
+class G4UIcmdWithABool;
+class G4UIcmdWithAString;
+class G4UIcmdWithADouble;
+
+#include "G4UImessenger.hh"
+#include "globals.hh"
+
+class DetectorConstructionMessenger: public G4UImessenger
+{
+  public:
+    DetectorConstructionMessenger(DetectorConstruction *detCon);
+    ~DetectorConstructionMessenger();
+
+    virtual void SetNewValue(G4UIcommand *command, G4String newValues);
+    virtual G4String GetCurrentValue(G4UIcommand * command);
+
+  private:
+    DetectorConstruction *fDetCon;
+
+    G4UIdirectory *fMuTomoDirectory;
+
+    G4UIcmdWithAString *fTargetMaterialCmd;
+    G4UIcmdWith3VectorAndUnit *fTargetSizeCmd;
+    G4UIcmdWith3VectorAndUnit *fTargetLocationCmd;
+};
+
+#endif

macros/batch_1M.mac

@@ -1,3 +1,8 @@
+# customize target parameters
+/muTomo/targetLocation 0 10 0 cm
+/muTomo/targetSize 10 10 20 cm
+/muTomo/targetMaterial G4_W
+
 /run/initialize
 /tracking/verbose 0
 /run/beamOn 1000000

macros/run1.mac

@@ -1,23 +0,0 @@
-#
-# 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

macros/run2.mac

@@ -1,22 +0,0 @@
-#/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

src/DetectorConstruction.cc

@@ -1,5 +1,6 @@
 #include "DetectorConstruction.hh"
 
+#include "DetectorConstructionMessenger.hh"
 #include "G4Box.hh"
 #include "G4LogicalVolume.hh"
 #include "G4MultiFunctionalDetector.hh"
@@ -20,49 +21,27 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 DetectorConstruction::DetectorConstruction()
-    : G4VUserDetectorConstruction(), fCheckOverlaps(true) {
+    : G4VUserDetectorConstruction(), targetMaterialStr("G4_Pb"),
+      fCheckOverlaps(false) {
+  targetX = 10 * cm;
+  targetY = 10 * cm;
+  targetZ = 5 * cm;
+  fDetConMessenger = new DetectorConstructionMessenger(this);
   // DefineMaterials();
 }
 
-DetectorConstruction::~DetectorConstruction() {}
+DetectorConstruction::~DetectorConstruction() { delete fDetConMessenger; }
 
 G4VPhysicalVolume *DetectorConstruction::Construct() {
   G4NistManager *nist = G4NistManager::Instance();
   G4Material *default_mat = nist->FindOrBuildMaterial("G4_AIR");
-  G4Material *targetMaterial = nist->FindOrBuildMaterial("G4_Pb");
+  G4Material *targetMaterial = nist->FindOrBuildMaterial(targetMaterialStr);
   G4Material *concreteMaterial = nist->FindOrBuildMaterial("G4_CONCRETE");
   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 concreteX = 30 * cm;
-  G4double concreteY = 30 * cm;
-  G4double concreteZ = 30 * cm;
-  G4double targetX = 10 * cm;
-  G4double targetY = 10 * cm;
-  G4double targetZ = 5 * cm;
-
-  G4double worldSizeX = 1.2 * layerSizeX;
-  G4double worldSizeY = 1.2 * layerSizeY;
-  G4double worldSizeZ =
-      1.2 * ((nLayerBottom + nLayerTop) * (layerSizeZ + layerGapZ) + topBotGap);
-
+  // ComputeGeometry first
+  ComputeGeometry();
   // world volume
   G4Box *solidWorld =
       new G4Box("World", 0.5 * worldSizeX, 0.5 * worldSizeY, 0.5 * worldSizeZ);
@@ -90,8 +69,9 @@ G4VPhysicalVolume *DetectorConstruction::Construct() {
       new G4Box("target", targetX / 2, targetY / 2, targetZ / 2);
   G4LogicalVolume *logTarget =
       new G4LogicalVolume(solidTarget, targetMaterial, "target");
-  new G4PVPlacement(0, G4ThreeVector(), logTarget, "target", logConcrete, false,
-                    0, fCheckOverlaps);
+  new G4PVPlacement(
+      0, G4ThreeVector(targetLocationX, targetLocationY, targetLocationZ),
+      logTarget, "target", logConcrete, false, 0, fCheckOverlaps);
 
   // scintillators
   G4Box *solidScX =
@@ -177,3 +157,30 @@ void DetectorConstruction::ConstructSDandField() {
   mfDetY->RegisterPrimitive(primitiv2);
   SetSensitiveDetector("logicScY", mfDetY);
 }
+
+void DetectorConstruction::ComputeGeometry() {
+  scSizeX = 5.0 * cm;
+  scSizeY = 50.0 * cm;
+  scSizeZ = 1. * cm;
+  scGapX = 0.2 * cm;
+  scGapZ = 0.2 * cm;
+  nScBar = 10;
+  nLayerTop = 2;
+  nLayerBottom = 2;
+
+  layerGapZ = 1 * cm;
+  topBotGap = 50 * cm;
+
+  layerSizeX = nScBar * scSizeX + nScBar * scGapX;
+  layerSizeY = nScBar * scSizeX + nScBar * scGapX;
+  layerSizeZ = scSizeZ + scGapZ;
+
+  concreteX = 30 * cm;
+  concreteY = 30 * cm;
+  concreteZ = 30 * cm;
+
+  worldSizeX = 1.2 * layerSizeX;
+  worldSizeY = 1.2 * layerSizeY;
+  worldSizeZ =
+      1.2 * ((nLayerBottom + nLayerTop) * (layerSizeZ + layerGapZ) + topBotGap);
+}

src/DetectorConstructionMessenger.cc

@@ -0,0 +1,91 @@
+#include "DetectorConstructionMessenger.hh"
+#include "DetectorConstruction.hh"
+#include "G4UIdirectory.hh"
+#include "G4UIcmdWithADoubleAndUnit.hh"
+#include "G4UIcmdWithABool.hh"
+#include "G4UIcmdWithAString.hh"
+#include "G4UIcmdWith3VectorAndUnit.hh"
+#include "G4UIcmdWith3Vector.hh"
+#include "G4UIcmdWithADouble.hh"
+#include "G4RunManager.hh"
+
+#include "G4ios.hh"
+
+DetectorConstructionMessenger::DetectorConstructionMessenger(
+    DetectorConstruction *detCon)
+    : fDetCon(detCon) {
+  fMuTomoDirectory = new G4UIdirectory("/muTomo/");
+  fMuTomoDirectory->SetGuidance("Customize muon tomography geometry commands.");
+
+  fTargetMaterialCmd = new G4UIcmdWithAString("/muTomo/targetMaterial", this);
+  fTargetMaterialCmd->SetGuidance("Set target material.");
+  fTargetMaterialCmd->SetParameterName("targetMaterial", true);
+  fTargetMaterialCmd->SetDefaultValue("G4_Pb");
+
+  fTargetSizeCmd = new G4UIcmdWith3VectorAndUnit("/muTomo/targetSize", this);
+  fTargetSizeCmd->SetGuidance("Set target size.");
+  fTargetSizeCmd->SetParameterName("targetSizeX", "targetSizeY", "targetSizeZ",
+                                   true);
+  fTargetSizeCmd->SetDefaultValue(G4ThreeVector(6., 2., 6.));
+  fTargetSizeCmd->SetDefaultUnit("mm");
+
+  fTargetLocationCmd =
+      new G4UIcmdWith3VectorAndUnit("/muTomo/targetLocation", this);
+  fTargetLocationCmd->SetGuidance("Set target location.");
+  fTargetLocationCmd->SetParameterName("targetLocationX", "targetLocationY",
+                                       "targetLocationZ", true);
+  fTargetLocationCmd->SetDefaultValue(G4ThreeVector(0., 0., 0.));
+  fTargetLocationCmd->SetDefaultUnit("mm");
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+DetectorConstructionMessenger::~DetectorConstructionMessenger() {
+  delete fTargetMaterialCmd;
+  delete fTargetSizeCmd;
+  delete fTargetLocationCmd;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void DetectorConstructionMessenger::SetNewValue(G4UIcommand *command,
+                                                G4String newValue) {
+  if (command == fTargetMaterialCmd) {
+    fDetCon->targetMaterialStr = newValue;
+  }
+  if (command == fTargetSizeCmd) {
+    G4ThreeVector newPos = fTargetSizeCmd->GetNew3VectorValue(newValue);
+    fDetCon->targetX = newPos.x();
+    fDetCon->targetY = newPos.y();
+    fDetCon->targetZ = newPos.z();
+  }
+  if (command == fTargetLocationCmd) {
+    G4ThreeVector newLoc = fTargetLocationCmd->GetNew3VectorValue(newValue);
+    fDetCon->targetLocationX = newLoc.x();
+    fDetCon->targetLocationY = newLoc.y();
+    fDetCon->targetLocationZ = newLoc.z();
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4String DetectorConstructionMessenger::GetCurrentValue(G4UIcommand *command) {
+  G4String cv;
+
+  if (command == fTargetMaterialCmd) {
+    cv = fDetCon->targetMaterialStr;
+  }
+  if (command == fTargetSizeCmd) {
+    cv = fTargetSizeCmd->ConvertToString(
+        G4ThreeVector(fDetCon->targetX, fDetCon->targetY, fDetCon->targetZ));
+  }
+  if (command == fTargetLocationCmd) {
+    cv = fTargetLocationCmd->ConvertToString(
+        G4ThreeVector(fDetCon->targetLocationX, fDetCon->targetLocationY,
+                      fDetCon->targetLocationZ));
+  }
+
+  return cv;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....