20#include "BDSAcceleratorComponent.hh"
21#include "BDSBeamline.hh"
22#include "BDSBeamlineElement.hh"
23#include "BDSException.hh"
24#include "BDSExtentGlobal.hh"
25#include "BDSGlobalConstants.hh"
27#include "BDSOutput.hh"
28#include "BDSSamplerPlane.hh"
29#include "BDSSimpleComponent.hh"
30#include "BDSTiltOffset.hh"
31#include "BDSTransform3D.hh"
32#include "BDSUtilities.hh"
33#include "BDSWarning.hh"
36#include "G4RotationMatrix.hh"
37#include "G4ThreeVector.hh"
38#include "G4Transform3D.hh"
40#include "CLHEP/Vector/AxisAngle.h"
50BDSBeamline::BDSBeamline(
const G4ThreeVector& initialGlobalPosition,
51 G4RotationMatrix* initialGlobalRotation,
58 previousReferencePositionEnd(initialGlobalPosition),
59 previousSPositionEnd(sInitial),
60 transformHasJustBeenApplied(false)
67 if (initialGlobalRotation)
78 G4cout << __METHOD_NAME__ <<
"with initial position and rotation" << G4endl;
79 G4cout <<
"Initial position: " << initialGlobalPosition << G4endl;
84BDSBeamline::BDSBeamline(G4Transform3D initialTransform,
87 new G4RotationMatrix(initialTransform.getRotation()),
91BDSBeamline::~BDSBeamline()
104 for (
const auto& element : *
this)
110 G4cout << __METHOD_NAME__ <<
"container size: " <<
sizeof(
beamline) << G4endl;
115std::ostream& operator<< (std::ostream& out,
BDSBeamline const &bl)
117 out <<
"Beamline with " << bl.
size() <<
" elements" << G4endl;
118 out <<
"Total arc length: " << bl.
totalArcLength <<
" mm" << G4endl;
130 G4String samplerName = samplerInfo ? samplerInfo->name :
"no_sampler_name_given";
131 throw BDSException(__METHOD_NAME__,
"invalid accelerator component " + samplerName);
137 G4String samplerName = samplerInfo->name;
140 G4cerr << __METHOD_NAME__ <<
"invalid sampler name \"" << samplerName <<
"\"" << G4endl;
153 G4int sizeLine = (G4int)line->size();
154 for (G4int i = 0; i < sizeLine; ++i)
170 if (samplerInfo->samplerType == BDSSamplerType::cylinder)
176 G4ThreeVector aaMidAxis;
183 aaMidAxis = aaFinish.axis();
184 aaMidAngle = 0.5 * aaFinish.delta();
188 aaMidAxis = aaStart.axis();
189 aaMidAngle = 0.5 * aaStart.delta();
193 aaMidAxis = (aaFinish.axis() + aaStart.axis()) / 2.0;
194 aaMidAngle = (aaFinish.delta() + aaStart.delta()) / 2.0;
196 auto aaCSampler = CLHEP::HepAxisAngle(aaMidAxis, aaMidAngle);
197 G4RotationMatrix rmCSampler = G4RotationMatrix(aaCSampler);
198 G4Transform3D trCSampler(rmCSampler, midRefPosition);
214 G4cout << G4endl << __METHOD_NAME__ <<
"adding component to beamline and calculating coordinates" << G4endl;
215 G4cout <<
"component name: " << component->
GetName() << G4endl;
231 G4double angle = component->
GetAngle();
234 G4bool hasFiniteTilt, hasFiniteOffset;
235 G4ThreeVector offset;
244 hasFiniteTilt =
false;
245 hasFiniteOffset =
false;
246 offset = G4ThreeVector();
251 G4bool hasFinitePlacementOffset =
BDS::IsFinite(placementOffset);
255 G4cout <<
"chord length " << chordLength <<
" mm" << G4endl;
256 G4cout <<
"angle " << angle <<
" rad" << G4endl;
258 {G4cout <<
"tilt offsetX offsetY " << *tiltOffset <<
" rad mm mm " << G4endl;}
260 {G4cout <<
"no tilt offset" << G4endl;}
261 G4cout <<
"has finite length " << hasFiniteLength << G4endl;
262 G4cout <<
"has finite angle " << hasFiniteAngle << G4endl;
263 G4cout <<
"has finite tilt " << hasFiniteTilt << G4endl;
264 G4cout <<
"has finite offset " << hasFiniteOffset << G4endl;
265 G4cout <<
"extent positive " << eP << G4endl;
266 G4cout <<
"extent negative " << eN << G4endl;
267 G4cout <<
"object placement offset " << placementOffset << G4endl;
268 G4cout <<
"has finite placement offset " << hasFinitePlacementOffset << G4endl;
269 G4cout <<
"output face normal " << oFNormal << G4endl;
274 if (!
empty() && (component->
GetType() !=
"drift") && (component->
GetType() !=
"thinmultipole"))
276 G4bool keepGoing =
true;
277 G4bool checkFaces =
true;
278 G4double zSeparation = 0;
281 G4ThreeVector iFNormal;
282 G4String clasherName =
"Unknown";
285 if (inspectedElement)
287 if ((inspectedElement->
GetType() ==
"drift")||(inspectedElement->
GetType() ==
"thinmultipole"))
306 G4cout <<
"input face normal " << iFNormal << G4endl;
315 G4bool willIntersect =
BDS::WillIntersect(iFNormal, oFNormal, zSeparation, extIF, extOF);
318 G4cout << __METHOD_NAME__ <<
"Error - angled faces of objects will cause overlap in beam line geometry" << G4endl;
319 G4cout <<
"\"" << component->
GetName() <<
"\" will overlap with \""
320 << clasherName <<
"\"" << G4endl;
330 G4RotationMatrix* referenceRotationStart;
340 G4RotationMatrix* referenceRotationMiddle =
new G4RotationMatrix(*referenceRotationStart);
341 G4RotationMatrix* referenceRotationEnd =
new G4RotationMatrix(*referenceRotationStart);
351 G4ThreeVector rotationAxisOfBend = G4ThreeVector(0,1,0);
352 G4ThreeVector rotationAxisOfBendEnd = rotationAxisOfBend;
355 G4double tilt = tiltOffset->
GetTilt();
356 G4RotationMatrix rotationAxisRM = G4RotationMatrix();
357 rotationAxisRM.rotateZ(tilt);
358 rotationAxisOfBend.transform(rotationAxisRM);
359 rotationAxisOfBendEnd.transform(rotationAxisRM);
366 G4RotationMatrix* rotationStart =
new G4RotationMatrix(*referenceRotationStart);
367 G4RotationMatrix* rotationMiddle =
new G4RotationMatrix(*referenceRotationMiddle);
368 G4RotationMatrix* rotationEnd =
new G4RotationMatrix(*referenceRotationEnd);
372 G4double tilt = tiltOffset->
GetTilt();
376 G4ThreeVector unitZ = G4ThreeVector(0,0,1);
377 rotationStart ->rotate(tilt, unitZ.transform(*referenceRotationStart));
378 unitZ = G4ThreeVector(0,0,1);
379 rotationMiddle->rotate(tilt, unitZ.transform(*referenceRotationMiddle));
380 unitZ = G4ThreeVector(0,0,1);
381 rotationEnd ->rotate(tilt, unitZ.transform(*referenceRotationEnd));
399 if (samplerType != BDSSamplerType::none)
408 G4ThreeVector componentGap = pad.transform(*previousReferenceRotationEnd2);
412 G4ThreeVector referencePositionStart, referencePositionMiddle, referencePositionEnd;
418 G4ThreeVector md = G4ThreeVector(0, 0, 0.5 * chordLength);
419 md.transform(*referenceRotationMiddle);
420 referencePositionMiddle = referencePositionStart + md;
423 G4ThreeVector delta = G4ThreeVector(0, 0, chordLength).transform(*referenceRotationMiddle);
424 referencePositionEnd = referencePositionStart + delta;
436 G4ThreeVector positionStart, positionMiddle, positionEnd;
437 if (hasFiniteOffset || hasFinitePlacementOffset)
439 if (hasFiniteOffset && hasFiniteAngle)
441 G4String name = component->
GetName();
442 G4String message =
"element has x offset, but this will cause geometry overlaps: " + name
443 +
" - omitting x offset";
444 BDS::Warning(__METHOD_NAME__, message);
449 G4ThreeVector total = offset + placementOffset;
450 G4ThreeVector displacement = total.transform(*referenceRotationMiddle);
451 positionStart = referencePositionStart + displacement;
452 positionMiddle = referencePositionMiddle + displacement;
453 positionEnd = referencePositionEnd + displacement;
457 positionStart = referencePositionStart;
458 positionMiddle = referencePositionMiddle;
459 positionEnd = referencePositionEnd;
478 sMaximum += arcLength;
485 G4cout <<
"calculated coordinates in mm and rad are " << G4endl;
486 G4cout <<
"reference position start: " << referencePositionStart << G4endl;
487 G4cout <<
"reference position middle: " << referencePositionMiddle << G4endl;
488 G4cout <<
"reference position end: " << referencePositionEnd << G4endl;
489 G4cout <<
"reference rotation start: " << *referenceRotationStart;
490 G4cout <<
"reference rotation middle: " << *referenceRotationMiddle;
491 G4cout <<
"reference rotation end: " << *referenceRotationEnd;
492 G4cout <<
"position start: " << positionStart << G4endl;
493 G4cout <<
"position middle: " << positionMiddle << G4endl;
494 G4cout <<
"position end: " << positionEnd << G4endl;
495 G4cout <<
"rotation start: " << *rotationStart;
496 G4cout <<
"rotation middle: " << *rotationMiddle;
497 G4cout <<
"rotation end: " << *rotationEnd;
513 referencePositionStart,
514 referencePositionMiddle,
515 referencePositionEnd,
516 referenceRotationStart,
517 referenceRotationMiddle,
518 referenceRotationEnd,
533 sEnd.push_back(sPositionEnd);
544 G4double dx = component->dx;
545 G4double dy = component->dy;
546 G4double dz = component->dz;
551 G4cerr << __METHOD_NAME__ <<
"Caution: Transform3d: " << component->
GetName() << G4endl;
552 G4cerr << __METHOD_NAME__ <<
"dz = " << dz <<
" < 0 -> could overlap previous element" << G4endl;
570 G4RotationMatrix trRotInverse = component->rotationMatrix.inverse();
571 (*previousReferenceRotationEnd) *= trRotInverse;
577 {
throw BDSException(__METHOD_NAME__,
"invalid BDSBeamlineElement");}
579 {
throw BDSException(__METHOD_NAME__,
"invalid BDSAcceleratorComponent");}
596 for (
int i=0; i<3; i++)
602 G4int* indexOfFoundElement)
const
608 G4String msg =
"s position " + std::to_string(s/CLHEP::m) +
" m is beyond length of accelerator (";
609 msg += std::to_string(
totalArcLength/CLHEP::m) +
" m)\nReturning identify transform";
610 BDS::Warning(__METHOD_NAME__, msg);
611 return G4Transform3D();
617 G4cout << __METHOD_NAME__ << G4endl;
618 G4cout <<
"S position requested: " << s << G4endl;
619 G4cout <<
"Index: " << indexOfFoundElement << G4endl;
620 G4cout <<
"Element: " << *element << G4endl;
632 G4double dS = s - element->GetSPositionMiddle();
633 G4double localZ = dS * (chordLength / arcLength);
634 G4double angle = component->
GetAngle();
635 G4RotationMatrix rotation;
636 G4RotationMatrix* rotMiddle = element->GetReferenceRotationMiddle();
642 G4ThreeVector localUnitY = G4ThreeVector(0,1,0);
643 localUnitY.transform(*(element->GetReferenceRotationStart()));
645 G4double partialAngle = angle * std::fabs(( (0.5*arcLength + dS) / arcLength));
646 rotation = G4RotationMatrix(*(element->GetReferenceRotationStart()));
647 rotation.rotate(partialAngle, localUnitY);
648 dx = localZ*tan(partialAngle);
651 {rotation = G4RotationMatrix(*rotMiddle);}
654 G4ThreeVector dLocal = G4ThreeVector(x + dx, y , localZ);
656 G4cout <<
"Local offset from middle: " << dLocal << G4endl;
659 G4ThreeVector globalPos = element->GetReferencePositionMiddle() + dLocal.transform(*rotMiddle);
661 G4Transform3D result = G4Transform3D(rotation, globalPos);
664 G4cout <<
"Global offset from middle: " << dLocal << G4endl;
665 G4cout <<
"Resultant global position: " << globalPos << G4endl;
671 G4int* indexOfFoundElement)
const
674 auto lower = std::lower_bound(
sEnd.begin(),
sEnd.end(), S);
675 G4int index = G4int(lower -
sEnd.begin());
676 if (indexOfFoundElement)
677 {*indexOfFoundElement = index;}
683 auto lower = std::lower_bound(
sEnd.begin(),
sEnd.end(), S);
685 std::advance(iter, std::distance(
sEnd.begin(), lower));
703 if (index < 1 || index > (G4int)(
beamline.size()-1))
715 return GetNext(G4int(result -
beamline.begin()));
723 if (index < 0 || index > (G4int)(
beamline.size()-2))
744 G4String suffix =
"_" + std::to_string(i);
745 G4String placementName = acceleratorComponentName + suffix;
746 const auto search =
components.find(placementName);
756 std::vector<const BDSBeamlineElement*> candidates;
757 std::for_each(this->
begin(),
762 if (candidates.empty())
766 auto foundItem = std::find_if(candidates.begin(),
769 {return BDS::EndsWith(el->GetPlacementName(), suffix);});
770 return foundItem != candidates.end() ? *foundItem :
nullptr;
774 {
return search->second;}
783 G4cerr << __METHOD_NAME__ <<
"No element named \""
784 << acceleratorComponentName <<
"\" found for placement number "
786 G4cout <<
"Note, this may be because the element is a bend and split into " << G4endl;
787 G4cout <<
"multiple sections with unique names." << G4endl;
802 for (
const auto& point : boundaryPoints)
804 for (
int i = 0; i < 3; ++i)
828 G4ThreeVector elPosStart = element->
GetPositionStart() - G4ThreeVector(0,0,2*pl).transform(*elRotStart);
829 G4ThreeVector positionMiddle = elPosStart - G4ThreeVector(0,0,endPieceLength*0.5).transform(*elRotStart);
830 G4ThreeVector positionStart = elPosStart - G4ThreeVector(0,0,endPieceLength).transform(*elRotStart);
840 new G4RotationMatrix(*elRotStart),
841 new G4RotationMatrix(*elRotStart),
842 new G4RotationMatrix(*elRotStart),
846 new G4RotationMatrix(*elRotStart),
847 new G4RotationMatrix(*elRotStart),
848 new G4RotationMatrix(*elRotStart),
849 elSPosStart - endPieceLength,
850 elSPosStart - 0.5*endPieceLength,
866 G4RotationMatrix* elRotEnd =
new G4RotationMatrix(*(element->
GetRotationMiddle()));
867 G4ThreeVector elPosEnd = element->
GetPositionEnd() + G4ThreeVector(0,0,pl).transform(*elRotEnd);
868 G4ThreeVector positionMiddle = elPosEnd + G4ThreeVector(0,0,endPieceLength*0.5).transform(*elRotEnd);
869 G4ThreeVector positionEnd = elPosEnd + G4ThreeVector(0,0,endPieceLength).transform(*elRotEnd);
872 G4ThreeVector localUnitY = G4ThreeVector(0,1,0).transform(*elRotEnd);
873 elRotEnd->rotate(CLHEP::pi, localUnitY);
884 new G4RotationMatrix(*elRotEnd),
885 new G4RotationMatrix(*elRotEnd),
886 new G4RotationMatrix(*elRotEnd),
890 new G4RotationMatrix(*elRotEnd),
891 new G4RotationMatrix(*elRotEnd),
892 new G4RotationMatrix(*elRotEnd),
894 elSPosEnd + 0.5*endPieceLength,
895 elSPosEnd + endPieceLength,
903 if (index < 0 || index > (G4int)(
beamline.size()-1))
911 std::vector<G4double> sPos;
917 if (sPos.size() == 1)
918 {sPos.push_back(1*CLHEP::m);}
924 return (std::abs(
totalAngle) > 0.99 * 2.0 * CLHEP::pi) and (std::abs(
totalAngle) < 1.01 * 2.0 * CLHEP::pi);
936 std::vector<G4int> result;
939 if (element->
GetType() == type)
940 {result.push_back(element->
GetIndex());}
947 std::vector<G4int> result;
950 G4String type = element->
GetType();
951 if (types.find(type) != types.end())
952 {result.push_back(element->
GetIndex());}
959 std::set<G4String> collimatorTypes = {
"ecol",
"rcol",
"jcol",
"crystalcol",
"element-collimator"};
Abstract class that represents a component of an accelerator.
G4double GetAngle() const
virtual G4String GetName() const
The name of the component without modification.
virtual G4double GetChordLength() const
G4ThreeVector OutputFaceNormal() const
virtual G4double GetArcLength() const
G4ThreeVector InputFaceNormal() const
G4String GetType() const
Get a string describing the type of the component.
A class that holds a fully constructed BDSAcceleratorComponent as well as any information relevant to...
G4double GetChordLength() const
Accessor.
G4ThreeVector GetPositionMiddle() const
Accessor.
G4ThreeVector GetReferencePositionEnd() const
Accessor.
G4ThreeVector GetPositionEnd() const
Accessor.
G4RotationMatrix * GetReferenceRotationStart() const
Accessor.
G4ThreeVector GetPositionStart() const
Accessor.
G4String GetPlacementName() const
Accessor.
BDSTiltOffset * GetTiltOffset() const
Accessor.
void UpdateSamplerPlacementTransform(const G4Transform3D &tranfsormIn)
G4RotationMatrix * GetReferenceRotationEnd() const
Accessor.
G4String GetName() const
Accessor.
BDSSamplerInfo * GetSamplerInfo() const
Accessor.
G4int GetIndex() const
Accessor.
G4double GetSPositionEnd() const
Accessor.
G4ThreeVector GetReferencePositionStart() const
Accessor.
BDSAcceleratorComponent * GetAcceleratorComponent() const
Accessor.
G4RotationMatrix * GetRotationMiddle() const
Accessor.
G4String GetType() const
Accessor.
BDSExtentGlobal GetExtentGlobal() const
Create a global extent object from the extent of the component.
G4double GetSPositionStart() const
Accessor.
A vector of BDSBeamlineElement instances - a beamline.
G4double totalChordLength
Sum of all chord lengths.
std::vector< G4double > GetEdgeSPositions() const
G4ThreeVector previousReferencePositionEnd
Current reference position at the end of the previous element.
const BDSBeamlineElement * GetPrevious(const BDSBeamlineElement *element) const
void PrintAllComponents(std::ostream &out) const
const BDSBeamlineElement * GetElement(G4String acceleratorComponentName, G4int i=0) const
Get the ith placement of an element in the beam line. Returns null pointer if not found.
BDSBeamlineElement * back() const
Return a reference to the last element.
G4ThreeVector maximumExtentPositive
maximum extent in the positive coordinates in each dimension
G4double previousSPositionEnd
Current s coordinate at the end of the previous element.
G4bool ElementAnglesSumToCircle() const
Check if the sum of the angle of all elements is two pi.
G4double totalAngle
Sum of all angles.
BDSBeamlineElement * ProvideEndPieceElementAfter(BDSSimpleComponent *endPiece, G4int index, G4bool flip=false) const
G4RotationMatrix * previousReferenceRotationEnd
Current reference rotation at the end of the previous element.
G4bool transformHasJustBeenApplied
void RegisterElement(BDSBeamlineElement *element)
void UpdateExtents(BDSBeamlineElement *element)
G4ThreeVector GetMaximumExtentAbsolute() const
Get the maximum extent absolute in each dimension.
void PrintMemoryConsumption() const
void ApplyTransform3D(BDSTransform3D *component)
static G4double paddingLength
The gap added for padding between each component.
const BDSBeamlineElement * at(int iElement) const
Return a reference to the element at i.
G4bool IndexOK(G4int index) const
Whether the supplied index will lie within the beam line vector.
BeamlineVector::size_type size() const
Get the number of elements.
void AddBeamlineElement(BDSBeamlineElement *element)
G4Transform3D GetTransformForElement(const G4String &acceleratorComponentName, G4int i=0) const
const_iterator FindFromS(G4double S) const
Returns an iterator to the beamline element at s.
G4double totalArcLength
Sum of all arc lengths.
BeamlineVector beamline
Vector of beam line elements - the data.
const BDSBeamlineElement * GetElementFromGlobalS(G4double S, G4int *indexOfFoundElement=nullptr) const
Return the element in this beam line according to a given s coordinate.
G4bool empty() const
Iterator mechanics.
void AddSingleComponent(BDSAcceleratorComponent *component, BDSTiltOffset *tiltOffset=nullptr, BDSSamplerInfo *samplerInfo=nullptr)
G4Transform3D GetGlobalEuclideanTransform(G4double s, G4double x=0, G4double y=0, G4int *indexOfFoundElement=nullptr) const
BDSExtentGlobal GetExtentGlobal() const
Get the global extents for this beamline.
iterator begin()
Iterator mechanics.
std::map< G4String, BDSBeamlineElement * > components
Map of objects by placement name stored in this beam line.
void AddComponent(BDSAcceleratorComponent *component, BDSTiltOffset *tiltOffset=nullptr, BDSSamplerInfo *samplerInfo=nullptr)
iterator end()
Iterator mechanics.
std::vector< G4int > GetIndicesOfElementsOfType(const G4String &type) const
Return vector of indices for this beam line where element of type name 'type' is found.
std::vector< G4double > sEnd
BeamlineVector::const_iterator const_iterator
Iterator mechanics.
std::vector< G4int > GetIndicesOfCollimators() const
Return indices in order of ecol, rcol, jcol and crystalcol elements.
G4ThreeVector maximumExtentNegative
maximum extent in the negative coordinates in each dimension
General exception with possible name of object and message.
Holder for +- extents in 3 dimensions with a rotation and translation.
std::vector< G4ThreeVector > AllBoundaryPointsGlobal() const
All 8 boundary points of the bounding box.
Holder for +- extents in 3 dimensions.
BDSExtent GetExtent() const
Accessor - see member for more info.
G4ThreeVector GetExtentPositive() const
Get the extent of the object in the positive direction in all dimensions.
G4ThreeVector GetExtentNegative() const
Get the extent of the object in the negative direction in all dimensions.
G4ThreeVector GetPlacementOffset() const
Accessor - see member for more info.
static BDSGlobalConstants * Instance()
Access method.
A class that hold multiple accelerator components.
static void PrintProtectedNames(std::ostream &out)
Feedback for protected names.
static G4bool InvalidSamplerName(const G4String &samplerName)
Test whether a sampler name is invalid or not.
All info required to build a sampler but not place it.
static G4double ChordLength()
Access the sampler plane length in other classes.
A BDSAcceleratorComponent wrapper for BDSGeometryComponent.
A holder for any placement offsets and rotations for a BDSAcceleratorComponent.
G4bool HasFiniteTilt() const
Inspector.
G4double GetTilt() const
Accessor.
G4bool HasFiniteOffset() const
Inspector.
G4ThreeVector GetOffset() const
More advance accessor for offset - only in x,y.
G4bool WillIntersect(const G4ThreeVector &incomingNormal, const G4ThreeVector &outgoingNormal, const G4double &zSeparation, const BDSExtent &incomingExtent, const BDSExtent &outgoingExtent)
G4bool StartsWith(const std::string &expression, const std::string &prefix)
Return true if a string "expression" starts with "prefix".
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())