Merge pull request #228 from next-exp/lab-generator

Lab generator

macros/NEXT100_muons.config.mac

@@ -1,7 +1,7 @@
 ## ----------------------------------------------------------------------------
 ## nexus | NEXT100_muons.config.mac
 ##
-## Configuration macro to simulate muons according in the NEXT-100 geometry.
+## Configuration macro to simulate muons in the NEXT-100 geometry.
 ##
 ## The NEXT Collaboration
 ## ----------------------------------------------------------------------------

macros/NEXT100_muons.init.mac

@@ -1,5 +1,5 @@
 ## ----------------------------------------------------------------------------
-## nexus | NEXT100_muons_lsc.init.mac
+## nexus | NEXT100_muons.init.mac
 ##
 ## Initialization macro to simulate muons in the NEXT-100 geometry.
 ##

source/generators/MuonGenerator.cc

@@ -34,7 +34,8 @@ REGISTER_CLASS(MuonGenerator, G4VPrimaryGenerator)
 MuonGenerator::MuonGenerator():
   G4VPrimaryGenerator(), msg_(0), particle_definition_(0),
   use_lsc_dist_(true), axis_rotation_(150), rPhi_(NULL), user_dir_{}, energy_min_(0.),
-  energy_max_(0.), dist_name_("za"), bInitialize_(false), geom_(0), geom_solid_(0)
+  energy_max_(0.), dist_name_("za"), bInitialize_(false), geom_(0), geom_solid_(0),
+  gen_rad_(223.33*cm)
 {
   msg_ = new G4GenericMessenger(this, "/Generator/MuonGenerator/",
 				"Control commands of muongenerator.");
@@ -71,7 +72,15 @@ MuonGenerator::MuonGenerator():
   rotation.SetParameterName("azimuth", false);
   rotation.SetRange("azimuth>0.");
 
-  DetectorConstruction* detconst = (DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction();
+  // defaults to length of the NEXT100 shielding diagonal if not present.
+  G4GenericMessenger::Command& generation_radius =
+    msg_->DeclareProperty("gen_rad", gen_rad_, "Set radius for generation disc");
+  generation_radius.SetUnitCategory("Length");
+  generation_radius.SetParameterName("gen_rad", false);
+  generation_radius.SetRange("gen_rad>0.");
+
+  DetectorConstruction* detconst =
+    (DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction();
   geom_ = detconst->GetGeometry();
 
 }
@@ -201,7 +210,7 @@ void MuonGenerator::GeneratePrimaryVertex(G4Event* event)
   // Particle properties
   G4double mass          = particle_definition_->GetPDGMass();
   G4double energy        = kinetic_energy + mass;
-  G4ThreeVector position = geom_->GenerateVertex(region_);
+  //G4ThreeVector position = geom_->GenerateVertex(region_);
 
   // Set default momentum and angular variables
   G4ThreeVector p_dir;
@@ -211,7 +220,7 @@ void MuonGenerator::GeneratePrimaryVertex(G4Event* event)
   // Momentum, zenith, azimuth (and energy) from angular distribution file
   if (use_lsc_dist_){
     GetDirection(p_dir, zenith, azimuth, energy, kinetic_energy, mass);
-    position = geom_->GenerateVertex(region_);
+    //    position = geom_->GenerateVertex(region_);
   }
   else {
 
@@ -236,9 +245,10 @@ void MuonGenerator::GeneratePrimaryVertex(G4Event* event)
       p_dir *= *rPhi_;
 
     }
-
   }
 
+  G4ThreeVector position = ProjectToVertex(p_dir);
+
   G4double pmod   = std::sqrt(energy*energy - mass*mass);
   G4double px = pmod * p_dir.x();
   G4double py = pmod * p_dir.y();
@@ -327,6 +337,36 @@ void MuonGenerator::GetDirection(G4ThreeVector& dir, G4double& zenith, G4double&
 }
 
 
+G4ThreeVector MuonGenerator::ProjectToVertex(const G4ThreeVector& dir)
+{
+  /////////////////////////////////////////////////////////////////////////
+  // This method of vertex generation decides the starting vertex
+  // of the primary particle in three steps:
+  // 1) A random point is generated on a disc centred on the main
+  //    volumes.
+  // 2) This disc is rotated so that it is perpendicular to the
+  //    direction vector which is the only argument of the function.
+  //    The new point is a point through which the vector passes.
+  // 3) The vertex is found by projecting backwards along the direction
+  //    vector from that point to find the point where the ray
+  //    (point - t*dir) intersects with the region configured as the
+  //    starting point for all vertices.
+  /////////////////////////////////////////////////////////////////////////
+  // Postion in disc
+  G4double radius = gen_rad_ * std::sqrt(G4UniformRand());
+  G4double ang = 2 * G4UniformRand() * pi;
+
+  // Continue assuming that Y is vertical and z drift,
+  // valid for NEW and NEXT-100 (at least).
+  // Rotate the disc (origin-point vector) to be perpendicular to dir.
+  G4ThreeVector point(radius * std::cos(ang), 0., radius * std::sin(ang));
+  point.rotate(pi / 2 - dir.angle(point), dir.cross(point));
+
+  // Now project back to the requested region intersection.
+  return geom_->ProjectToRegion(region_, point, -dir);
+}
+
+
 G4double MuonGenerator::GetZenith() const
 {
   return fRandomGeneral_->fire()*pi/2;

source/generators/MuonGenerator.h

@@ -51,6 +51,8 @@ namespace nexus {
     void GetDirection(G4ThreeVector& dir, G4double& zenith, G4double& azimuth,
                       G4double& energy, G4double& kinetic_energy, G4double mass);
 
+    G4ThreeVector ProjectToVertex(const G4ThreeVector& dir);
+
     G4bool CheckOverlap(const G4ThreeVector& vtx, const G4ThreeVector& dir);
 
     /// Load in the Muon Angular/Energy Distribution from CSV file
@@ -95,6 +97,8 @@ namespace nexus {
     std::vector<G4double> energy_smear_;  ///< List of Energy bin smear values
     G4RandGeneral *fRandomGeneral_; ///< Pointer to the RNG flux distribution
 
+    G4double gen_rad_; ///< Radius of disc for generation
+
   };
 
 } // end namespace nexus

source/geometries/GeometryBase.h

@@ -34,6 +34,12 @@ namespace nexus {
     /// Returns a point within a given region of the geometry
     virtual G4ThreeVector GenerateVertex(const G4String&) const;
 
+    /// Returns a point within a region projecting from a
+    /// given point backwards along a line.
+    virtual G4ThreeVector ProjectToRegion(const G4String&,
+					  const G4ThreeVector&,
+					  const G4ThreeVector&) const;
+
     /// Returns the span (maximum dimension) of the geometry
     G4double GetSpan();
 
@@ -102,6 +108,11 @@ namespace nexus {
   inline G4ThreeVector GeometryBase::GenerateVertex(const G4String&) const
   { return G4ThreeVector(0., 0., 0.); }
 
+  inline G4ThreeVector GeometryBase::ProjectToRegion(const G4String&,
+						     const G4ThreeVector&,
+						     const G4ThreeVector&) const
+  { return G4ThreeVector(0., 0., 0.); }
+
   inline void GeometryBase::SetSpan(G4double s) { span_ = s; }
 
   inline G4double GeometryBase::GetSpan() { return span_; }

source/geometries/LSCHallA.cc

@@ -119,4 +119,23 @@ namespace nexus {
     return vertex;
   }
 
+  G4ThreeVector LSCHallA::ProjectToRegion(const G4String& region,
+					  const G4ThreeVector& point,
+					  const G4ThreeVector& dir) const
+  {
+    // Project along dir from point to find the first intersection
+    // with region.
+    G4ThreeVector vertex(0., 0., 0.);
+    if (region == "HALLA_INNER")
+      return hallA_vertex_gen_->GetIntersect(point, dir);
+    else if (region == "HALLA_OUTER")
+      return hallA_outer_gen_->GetIntersect(point, dir);
+    else {
+      G4Exception("[LSCHallA]", "ProjectToRegion()", FatalException,
+		  "Unknown vertex generation region!");
+    }
+
+    return vertex;
+  }
+
 }

source/geometries/LSCHallA.h

@@ -30,6 +30,12 @@ namespace nexus {
     /// Generate a vertex within a given region of the geometry
     G4ThreeVector GenerateVertex(const G4String& region) const;
 
+    /// Returns a point within a region projecting from a
+    /// given point backwards along a line.
+    G4ThreeVector ProjectToRegion(const G4String& region,
+				  const G4ThreeVector& point,
+				  const G4ThreeVector& dir) const;
+
     /// Builder
     void Construct();
 

source/geometries/Next100.cc

@@ -271,4 +271,29 @@ namespace nexus {
     return vertex;
   }
 
+
+  G4ThreeVector Next100::ProjectToRegion(const G4String& region,
+					 const G4ThreeVector& point,
+					 const G4ThreeVector& dir) const
+  {
+    // Project along dir from point to find the first intersection
+    // with region.
+    G4ThreeVector vertex(0., 0., 0.);
+    if (region == "EXTERNAL"){
+      return shielding_->ProjectToRegion(region, point, dir);
+    }
+    else if ((region == "HALLA_OUTER") || (region == "HALLA_INNER")){
+      if (!lab_walls_)
+	G4Exception("[Next100]", "ProjectToRegion()", FatalException,
+                    "To project to this region you need lab_walls == true!");
+      return hallA_walls_->ProjectToRegion(region, point, dir);
+    }
+    else {
+      G4Exception("[Next100]", "ProjectToRegion()", FatalException,
+		  "Unknown vertex generation region!");
+    }
+
+    return vertex + G4ThreeVector(0., 0., -gate_zpos_in_vessel_);
+  }
+
 } //end namespace nexus

source/geometries/Next100.h

@@ -38,6 +38,12 @@ namespace nexus {
     /// Generate a vertex within a given region of the geometry
     G4ThreeVector GenerateVertex(const G4String& region) const;
 
+    /// Returns a point within a region projecting from a
+    /// given point backwards along a line.
+    G4ThreeVector ProjectToRegion(const G4String& region,
+				  const G4ThreeVector& point,
+				  const G4ThreeVector& dir) const;
+
 
   private:
     void BuildLab();