confinement: add option to weigh by mass instead of volume

docs/manual/confinement.md

@@ -168,6 +168,19 @@ All sampling modes described above are available, with few notes/limitations:
   with a small but usually acceptable performance cost, however very large
   values can significantly slow down the simulation.
 
+## Weighting modes
+
+For the sampling in volumes, _remage_ by default weights the probability to find
+a sampled vertex in one of the volumes by the ratio of their volume to the total
+sampling volume. Other weighting modes can be selected by the user with a macro
+command:
+
+- by the mass of the volumes can with the command
+  <project:../rmg-commands.md#rmggeneratorconfinementsampleweightbymass>.
+- by the mass of a given isotope (specified by the pair of proton and neutron
+  number) in the materials with the command
+  <project:../rmg-commands.md#rmggeneratorconfinementsampleweightbymassisotope>.
+
 ## Vertices from external files
 
 For more complicated vertex confinement _remage_ supports the possibility to
@@ -202,6 +215,14 @@ capabilities described above:
 /RMG/Generator/Confinement/Physical/AddVolume [BCPV]\w+
 ```
 
+This can now be used as an example to weigh the volumes not by their volume: It
+could be useful to weigh the different volumes by the mass fraction of the
+isotope $^{76}$Ge, i.e. for simulating double beta decays with bxdecay0:
+
+```geant4
+/RMG/Generator/Confinement/SampleWeightByMassIsotope 32 76
+```
+
 ### Active LAr region of LEGEND-200
 
 This is a typical example of the sampling from a volume intersection. This will

include/RMGPrimaryTransformer.hh

@@ -21,7 +21,7 @@
 class RMGPrimaryTransformer : public G4PrimaryTransformer {
   public:
 
-    RMGPrimaryTransformer() : G4PrimaryTransformer() {
+    RMGPrimaryTransformer() {
       // this suppresses the ZeroPolarization warning emitted for optical photons.
       nWarn = 11;
     };

include/RMGVertexConfinement.hh

@@ -113,6 +113,12 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
 
     void SetSamplingMode(SamplingMode mode) { fSamplingMode = mode; }
     void SetFirstSamplingVolumeType(VolumeType type) { fFirstSamplingVolumeType = type; }
+    void SetWeightByMass(bool mode) { fWeightByMass = mode; }
+    void SetWeightByMassIsotope(int z, int n) {
+      fWeightByMass = true;
+      fWeightByMassIsotopeZ = z;
+      fWeightByMassIsotopeN = n;
+    }
 
     std::vector<GenericGeometricalSolidData>& GetGeometricalSolidDataList() {
       return fGeomVolumeData;
@@ -238,12 +244,15 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
          */
         void GetDirection(G4ThreeVector& dir, G4ThreeVector& pos) const;
 
+        void RecalcMass(int z, int n);
+
         G4VPhysicalVolume* physical_volume = nullptr;
         G4VSolid* sampling_solid = nullptr;
         G4RotationMatrix rotation;
         G4ThreeVector translation;
 
         double volume = -1;
+        double mass = -1;
         double surface = -1;
 
         bool surface_sample = false;
@@ -266,8 +275,9 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
 
         /** @brief Select a @c SampleableObject from the collection, weighted by volume.
          *  @returns a reference to the chosen @c SampleableObject .
+         *  @param weight_by_mass A flag of whether the volume weighting should be done by mass and not by volume.
          */
-        [[nodiscard]] const SampleableObject& VolumeWeightedRand() const;
+        [[nodiscard]] const SampleableObject& VolumeWeightedRand(bool weight_by_mass) const;
         [[nodiscard]] bool IsInside(const G4ThreeVector& vertex) const;
 
         // emulate @c std::vector
@@ -279,10 +289,16 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
         void clear() { data.clear(); }
         void insert(SampleableObjectCollection& other) {
           for (size_t i = 0; i < other.size(); ++i) this->emplace_back(other.at(i));
+          this->total_volume += other.total_volume;
+          this->total_mass += other.total_mass;
+          this->total_surface += other.total_surface;
         }
 
+        void recalc_total(bool weigh_by_mass, int mass_isotope_z, int mass_istotope_n);
+
         std::vector<SampleableObject> data;
         double total_volume = 0;
+        double total_mass = 0;
         double total_surface = 0;
     };
 
@@ -311,6 +327,9 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
 
     SamplingMode fSamplingMode = SamplingMode::kUnionAll;
     VolumeType fFirstSamplingVolumeType = VolumeType::kUnset;
+    bool fWeightByMass = false;
+    int fWeightByMassIsotopeZ = 0;
+    int fWeightByMassIsotopeN = 0;
 
     bool fOnSurface = false;
     bool fForceContainmentCheck = false;