Merge pull request #10 from carmenromo/nest-det

Make NEST use more general
petalo-project · Oct 5, 2022 · 27bfae4 · 27bfae4
2 parents 22e99f8 + 74c396a
commit 27bfae4
Show file tree

Hide file tree

Showing 5 changed files with 91 additions and 5 deletions.
diff --git a/macros/PET_full_body_nest.config.mac b/macros/PET_full_body_nest.config.mac
@@ -26,6 +26,7 @@
 /process/optical/processActivation Scintillation false
 /PhysicsList/Petalo/nest true
 #/PhysicsList/Petalo/thermal_electrons false
+/PhysicsList/Petalo/petalo_detector FullRing
 
 
 ### VERBOSITIES

diff --git a/macros/PetBox_nest.config.mac b/macros/PetBox_nest.config.mac
@@ -0,0 +1,42 @@
+### GEOMETRY
+/Geometry/PetBox/visibility true
+/Geometry/PetBox/tile_type_d HamamatsuVUV
+/Geometry/PetBox/tile_type_c HamamatsuVUV
+/Geometry/PetBox/single_tile_coinc_plane 0
+/Geometry/PetBox/add_teflon_block 1
+
+/Geometry/PetBox/tile_vis true
+/Geometry/PetBox/tile_refl 0.
+/Geometry/PetBox/sipm_time_binning 5. picosecond
+/Geometry/PetBox/sipm_pde 0.3
+
+### GENERATOR
+/Generator/IonGenerator/region SOURCE
+/Generator/IonGenerator/atomic_number 11
+/Generator/IonGenerator/mass_number 22
+#/Generator/SingleParticle/particle geantino
+#/Generator/SingleParticle/min_energy 1  eV
+#/Generator/SingleParticle/max_energy 1 keV
+#/Generator/SingleParticle/region CENTER
+
+/Actions/PetSensorsEventAction/min_charge 50
+
+### PHYSICS
+/process/optical/processActivation Scintillation false
+/process/optical/processActivation Cerenkov      true
+/PhysicsList/Petalo/nest true
+/PhysicsList/Petalo/thermal_electrons false
+/PhysicsList/Petalo/petalo_detector Petit
+
+### VERBOSITIES
+/run/verbose 0
+/event/verbose 0
+/tracking/verbose 0
+
+/process/em/verbose 0
+
+/nexus/random_seed 132323
+
+### OUTPUT FILE
+/nexus/persistency/start_id 0
+/nexus/persistency/outputFile petbox_nest.pet
diff --git a/macros/PetBox_nest.init.mac b/macros/PetBox_nest.init.mac
@@ -0,0 +1,24 @@
+### PHYSICS
+/PhysicsList/RegisterPhysics G4EmStandardPhysics_option4
+/PhysicsList/RegisterPhysics G4OpticalPhysics
+/PhysicsList/RegisterPhysics PetaloPhysics
+/PhysicsList/RegisterPhysics G4DecayPhysics
+/PhysicsList/RegisterPhysics G4RadioactiveDecayPhysics
+/PhysicsList/RegisterPhysics G4StepLimiterPhysics
+
+### GEOMETRY
+/nexus/RegisterGeometry PetBox
+
+### GENERATOR
+/nexus/RegisterGenerator IonGenerator
+#/nexus/RegisterGenerator SingleParticleGenerator
+
+### ACTIONS
+/nexus/RegisterRunAction DefaultRunAction
+/nexus/RegisterEventAction PetSensorsEventAction
+/nexus/RegisterTrackingAction PetaloTrackingAction
+/nexus/RegisterStackingAction PetNESTStackingAction
+
+/nexus/RegisterPersistencyManager PetaloPersistencyManager
+
+/nexus/RegisterMacro macros/PetBox_nest.config.mac
diff --git a/source/physics_lists/PetaloPhysics.cc b/source/physics_lists/PetaloPhysics.cc
@@ -12,7 +12,9 @@
 #include "PetaloPersistencyManager.h"
 
 #include <NESTProc.hh>
+#include <VDetector.hh>
 #include <PetaloDetector.hh>
+#include <PetitDetector.hh>
 
 #include <G4Scintillation.hh>
 #include <G4GenericMessenger.hh>
@@ -32,7 +34,8 @@ G4_DECLARE_PHYSCONSTR_FACTORY(PetaloPhysics);
 
 PetaloPhysics::PetaloPhysics() : G4VPhysicsConstructor("PetaloPhysics"),
                                  risetime_(false), noCompt_(false),
-                                 nest_(false), prod_th_el_(false)
+                                 nest_(false), prod_th_el_(false),
+                                 petalo_detector_("FullRing")
 {
   msg_ = new G4GenericMessenger(this, "/PhysicsList/Petalo/",
                                 "Control commands of the nexus physics list.");
@@ -48,6 +51,9 @@ PetaloPhysics::PetaloPhysics() : G4VPhysicsConstructor("PetaloPhysics"),
 
   msg_->DeclareProperty("thermal_electrons", prod_th_el_,
                         "If true, NEST thermal electrons are produced.");
+
+  msg_->DeclareProperty("petalo_detector", petalo_detector_,
+                        "Detector geometry chosen.");
 }
 
 PetaloPhysics::~PetaloPhysics()
@@ -108,14 +114,22 @@ void PetaloPhysics::ConstructProcess()
 
 
   if (nest_) {
-    PetaloDetector* petalo = new PetaloDetector();
-    G4double e_field = petalo->FitEF(0, 0, 0);
+    if (petalo_detector_ == "FullRing"){
+      petalo_ = new PetaloDetector();
+    } else if (petalo_detector_ == "Petit") {
+      petalo_ = new PetitDetector();
+    } else {
+      G4Exception("[PetaloPhysics]", "ConstructProcess()", FatalException,
+                  "Unknown PETALO detector!");
+    }
+
+    G4double e_field = petalo_->FitEF(0, 0, 0);
     PetaloPersistencyManager* pm =
       dynamic_cast<PetaloPersistencyManager*>(G4VPersistencyManager::GetPersistencyManager());
     pm->SetElectricField(e_field);
-    NEST::NESTcalc* petaloCalc = new NEST::NESTcalc(petalo);
+    NEST::NESTcalc* petaloCalc = new NEST::NESTcalc(petalo_);
     NEST::NESTProc* theNESTScintillationProcess =
-      new NEST::NESTProc("S1", fElectromagnetic, petaloCalc, petalo);
+      new NEST::NESTProc("S1", fElectromagnetic, petaloCalc, petalo_);
     theNESTScintillationProcess->SetDetailedSecondaries(true); // this is to use the full scintillation spectrum of LXe.
     theNESTScintillationProcess->SetStackElectrons(prod_th_el_); //false if only light is collected
 

diff --git a/source/physics_lists/PetaloPhysics.h b/source/physics_lists/PetaloPhysics.h
@@ -10,6 +10,7 @@
 #define PETALO_PHYSICS_H
 
 #include <G4VPhysicsConstructor.hh>
+#include <VDetector.hh>
 
 class G4GenericMessenger;
 class WavelengthShifting;
@@ -28,6 +29,8 @@ class PetaloPhysics : public G4VPhysicsConstructor
   /// Construct all required physics processes (Geant4 mandatory method)
   virtual void ConstructProcess();
 
+  VDetector *petalo_;
+
 private:
   G4bool risetime_; ///< Rise time for LYSO
 
@@ -37,6 +40,8 @@ class PetaloPhysics : public G4VPhysicsConstructor
 
   G4bool prod_th_el_; ///< If true, NEST thermal electrons are produced
 
+  G4String petalo_detector_;
+
   G4GenericMessenger *msg_;
   WavelengthShifting *wls_;
   PositronAnnihilation* pos_annihil_;