Allow Rp Distribution in ring Bunch Distribution.

examples/features/beam/CMakeLists.txt

@@ -60,6 +60,8 @@ simple_testing(bunch-reference-backwards      "--file=reference-backwards.gmad"
 simple_testing(bunch-reference-p0             "--file=reference-p0.gmad"            "")
 simple_testing(bunch-reference-ek0            "--file=reference-ek0.gmad"           "")
 simple_testing(bunch-ring                     "--file=ring.gmad"                    "")
+simple_testing(bunch-ring-rp                  "--file=ring-rp.gmad"                 "")
+simple_testing(bunch-ring-rp-range            "--file=ring-rp-range.gmad"           "")
 simple_testing(bunch-sphere                   "--file=sphere.gmad"                  "")
 simple_testing(bunch-sixtrack                 "--file=sixtrack.gmad"                "")
 simple_testing(bunch-slow-ext                 "--file=slow-ext.gmad"                "")

examples/features/beam/bad-params/CMakeLists.txt

@@ -24,7 +24,8 @@ simple_fail(bunch-bad-no-energy               "--file=no-energy.gmad"
 simple_fail(bunch-bad-no-bunch-frequency      "--file=no-bunch-frequency.gmad"       "")
 simple_fail(bunch-bad-sigmaE-and-sigmaP       "--file=sigmaE-and-sigmaP.gmad"        "")
 simple_fail(bunch-ptc-bad-parms               "--file=ptc-bad.gmad"                  "")
-simple_fail(bunch-ring-bad-parms              "--file=ring-bad.gmad"                 "")
+simple_fail(bunch-ring-bad-params1            "--file=ring-bad.gmad"                 "")
+simple_fail(bunch-ring-bad-params2            "--file=ring-rp-bad1.gmad"             "")
 simple_fail(bunch-sixtrack-bad-parms          "--file=sixtrack-bad.gmad"             "")
 simple_fail(bunch-slowext-bad                 "--file=slowext-bad.gmad"              "")
 simple_fail(bunch-square-bad-parms            "--file=square-bad.gmad"               "")

examples/features/beam/bad-params/ring-rp-bad1.gmad

@@ -0,0 +1,16 @@
+beam,  particle="e-",
+       energy=1.0 * GeV,
+       distrType="ring",
+       Rmin  = 5.0*mm, 
+       Rmax  = 5.1*mm,
+       Rpmin = 60*mrad,
+       Rpmax = 50*mrad,
+       sigmaT  = 1e-9,
+       sigmaE  = 1e-5;
+
+! rpmin > rpmax
+
+include options.gmad
+include fodo.gmad;
+
+option, ngenerate=10;

examples/features/beam/ring-rp-range.gmad

@@ -0,0 +1,14 @@
+beam,  particle="e-",
+       energy=1.0 * GeV,
+       distrType="ring",
+       Rmin  = 5.0*mm, 
+       Rmax  = 6.0*mm,
+       Rpmin = 10*mrad,
+       Rpmax = 20*mrad,
+       sigmaT  = 1e-9,
+       sigmaE  = 1e-5;
+
+include options.gmad
+include fodo.gmad;
+
+option, ngenerate=10;

examples/features/beam/ring-rp.gmad

@@ -0,0 +1,14 @@
+beam,  particle="e-",
+       energy=1.0 * GeV,
+       distrType="ring",
+       Rmin  = 5.0*mm, 
+       Rmax  = 5.1*mm,
+       Rpmin = 50*mrad,
+       Rpmax = 50*mrad,
+       sigmaT  = 1e-9,
+       sigmaE  = 1e-5;
+
+include options.gmad
+include fodo.gmad;
+
+option, ngenerate=10;

include/BDSBunchRing.hh

@@ -42,16 +42,20 @@ public:
   BDSBunchRing(BDSBunchRing&) = delete;
   /// @}
   virtual void SetOptions(const BDSParticleDefinition* beamParticle,
-			  const GMAD::Beam& beam,
-			  const BDSBunchType& distrType,
-			  G4Transform3D beamlineTransformIn = G4Transform3D::Identity,
-			  const G4double beamlineS = 0);
+                          const GMAD::Beam& beam,
+                          const BDSBunchType& distrType,
+                          G4Transform3D beamlineTransformIn = G4Transform3D::Identity,
+                          const G4double beamlineS = 0);
   virtual void CheckParameters();
   virtual BDSParticleCoordsFull GetNextParticleLocal();
 
 protected:
   G4double rMin;
   G4double rMax;
+  G4double rpMin;
+  G4double rpMax;
+  G4bool nonZeroRp;
+  G4bool generateRp;
 };
 
 #endif

include/BDSOutputROOTEventBeam.hh

@@ -40,7 +40,7 @@ public:
 #ifndef __ROOTBUILD__   
   void Fill();
 #endif
-  ClassDef(BDSOutputROOTEventBeam,6);
+  ClassDef(BDSOutputROOTEventBeam,7);
 };
 
 #endif

manual/source/model_control.rst

@@ -773,22 +773,32 @@ ring
 
 The ring distribution randomly and uniformly distributes particles around a circle in `x` and `y`. Then,
 for a given x,y the radius is randomly and uniformly in density distributed in that annulus. For
-all other parameters, the `reference`_ coordinates are used, i.e. `xp`, `yp` etc.
+all other parameters, the `reference`_ coordinates are used. Optionally, a similar distribution in
+`xp` and `yp` can be provided between two radii. If not specified, then no angle is given.
 
 
 * All parameters from `reference`_ distribution are used as centroids.
 
 .. tabularcolumns:: |p{5cm}|p{10cm}|
 
-+----------------------------------+-------------------------------------------------------+
-| **Variable**                     | **Description**                                       |
-+==================================+=======================================================+
-| `Rmin`                           | Minimum radius in `x` and `y` [m]                     |
-+----------------------------------+-------------------------------------------------------+
-| `Rmax`                           | Maximum radius in `x` and `y` [m]                     |
-+----------------------------------+-------------------------------------------------------+
-
-* No variation in `z`, `xp`, `yp`, `t`, `s` and total energy. Only central values.
++--------------------+------------------------------------+------------------+
+| **Variable**       | **Description**                    | **Range**        |
++====================+====================================+==================+
+| `Rmin`             | Minimum radius in `x` and `y` [m]  | 0 - inf.         |
++--------------------+------------------------------------+------------------+
+| `Rmax`             | Maximum radius in `x` and `y` [m]  | 0 - inf.         |
++--------------------+------------------------------------+------------------+
+| `Rpmin`            | Minimum radius in `rp` (optional)  | 0 - 1            |
++--------------------+------------------------------------+------------------+
+| `Rpmax`            | Maximum radius in `rp` (optional)  | 0 - 1            |
++--------------------+------------------------------------+------------------+
+
+* No variation in `z`, `t`, `s` and total energy. Only central values.
+* `Rpmin` and `Rpmax` are the radial component of the unit momentum vector, akin to `xp` and `yp`
+  and should be < 1.
+* `Rpmin` must be less than or equal to `Rpmax`.
+* If only one single, non-zero, value of `Rp` is required, both `Rpmin` and `Rpmax`
+  should be set to the same value.
 
 
 eshell

manual/source/version_history.rst

@@ -46,6 +46,8 @@ New Features
   other distributions with each type of Gaussian have been added: :code:`gauss-slow-ext`,
   :code:`gauss-slow-ext`, :code:`gauss-twiss-slow-ext` where the Gaussian part is first generated
   and the slow extraction sweep added on top.
+* The :code:`ring` distribution can now make a circular ring also in `rp` (i.e. in `xp`, `yp`). By
+  default the old behaviour is retained and no angle is produced.
 
 **Fields**
 
@@ -163,7 +165,7 @@ Output Class Versions
 +===================================+=============+=================+=================+
 | BDSOutputROOTEventAperture        | N           | 1               | 1               |
 +-----------------------------------+-------------+-----------------+-----------------+
-| BDSOutputROOTEventBeam            | N           | 6               | 6               |
+| BDSOutputROOTEventBeam            | N           | 6               | 7               |
 +-----------------------------------+-------------+-----------------+-----------------+
 | BDSOutputROOTEventCavityInfo      | N           | 1               | 1               |
 +-----------------------------------+-------------+-----------------+-----------------+

parser/beam.cc

@@ -288,8 +288,10 @@ void Beam::PublishMembers()
   publish("shellYpWidth",&Beam::shellYpWidth);
 
   // options for beam distrType="ring"
-  publish("Rmin",&Beam::Rmin);
-  publish("Rmax",&Beam::Rmax);
+  publish("Rmin", &Beam::Rmin);
+  publish("Rmax", &Beam::Rmax);
+  publish("Rpmin", &Beam::Rpmin);
+  publish("Rpmax", &Beam::Rpmax);
 
   // options for beam distrType="halo"
   publish("haloNSigmaXInner",      &Beam::haloNSigmaXInner);

parser/beamBase.cc

@@ -114,6 +114,8 @@ BeamBase::BeamBase()
 
   Rmin = 0.0;
   Rmax = 0.0;
+  Rpmin = 0.0;
+  Rpmax = 0.0;
 
   haloNSigmaXInner      = 0.0;
   haloNSigmaXOuter      = 1e9;
@@ -162,5 +164,4 @@ BeamBase::BeamBase()
   polarization1 = 0;
   polarization2 = 0;
   polarization3 = 0;
-
 }

parser/beamBase.h

@@ -118,6 +118,7 @@ namespace GMAD
 
       ///@{ for the ring beam distribution
       double Rmin, Rmax;
+      double Rpmin, Rpmax;
       ///@}
 
       ///@{ for the halo distribution
@@ -169,15 +170,13 @@ namespace GMAD
       double dPStart;
       double dPStop;
       /// @}
-
 
       /// A list of all the keys that have been set in this instance.
       std::vector<std::string> setKeys;
 
       double polarization1;
       double polarization2;
       double polarization3;
-
     };
 }
 

src/BDSBunchRing.cc

@@ -20,6 +20,7 @@ along with BDSIM.  If not, see <http://www.gnu.org/licenses/>.
 #include "BDSDebug.hh"
 #include "BDSException.hh"
 #include "BDSParticleCoordsFull.hh"
+#include "BDSUtilities.hh"
 
 #include "parser/beam.h"
 
@@ -33,7 +34,11 @@ along with BDSIM.  If not, see <http://www.gnu.org/licenses/>.
 BDSBunchRing::BDSBunchRing():
   BDSBunch("ring"),
   rMin(0),
-  rMax(0)
+  rMax(0),
+  rpMin(0),
+  rpMax(0),
+  nonZeroRp(false),
+  generateRp(false)
 {;}
 
 BDSBunchRing::~BDSBunchRing()
@@ -48,17 +53,31 @@ void BDSBunchRing::SetOptions(const BDSParticleDefinition* beamParticle,
   BDSBunch::SetOptions(beamParticle, beam, distrType, beamlineTransformIn, beamlineSIn);
   rMin = beam.Rmin * CLHEP::m;
   rMax = beam.Rmax * CLHEP::m;
+  rpMin = beam.Rpmin * CLHEP::rad;
+  rpMax = beam.Rpmax * CLHEP::rad;
+  nonZeroRp = rpMin > 0 || rpMax > 0;
+  generateRp = rpMin != rpMax;
 }
 
 void BDSBunchRing::CheckParameters()
 {
   BDSBunch::CheckParameters();
   if (rMin < 0)
-    {throw BDSException(__METHOD_NAME__, "rMin: " + std::to_string(rMin) + " < 0");}
+    {throw BDSException(__METHOD_NAME__, "Rmin: " + std::to_string(rMin) + " < 0");}
   if (rMax < 0)
-    {throw BDSException(__METHOD_NAME__, "rMax: " + std::to_string(rMin) + " < 0");}
+    {throw BDSException(__METHOD_NAME__, "Rmax: " + std::to_string(rMin) + " < 0");}
   if (rMax <= rMin)
-    {throw BDSException(__METHOD_NAME__, "rMax: " + std::to_string(rMax) + " < rMin: " + std::to_string(rMin));}
+    {throw BDSException(__METHOD_NAME__, "Rmax: " + std::to_string(rMax) + " < Rmin: " + std::to_string(rMin));}
+  if (rpMin < 0)
+    {throw BDSException(__METHOD_NAME__, "Rpmin: " + std::to_string(rpMin) + " < 0");}
+  if (rpMax < 0)
+    {throw BDSException(__METHOD_NAME__, "Rpmax: " + std::to_string(rpMax) + " < 0");}
+  if (rpMax < rpMin)
+    {throw BDSException(__METHOD_NAME__, "Rpmax: " + std::to_string(rpMin) + " < Rpmin: " + std::to_string(rpMin));}
+  if (rpMin > 1)
+    {throw BDSException(__METHOD_NAME__, "Rpmin: " + std::to_string(rpMin) + " > 1");}
+  if (rpMax > 1)
+    {throw BDSException(__METHOD_NAME__, "Rpmax: " + std::to_string(rpMin) + " > 1");}
 }
 
 BDSParticleCoordsFull BDSBunchRing::GetNextParticleLocal()
@@ -67,6 +86,20 @@ BDSParticleCoordsFull BDSBunchRing::GetNextParticleLocal()
   G4double phi = 2 * CLHEP::pi * G4RandFlat::shoot();
   G4double x   = X0 + r * std::sin(phi);
   G4double y   = Y0 + r * std::cos(phi);
-
-  return BDSParticleCoordsFull(x,y,Z0,Xp0,Yp0,Zp0,T0,S0,E0,/*weight=*/1.0);
+  G4double xp = Xp0;
+  G4double yp = Yp0;
+  if (generateRp)
+    {
+      G4double rp   = std::sqrt(G4RandFlat::shoot(std::pow(rpMin,2), std::pow(rpMax,2)));
+      xp += rp * std::sin(phi);
+      yp += rp * std::cos(phi);
+    }
+  else if (nonZeroRp)
+    {
+      G4double rp = BDS::IsFinite(rpMin) ? rpMin : rpMax;
+      xp += rp * std::sin(phi);
+      yp += rp * std::cos(phi);
+    }
+  G4double zp = CalculateZp(xp,yp,Zp0);
+  return BDSParticleCoordsFull(x,y,Z0,xp,yp,zp,T0,S0,E0,/*weight=*/1.0);
 }