Skip to content

Muon Cooler solenoid rotations #111

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 14 commits into
base: develop
Choose a base branch
from
Open

Muon Cooler solenoid rotations #111

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
ce8c880
added hard coded rotations
kamathrohan Sep 29, 2025
39420b7
added tilts
kamathrohan Sep 30, 2025
5e0030a
added new gmad file
kamathrohan Oct 2, 2025
606b051
removed rogue cout from debugging
kamathrohan Oct 6, 2025
027a818
added offsets
kamathrohan Oct 9, 2025
cfdb52f
rewrote offsets to pass them globally
kamathrohan Oct 13, 2025
53ab05a
added offsets to volumes
kamathrohan Oct 15, 2025
f83b754
added rotations to the volumes
kamathrohan Oct 16, 2025
3c62a43
added changes to examples
kamathrohan Oct 20, 2025
5b54562
changed gmad file in one more location
kamathrohan Oct 20, 2025
486764b
changed order - cosmetic change
kamathrohan Oct 20, 2025
10c740d
another change to the components
kamathrohan Oct 20, 2025
9a66aa0
added changes to manual
kamathrohan Oct 21, 2025
4adcba1
added rotation convention to the manual
kamathrohan Oct 22, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion examples/features/components/muoncooler.gmad
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,14 @@ cooldef1: coolingchannel,
nCoils=6,
coilInnerRadius = {0.3*m},
coilRadialThickness = {0.1*m},
coilLengthZ = {0.1333333*m},
coilOffsetX = {0.0*m},
coilOffsetY = {0.0*m},
coilOffsetZ ={-1.3153, -0.6847, -0.3153, 0.3153, 0.6847, 1.3153},

coilLengthZ = {0.1333333*m},
coilTiltX ={0.0*degrees},
coilTiltY = {0.0*degrees},
coilTiltZ = {0.0*degrees},
coilCurrent = {4900000, -4900000, 4900000, -4900000, 4900000, -4900000},
coilMaterial = {"G4_Cu"},
onAxisTolerance=1e-7,
Expand Down
5 changes: 5 additions & 0 deletions examples/muoncooling/muoncooler.gmad
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,12 @@ cooldef1: coolingchannel,
coilInnerRadius = {0.3*m},
coilRadialThickness = {0.1*m},
coilLengthZ = {0.1333333*m},
coilOffsetX = {0.0*m},
coilOffsetY = {0.0*m},
coilOffsetZ ={-1.3153, -0.6847, -0.3153, 0.3153, 0.6847, 1.3153},
coilTiltX ={0.0*degrees},
coilTiltY = {0.0*degrees},
coilTiltZ = {0.0*degrees},
coilCurrent = {4900000, -4900000, 4900000, -4900000, 4900000, -4900000},
coilMaterial = {"G4_Cu"},
onAxisTolerance=1e-7,
Expand Down
6 changes: 6 additions & 0 deletions include/BDSFieldMagSolenoidBlock.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,9 @@ public:
G4double innerRadiusIn,
G4double radialThicknessIn,
G4double fullLengthZIn,
G4double tiltXIn,
G4double tiltYIn,
G4double tiltZIn,
G4double toleranceIn,
G4int nSheetsIn);
/// Alternative constructor for field factory that uses "field" (i.e. B) strength
Expand All @@ -71,6 +74,9 @@ private:
G4double fullLengthZ;
G4double B0;
G4double I;
G4double tiltX;
G4double tiltY;
G4double tiltZ;
G4double coilTolerance;
G4int nSheetsBlock;
G4double currentDensity;
Expand Down
6 changes: 6 additions & 0 deletions include/BDSFieldMagSolenoidSheet.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,9 @@ public:
G4bool strengthIsCurrent,
G4double sheetRadius,
G4double fullLength,
G4double tiltX,
G4double tiltY,
G4double tiltZ,
G4double toleranceIn = 0.0);
virtual ~BDSFieldMagSolenoidSheet(){;}

Expand All @@ -79,6 +82,9 @@ private:
G4double I;
G4double spatialLimit;
G4double normalisation;
G4double rotateX;
G4double rotateY;
G4double rotateZ;
G4double coilTolerance;
};

Expand Down
5 changes: 5 additions & 0 deletions include/BDSMuonCoolerStructs.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,12 @@ namespace BDS
G4double radialThickness;
G4double fullLengthZ;
G4double current;
G4double offsetX;
G4double offsetY;
G4double offsetZ;
G4double tiltX;
G4double tiltY;
G4double tiltZ;
G4Material* material;
G4double onAxisTolerance;
G4int nSheets;
Expand Down
18 changes: 18 additions & 0 deletions manual/source/model_construction.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2311,6 +2311,9 @@ Parameters for these components can be specified as either:
- For dipoles, two models exist currently: `dipole` and `dipoleenge`. The `dipole` model is a simple hard-edge dipole field, while the `dipoleenge` model includes Enge-type fringe fields and follows the treatment outlined in: Muratori, B.D. et al (2015) ‘Analytical expressions for fringe fields in multipole magnets’, *Physical Review Special Topics - Accelerators and Beams*, 18(6). https://doi.org/10.1103/physrevstab.18.064001
- For the RF cavities, a simple RF pillbox (`rfpillbox`) model has been implemented.

**Rotations and Offsets**

Rotations follow the right-hand rule using axis-angle representation applied in XYZ order about the X, Y, and Z axes respectively. The reference frame assumes the solenoid is initially centered at the origin with its axis aligned along the z-direction. The position offsets then translate this solenoid center from the origin to an arbitrary location in the global coordinate system of the cooling channel.
Copy link
Contributor

@lnevay lnevay Oct 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"offsets then translate(s)".. just translates. I think not in the global frame but in the coordinate frame of the cooling channel element.


**Table of Parameters**

Expand All @@ -2328,9 +2331,24 @@ Parameters for these components can be specified as either:
+------------------------------+-------------------------------+--------------+
| `coilLengthZ` | Lengths of coils along Z [m] | List[Float] |
+------------------------------+-------------------------------+--------------+
| `coilOffsetX` | X-positions of coil centers | List[Float] |
| | [m] | |
+------------------------------+-------------------------------+--------------+
| `coilOffsetY` | Y-positions of coil centers | List[Float] |
| | [m] | |
+------------------------------+-------------------------------+--------------+
| `coilOffsetZ` | Z-positions of coil centers | List[Float] |
| | [m] | |
+------------------------------+-------------------------------+--------------+
| `coilTiltX` | Tilt angles of coil centers | List[Float] |
| | about X axis [rad] | |
+------------------------------+-------------------------------+--------------+
| `coilTiltY` | Tilt angles of coil centers | List[Float] |
| | about Y axis [rad] | |
+------------------------------+-------------------------------+--------------+
| `coilTiltZ` | Tilt angles of coil centers | List[Float] |
| | about Z axis [rad] | |
+------------------------------+-------------------------------+--------------+
| `coilCurrent` | Currents in [A] (sheet model) | List[Float] |
| | or densities | |
| | [A/m^2] (block model) | |
Expand Down
15 changes: 15 additions & 0 deletions parser/coolingchannel.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,12 @@ void CoolingChannel::clear()
coilRadialThickness.clear();
coilLengthZ.clear();
coilCurrent.clear();
coilOffsetX.clear();
coilOffsetY.clear();
coilOffsetZ.clear();
coilTiltX.clear();
coilTiltY.clear();
coilTiltZ.clear();
coilMaterial.clear();
mirrorCoils = false;
onAxisTolerance = 0;
Expand Down Expand Up @@ -99,7 +104,12 @@ void CoolingChannel::PublishMembers()
publish("coilRadialThickness", &CoolingChannel::coilRadialThickness);
publish("coilLengthZ", &CoolingChannel::coilLengthZ);
publish("coilCurrent", &CoolingChannel::coilCurrent);
publish("coilOffsetX", &CoolingChannel::coilOffsetX);
publish("coilOffsetY", &CoolingChannel::coilOffsetY);
publish("coilOffsetZ", &CoolingChannel::coilOffsetZ);
publish("coilTiltX", &CoolingChannel::coilTiltX);
publish("coilTiltY", &CoolingChannel::coilTiltY);
publish("coilTiltZ", &CoolingChannel::coilTiltZ);
publish("coilMaterial", &CoolingChannel::coilMaterial);
publish("mirrorCoils", &CoolingChannel::mirrorCoils);
publish("onAxisTolerance", &CoolingChannel::onAxisTolerance);
Expand Down Expand Up @@ -149,7 +159,12 @@ void CoolingChannel::PublishMembers()
attribute_map_list_double["coilRadialThickness"] = &coilRadialThickness;
attribute_map_list_double["coilLengthZ"] = &coilLengthZ;
attribute_map_list_double["coilCurrent"] = &coilCurrent;
attribute_map_list_double["coilOffsetX"] = &coilOffsetX;
attribute_map_list_double["coilOffsetY"] = &coilOffsetY;
attribute_map_list_double["coilOffsetZ"] = &coilOffsetZ;
attribute_map_list_double["coilTiltX"] = &coilTiltX;
attribute_map_list_double["coilTiltY"] = &coilTiltY;
attribute_map_list_double["coilTiltZ"] = &coilTiltZ;
attribute_map_list_string["coilMaterial"] = &coilMaterial;
attribute_map_list_double["dipoleAperture"] = &dipoleAperture;
attribute_map_list_double["dipoleLengthZ"] = &dipoleLengthZ;
Expand Down
5 changes: 5 additions & 0 deletions parser/coolingchannel.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,12 @@ namespace GMAD
std::list<double> coilRadialThickness;
std::list<double> coilLengthZ;
std::list<double> coilCurrent;
std::list<double> coilOffsetX;
std::list<double> coilOffsetY;
std::list<double> coilOffsetZ;
std::list<double> coilTiltX;
std::list<double> coilTiltY;
std::list<double> coilTiltZ;
std::list<std::string> coilMaterial;
bool mirrorCoils;
double onAxisTolerance;
Expand Down
5 changes: 5 additions & 0 deletions parser/test/coolingchannel.gmad
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,12 @@ extracool: coolingchannel, surroundingMaterial="G4_Galactic",
coilInnerRadius = {0.3*m},
coilRadialThickness = {0.1*m},
coilLengthZ = {0.1333333*m},
coilOffsetX = {0.0*m},
coilOffsetY = {0.0*m},
coilOffsetZ ={-1.3153, -0.6847, -0.3153, 0.3153, 0.6847, 1.3153},
coilTiltX ={0.0*degrees},
coilTiltY = {0.0*degrees},
coilTiltZ = {0.0*degrees},
coilCurrent = {4900000, -4900000, 4900000, -4900000, 4900000, -4900000},
coilMaterial = {"G4_Cu"},
onAxisTolerance=1e-7,
Expand Down
10 changes: 8 additions & 2 deletions src/BDSFieldEMMuonCooler.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,9 +63,12 @@ void BDSFieldEMMuonCooler::BuildMagnets(const BDSFieldInfoExtraMuonCooler* info)
ci.innerRadius,
ci.radialThickness,
ci.fullLengthZ,
ci.tiltX,
ci.tiltY,
ci.tiltZ,
ci.onAxisTolerance,
ci.nSheets));
fieldOffsets.emplace_back(0,0,ci.offsetZ);
fieldOffsets.emplace_back(ci.offsetX, ci.offsetY, ci.offsetZ);
}
coilField = new BDSFieldMagVectorSum(fields, fieldOffsets);
break;
Expand All @@ -81,8 +84,11 @@ void BDSFieldEMMuonCooler::BuildMagnets(const BDSFieldInfoExtraMuonCooler* info)
true,
ci.innerRadius + 0.5*ci.radialThickness,
ci.fullLengthZ,
ci.tiltX,
ci.tiltY,
ci.tiltZ,
ci.onAxisTolerance));
fieldOffsets.emplace_back(0,0,ci.offsetZ);
fieldOffsets.emplace_back(ci.offsetX, ci.offsetY, ci.offsetZ);
}
coilField = new BDSFieldMagVectorSum(fields, fieldOffsets);
break;
Expand Down
11 changes: 10 additions & 1 deletion src/BDSFieldMagSolenoidBlock.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ along with BDSIM. If not, see <http://www.gnu.org/licenses/>.

BDSFieldMagSolenoidBlock::BDSFieldMagSolenoidBlock(BDSMagnetStrength const* strength,
G4double innerRadiusIn):
BDSFieldMagSolenoidBlock((*strength)["field"], false, innerRadiusIn, (*strength)["coilRadialThickness"], (*strength)["length"], 0, 1)
BDSFieldMagSolenoidBlock((*strength)["field"], false, innerRadiusIn, (*strength)["coilRadialThickness"], (*strength)["length"], 0.0, 0.0, 0.0, 0, 1)
{;}


Expand All @@ -44,13 +44,19 @@ BDSFieldMagSolenoidBlock::BDSFieldMagSolenoidBlock(G4double strength,
G4double innerRadiusIn,
G4double radialThicknessIn,
G4double fullLengthZIn,
G4double tiltXIn,
G4double tiltYIn,
G4double tiltZIn,
G4double toleranceIn,
G4int nSheetsIn):
a(innerRadiusIn),
radialThickness(radialThicknessIn),
fullLengthZ(fullLengthZIn),
B0(0),
I(0),
tiltX(tiltXIn),
tiltY(tiltYIn),
tiltZ(tiltZIn),
coilTolerance(toleranceIn/(nSheetsIn*2)), //double the tolerance for each sheet
nSheetsBlock(nSheetsIn)
{
Expand Down Expand Up @@ -86,6 +92,9 @@ G4ThreeVector BDSFieldMagSolenoidBlock::GetField(const G4ThreeVector& position,
true,
a + (sheet * dr) + dr / 2,
fullLengthZ,
tiltX,
Copy link
Collaborator

@stewartboogert stewartboogert Oct 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fix indenting

tiltY,
tiltZ,
coilTolerance);
sheetField = field->GetField(position);
if (sheetField == G4ThreeVector(0, 0, 0))
Expand Down
44 changes: 35 additions & 9 deletions src/BDSFieldMagSolenoidSheet.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,21 +32,27 @@ along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
#include <cmath>

BDSFieldMagSolenoidSheet::BDSFieldMagSolenoidSheet(BDSMagnetStrength const* strength,
G4double radiusIn, G4double toleranceIn):
BDSFieldMagSolenoidSheet((*strength)["field"], false, radiusIn, (*strength)["length"], toleranceIn)
G4double radiusIn,G4double toleranceIn):
BDSFieldMagSolenoidSheet((*strength)["field"], false, radiusIn, (*strength)["length"], 0.0, 0.0, 0.0, toleranceIn)
{;}

BDSFieldMagSolenoidSheet::BDSFieldMagSolenoidSheet(G4double strength,
G4bool strengthIsCurrent,
G4double sheetRadius,
G4double fullLength,
G4double tiltX,
Copy link
Collaborator

@stewartboogert stewartboogert Oct 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fix indenting

G4double tiltY,
G4double tiltZ,
G4double toleranceIn):
a(sheetRadius),
halfLength(0.5*fullLength),
B0(0.0),
I(0.0),
spatialLimit(std::min(1e-5*sheetRadius, 1e-5*fullLength)),
normalisation(1.0) ,
rotateX(tiltX),
Copy link
Collaborator

@stewartboogert stewartboogert Oct 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Try to be more consistent in variable names. Could be rotateX(rotateXIn)

Copy link
Contributor

@lnevay lnevay Oct 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agree - should be consistent.

rotateY(tiltY),
rotateZ(tiltZ),
coilTolerance(toleranceIn)
{
finiteStrength = BDS::IsFinite(std::abs(strength));
Expand All @@ -72,10 +78,22 @@ BDSFieldMagSolenoidSheet::BDSFieldMagSolenoidSheet(G4double strength,
G4ThreeVector BDSFieldMagSolenoidSheet::GetField(const G4ThreeVector& position,
const G4double /*t*/) const
{
G4double z = position.z();
G4double rho = position.perp();
G4double phi = position.phi(); // angle about z axis

// Rotation angles - to be moved to struct later
//G4double rotationX = CLHEP::pi/4; // 45 degrees to put solenoid along z axis
//G4double rotationY = 0.0;
//G4double rotationZ = 0.0;

// Transform position from global to local coordinates (inverse rotation)
G4ThreeVector localPosition = position;

// Apply inverse rotations in reverse order (Z -> Y -> X)
// to transform from global frame to solenoid's local frame
localPosition.rotateZ(-rotateZ);
Copy link
Collaborator

@stewartboogert stewartboogert Oct 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not formulate the transformation and just call the inverse method?

localPosition.rotateY(-rotateY);
Copy link
Contributor

@lnevay lnevay Oct 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Like Stewart mentioned, I would perhaps use a G4RotationMatrix. You can prepare it in the constructor and the inverse as well. Then you can just multiply by the appropriate one each time.

localPosition.rotateX(-rotateX);
G4double z = localPosition.z();
G4double rho = localPosition.perp();
G4double phi = localPosition.phi(); // angle about z axis
// check if close to current source - function not well-behaved at exactly the rho of
// the current source or at the boundary of +- halfLength
if (std::abs(rho - a) < spatialLimit && (std::abs(z) < halfLength+2*spatialLimit))
Expand Down Expand Up @@ -129,9 +147,17 @@ G4ThreeVector BDSFieldMagSolenoidSheet::GetField(const G4ThreeVector& position,
}
// we have to be consistent with the phi we calculated at the beginning,
// so unit rho is in the x direction.
G4ThreeVector result = G4ThreeVector(Brho,0,Bz)* normalisation;
result = result.rotateZ(phi);
return result;
G4ThreeVector localField = G4ThreeVector(Brho,0,Bz) * normalisation;
localField = localField.rotateZ(phi);

// Transform field from local back to global coordinates (forward rotation)
// Apply rotations in forward order (X -> Y -> Z)
G4ThreeVector globalField = localField;
globalField.rotateX(rotateX);
globalField.rotateY(rotateY);
globalField.rotateZ(rotateZ);

return globalField;
Copy link
Contributor

@lnevay lnevay Oct 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please can you change this to some other name. The "globalField" is used in BDSIM to mean in the global frame, whereas this is still in the local frame of the element. You usually use one mega-cooling channel but conceptually it's 'local' as the global transform for the beamline is done in a wrapper class. How about unrotatedField and localField.

}

G4double BDSFieldMagSolenoidSheet::OnAxisBz(G4double zp,
Expand Down
10 changes: 7 additions & 3 deletions src/BDSMuonCooler.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,9 +148,13 @@ void BDSMuonCooler::BuildCoils()
auto coilLV = new G4LogicalVolume(coilSolid, info.material, baseName + "_lv");
RegisterLogicalVolume(coilLV);
coilLV->SetVisAttributes(coilVises[info.material]);

auto coilPV = new G4PVPlacement(nullptr,
G4ThreeVector(0,0,info.offsetZ),
G4RotationMatrix* rmCoil = new G4RotationMatrix();
rmCoil->rotateX(info.tiltX);
rmCoil->rotateY(info.tiltY);
rmCoil->rotateZ(info.tiltZ);
RegisterRotationMatrix(rmCoil);
auto coilPV = new G4PVPlacement(rmCoil,
G4ThreeVector(info.offsetX,info.offsetY,info.offsetZ),
coilLV,
baseName + "_pv",
containerLogicalVolume,
Expand Down
Loading
Loading