BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSBeamline.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2023.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "BDSDebug.hh"
20#include "BDSAcceleratorComponent.hh"
21#include "BDSBeamline.hh"
22#include "BDSBeamlineElement.hh"
23#include "BDSException.hh"
24#include "BDSExtentGlobal.hh"
25#include "BDSGlobalConstants.hh"
26#include "BDSLine.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"
34
35#include "globals.hh" // geant4 globals / types
36#include "G4RotationMatrix.hh"
37#include "G4ThreeVector.hh"
38#include "G4Transform3D.hh"
39
40#include "CLHEP/Vector/AxisAngle.h"
41
42#include <algorithm>
43#include <iterator>
44#include <ostream>
45#include <set>
46#include <vector>
47
48G4double BDSBeamline::paddingLength = -1;
49
50BDSBeamline::BDSBeamline(const G4ThreeVector& initialGlobalPosition,
51 G4RotationMatrix* initialGlobalRotation,
52 G4double initialSIn):
53 sInitial(initialSIn),
54 sMaximum(initialSIn),
55 totalChordLength(0),
56 totalArcLength(0),
57 totalAngle(0),
58 previousReferencePositionEnd(initialGlobalPosition),
59 previousSPositionEnd(sInitial),
60 transformHasJustBeenApplied(false)
61{
62 // initialise extents
63 maximumExtentPositive = G4ThreeVector(0,0,0);
64 maximumExtentNegative = G4ThreeVector(0,0,0);
65
66 // initial rotation matrix
67 if (initialGlobalRotation) // default is null
68 {previousReferenceRotationEnd = initialGlobalRotation;}
69 else
70 {previousReferenceRotationEnd = new G4RotationMatrix();}
71
72 // gap between each element added to the beam line
73 if (paddingLength <= 0)
74 {paddingLength = 3 * BDSGlobalConstants::Instance()->LengthSafety();}
75 //paddingLength = 3*CLHEP::mm;
76
77#ifdef BDSDEBUG
78 G4cout << __METHOD_NAME__ << "with initial position and rotation" << G4endl;
79 G4cout << "Initial position: " << initialGlobalPosition << G4endl;
80 G4cout << "Initial rotation: " << *previousReferenceRotationEnd << G4endl;
81#endif
82}
83
84BDSBeamline::BDSBeamline(G4Transform3D initialTransform,
85 G4double initialSIn):
86 BDSBeamline(initialTransform.getTranslation(),
87 new G4RotationMatrix(initialTransform.getRotation()),
88 initialSIn)
89{;}
90
91BDSBeamline::~BDSBeamline()
92{
93 for (auto it : *this)
94 {delete it;}
95 // special case, if empty then previousReferenceRotationEnd is not used in the first element
96 if (size()==0)
98 // components map goes out of scope - elements are already deleted - no need to
99 // explicitly delete
100}
101
102void BDSBeamline::PrintAllComponents(std::ostream& out) const
103{
104 for (const auto& element : *this)
105 {out << element;}
106}
107
109{
110 G4cout << __METHOD_NAME__ << "container size: " << sizeof(beamline) << G4endl;
111 G4cout << __METHOD_NAME__ << "beamline element cumulative size: " << sizeof(BDSBeamlineElement) * beamline.size() << G4endl;
112 G4cout << __METHOD_NAME__ << "full usage including components: " << (sizeof(BDSBeamlineElement) + sizeof(BDSAcceleratorComponent)) * beamline.size() << G4endl;
113}
114
115std::ostream& operator<< (std::ostream& out, BDSBeamline const &bl)
116{
117 out << "Beamline with " << bl.size() << " elements" << G4endl;
118 out << "Total arc length: " << bl.totalArcLength << " mm" << G4endl;
119 out << "Total chord length: " << bl.totalChordLength << " mm" << G4endl;
120
121 return out;
122}
123
125 BDSTiltOffset* tiltOffset,
126 BDSSamplerInfo* samplerInfo)
127{
128 if (!component)
129 {
130 G4String samplerName = samplerInfo ? samplerInfo->name : "no_sampler_name_given";
131 throw BDSException(__METHOD_NAME__, "invalid accelerator component " + samplerName);
132 }
133
134 // check the sampler name is allowed in the output
135 if (samplerInfo)
136 {
137 G4String samplerName = samplerInfo->name;
138 if (BDSOutput::InvalidSamplerName(samplerName))
139 {
140 G4cerr << __METHOD_NAME__ << "invalid sampler name \"" << samplerName << "\"" << G4endl;
142 throw BDSException(__METHOD_NAME__, "");
143 }
144 }
145
146 if (BDSLine* line = dynamic_cast<BDSLine*>(component))
147 {
148 // in the case a single component has become a line, when we have a cylindrical
149 // sample we should flag that it should cover the full 'line', however with a plane
150 // one it is only attached to the last element.
151 BDSBeamlineElement* first = nullptr;
152 BDSBeamlineElement* last = nullptr;
153 G4int sizeLine = (G4int)line->size();
154 for (G4int i = 0; i < sizeLine; ++i)
155 {
156 if (i < sizeLine-1)
157 {
158 AddSingleComponent((*line)[i], tiltOffset);
159 if (i == 0)
160 {first = back();}
161 }
162 else
163 {// only attach the desired sampler to the last one in the line
164 AddSingleComponent((*line)[i], tiltOffset, samplerInfo);
165 last = back();
166 }
167 }
168 if (samplerInfo) // could be nullptr
169 {
170 if (samplerInfo->samplerType == BDSSamplerType::cylinder) // only a cylinder or plane can be attached to an element
171 {// cache the range it should cover as a cylinder
172 last->GetSamplerInfo()->startElement = first;
173 last->GetSamplerInfo()->finishElement = last;
174 // calculate the mid (i.e. mean) position and rotation
175 G4ThreeVector midRefPosition = (last->GetReferencePositionEnd() + first->GetReferencePositionStart()) / 2.0;
176 G4ThreeVector aaMidAxis;
177 G4double aaMidAngle;
178 auto aaStart = first->GetReferenceRotationStart()->axisAngle();
179 auto aaFinish = last->GetReferenceRotationEnd()->axisAngle();
180 // careful of identity rotations in AA form (axis=(0,0,1),angle=0) as our average of these would be wrong
181 if (first->GetReferenceRotationStart()->isIdentity())
182 {
183 aaMidAxis = aaFinish.axis();
184 aaMidAngle = 0.5 * aaFinish.delta();
185 }
186 else if (last->GetReferenceRotationEnd()->isIdentity())
187 {
188 aaMidAxis = aaStart.axis();
189 aaMidAngle = 0.5 * aaStart.delta();
190 }
191 else
192 {
193 aaMidAxis = (aaFinish.axis() + aaStart.axis()) / 2.0;
194 aaMidAngle = (aaFinish.delta() + aaStart.delta()) / 2.0;
195 }
196 auto aaCSampler = CLHEP::HepAxisAngle(aaMidAxis, aaMidAngle);
197 G4RotationMatrix rmCSampler = G4RotationMatrix(aaCSampler);
198 G4Transform3D trCSampler(rmCSampler, midRefPosition);
199 last->UpdateSamplerPlacementTransform(trCSampler);
200 }
201 }
202 }
203 else
204 {AddSingleComponent(component, tiltOffset, samplerInfo);}
205 // free memory - as once the rotations are calculated, this is no longer needed
206 delete tiltOffset;
207}
208
210 BDSTiltOffset* tiltOffset,
211 BDSSamplerInfo* samplerInfo)
212{
213#ifdef BDSDEBUG
214 G4cout << G4endl << __METHOD_NAME__ << "adding component to beamline and calculating coordinates" << G4endl;
215 G4cout << "component name: " << component->GetName() << G4endl;
216#endif
217 // Test if it's a BDSTransform3D instance - this is a unique component that requires
218 // rotation in all dimensions and can skip normal addition as isn't a real volume
219 // that can be placed. Apply the transform and skip the rest of this function by returning
220 // This modifies the "end" coordinates, rotation and axes of the last element in the beamline
221 if (BDSTransform3D* transform = dynamic_cast<BDSTransform3D*>(component))
222 {
223 ApplyTransform3D(transform);
225 return;
226 }
227
228 // if it's not a transform3d instance, continue as normal
229 // interrogate the item
230 G4double chordLength = component->GetChordLength();
231 G4double angle = component->GetAngle();
232 G4bool hasFiniteLength = BDS::IsFinite(chordLength);
233 G4bool hasFiniteAngle = BDS::IsFinite(angle);
234 G4bool hasFiniteTilt, hasFiniteOffset;
235 G4ThreeVector offset;
236 if (tiltOffset)
237 {
238 hasFiniteTilt = tiltOffset->HasFiniteTilt();
239 hasFiniteOffset = tiltOffset->HasFiniteOffset();
240 offset = tiltOffset->GetOffset(); // returns 3Vector
241 }
242 else
243 {
244 hasFiniteTilt = false;
245 hasFiniteOffset = false;
246 offset = G4ThreeVector();
247 }
248 G4ThreeVector eP = component->GetExtentPositive() + offset;
249 G4ThreeVector eN = component->GetExtentNegative() + offset;
250 G4ThreeVector placementOffset = component->GetPlacementOffset();
251 G4bool hasFinitePlacementOffset = BDS::IsFinite(placementOffset);
252 G4ThreeVector oFNormal = component->InputFaceNormal();
253
254#ifdef BDSDEBUG
255 G4cout << "chord length " << chordLength << " mm" << G4endl;
256 G4cout << "angle " << angle << " rad" << G4endl;
257 if (tiltOffset)
258 {G4cout << "tilt offsetX offsetY " << *tiltOffset << " rad mm mm " << G4endl;}
259 else
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;
270#endif
271
272 // Check this won't overlap with any previous geometry. This is only done for elements
273 // that aren't drifts as they should be built by the component factory to match any angles.
274 if (!empty() && (component->GetType() != "drift") && (component->GetType() != "thinmultipole"))
275 {// can only look back if there is an element - won't clash if no element; also add drifts always
276 G4bool keepGoing = true;
277 G4bool checkFaces = true;
278 G4double zSeparation = 0;
279 const BDSBeamlineElement* inspectedElement = back(); // remember we haven't added this new element yet
280 // find previous non drift output face.
281 G4ThreeVector iFNormal;
282 G4String clasherName = "Unknown";
283 while (keepGoing)
284 {
285 if (inspectedElement) // valid element
286 {// decrement could return nullptr so have to check if valid element
287 if ((inspectedElement->GetType() == "drift")||(inspectedElement->GetType() == "thinmultipole")) // leave keepGoing true
288 {
289 zSeparation += inspectedElement->GetChordLength();
290 inspectedElement = GetPrevious(inspectedElement); // decrement
291 }
292 else
293 {
294 keepGoing = false; // found a non drift - stop here
295 iFNormal = inspectedElement->GetAcceleratorComponent()->OutputFaceNormal();
296 clasherName = inspectedElement->GetAcceleratorComponent()->GetName();
297 }
298 }
299 else
300 {
301 keepGoing = false;
302 checkFaces = false; // got to the beginning with only drifts - don't check
303 }
304 }
305#ifdef BDSDEBUG
306 G4cout << "input face normal " << iFNormal << G4endl; // matches above debug formatting
307#endif
308
309 if (checkFaces)
310 {
311 // now do checks
312 BDSExtent extOF = inspectedElement->GetAcceleratorComponent()->GetExtent(); // output face
313 BDSExtent extIF = component->GetExtent(); // input face
314
315 G4bool willIntersect = BDS::WillIntersect(iFNormal, oFNormal, zSeparation, extIF, extOF);
316 if (willIntersect)
317 {
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;
321 throw BDSException(__METHOD_NAME__, "");
322 }
323 }
324 }
325
326 // Calculate the reference placement rotation
327 // rotations are done first as they're required to transform the spatial displacements.
328 // if not the first element in the beamline, copy the rotation matrix (cumulative along line)
329 // from end of last component, else use initial rotation matrix (no copy to prevent memory leak)
330 G4RotationMatrix* referenceRotationStart;
331 if (empty())
332 {referenceRotationStart = previousReferenceRotationEnd;}
333 else
334 {
337 referenceRotationStart = new G4RotationMatrix(*previousReferenceRotationEnd); // always create a new copy
338 }
339
340 G4RotationMatrix* referenceRotationMiddle = new G4RotationMatrix(*referenceRotationStart);
341 G4RotationMatrix* referenceRotationEnd = new G4RotationMatrix(*referenceRotationStart);
342
343 // if the component induces an angle in the reference trajectory, rotate the mid and end point
344 // rotation matrices appropriately
345 if (hasFiniteAngle)
346 {
347 // remember our definition of angle - +ve angle bends in -ve x direction in right
348 // handed coordinate system
349 // rotate about cumulative local y axis of beamline
350 // middle rotated by half angle in local x,z plane
351 G4ThreeVector rotationAxisOfBend = G4ThreeVector(0,1,0); // nominally about local unit Y
352 G4ThreeVector rotationAxisOfBendEnd = rotationAxisOfBend; // a copy
353 if (hasFiniteTilt)
354 {
355 G4double tilt = tiltOffset->GetTilt();
356 G4RotationMatrix rotationAxisRM = G4RotationMatrix();
357 rotationAxisRM.rotateZ(tilt);
358 rotationAxisOfBend.transform(rotationAxisRM);
359 rotationAxisOfBendEnd.transform(rotationAxisRM);
360 }
361 referenceRotationMiddle->rotate(angle*0.5, rotationAxisOfBend.transform(*previousReferenceRotationEnd));
362 // end rotated by full angle in local x,z plane
363 referenceRotationEnd->rotate(angle, rotationAxisOfBendEnd.transform(*previousReferenceRotationEnd));
364 }
365
366 G4RotationMatrix* rotationStart = new G4RotationMatrix(*referenceRotationStart);
367 G4RotationMatrix* rotationMiddle = new G4RotationMatrix(*referenceRotationMiddle);
368 G4RotationMatrix* rotationEnd = new G4RotationMatrix(*referenceRotationEnd);
369 // add the tilt to the rotation matrices (around z axis)
370 if (hasFiniteTilt)
371 {
372 G4double tilt = tiltOffset->GetTilt();
373
374 // transform a unit z vector with the rotation matrices to get the local axes
375 // of rotation to apply the tilt.
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));
382 }
383
384 // calculate the reference placement position
385 // if not the first item in the beamline, get the reference trajectory global position
386 // at the end of the previous element
387 if (!empty())
388 {
389 // if a transform has been applied, the previousReferencePositionEnd variable is already calculated
392 // leave a small gap for unambiguous geometry navigation. Transform that length
393 // to a unit z vector along the direction of the beam line before this component.
394 // increase it by sampler length if we're placing a sampler there.
395 G4ThreeVector pad = G4ThreeVector(0,0,paddingLength);
396 if (samplerInfo)
397 {
398 BDSSamplerType samplerType = samplerInfo->samplerType;
399 if (samplerType != BDSSamplerType::none)
400 {pad += G4ThreeVector(0,0,BDSSamplerPlane::ChordLength());}
401 }
402
403 // even if a transform has been applied that might induce a rotation, we introduce
404 // the padding length along the outgoing vector of the previous component to ensure
405 // the padding length is respected - hence we get the rotation from back() and not
406 // from the previousReferenceRotationEnd member variable
407 auto previousReferenceRotationEnd2 = back()->GetReferenceRotationEnd();
408 G4ThreeVector componentGap = pad.transform(*previousReferenceRotationEnd2);
409 previousReferencePositionEnd += componentGap;
410 }
411
412 G4ThreeVector referencePositionStart, referencePositionMiddle, referencePositionEnd;
413 if (hasFiniteLength)
414 {
415 referencePositionStart = previousReferencePositionEnd;
416
417 // calculate delta to mid point
418 G4ThreeVector md = G4ThreeVector(0, 0, 0.5 * chordLength);
419 md.transform(*referenceRotationMiddle);
420 referencePositionMiddle = referencePositionStart + md;
421 // remember the end position is the chord length along the half angle, not the full angle
422 // the particle achieves the full angle though by the end position.
423 G4ThreeVector delta = G4ThreeVector(0, 0, chordLength).transform(*referenceRotationMiddle);
424 referencePositionEnd = referencePositionStart + delta;
425 }
426 else
427 {
428 // element has no finite size so all positions are previous end position
429 // likely this is a transform3d or similar - but not hard coded just for transform3d
430 referencePositionStart = previousReferencePositionEnd;
431 referencePositionMiddle = previousReferencePositionEnd;
432 referencePositionEnd = previousReferencePositionEnd;
433 }
434
435 // add the placement offset
436 G4ThreeVector positionStart, positionMiddle, positionEnd;
437 if (hasFiniteOffset || hasFinitePlacementOffset)
438 {
439 if (hasFiniteOffset && hasFiniteAngle)
440 {// do not allow x offsets for bends as this will cause overlaps
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);
445 offset.setX(0.0);
446 }
447 // note the displacement is applied in the accelerator x and y frame so use
448 // the reference rotation rather than the one with tilt already applied
449 G4ThreeVector total = offset + placementOffset;
450 G4ThreeVector displacement = total.transform(*referenceRotationMiddle);
451 positionStart = referencePositionStart + displacement;
452 positionMiddle = referencePositionMiddle + displacement;
453 positionEnd = referencePositionEnd + displacement;
454 }
455 else
456 {
457 positionStart = referencePositionStart;
458 positionMiddle = referencePositionMiddle;
459 positionEnd = referencePositionEnd;
460 }
461
462 // calculate the s position
463 // if not the first element in the beamline, get the s position at the end of the previous element
464 if (!empty())
466
467 // chord length set earlier
468 G4double arcLength = component->GetArcLength();
469
470 // integrate lengths
471 totalChordLength += chordLength;
472 totalArcLength += arcLength;
473
474 // advance s coordinate
475 G4double sPositionStart = previousSPositionEnd;
476 G4double sPositionMiddle = previousSPositionEnd + 0.5 * arcLength;
477 G4double sPositionEnd = previousSPositionEnd + arcLength;
478 sMaximum += arcLength;
479
480 // integrate angle
481 totalAngle += component->GetAngle();
482
483#ifdef BDSDEBUG
484 // feedback about calculated coordinates
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;
498#endif
499
500 // construct beamline element
501 BDSTiltOffset* tiltOffsetToStore = nullptr;
502 if (tiltOffset)
503 {tiltOffsetToStore = new BDSTiltOffset(*tiltOffset);} // copy as can be used multiple times
504
505 BDSBeamlineElement* element;
506 element = new BDSBeamlineElement(component,
507 positionStart,
508 positionMiddle,
509 positionEnd,
510 rotationStart,
511 rotationMiddle,
512 rotationEnd,
513 referencePositionStart,
514 referencePositionMiddle,
515 referencePositionEnd,
516 referenceRotationStart,
517 referenceRotationMiddle,
518 referenceRotationEnd,
519 sPositionStart,
520 sPositionMiddle,
521 sPositionEnd,
522 tiltOffsetToStore,
523 samplerInfo,
524 (G4int)beamline.size());
525
526 // calculate extents for world size determination
527 UpdateExtents(element);
528
529 // append it to the beam line
530 beamline.push_back(element);
531
532 // register the s position at the end for curvilinear transform
533 sEnd.push_back(sPositionEnd);
534
535 // register it by name
536 RegisterElement(element);
537
538 // reset flag for transform since we've now added a component
540}
541
543{
544 G4double dx = component->dx;
545 G4double dy = component->dy;
546 G4double dz = component->dz;
547
548 // test validity for potential overlaps
549 if (dz < 0)
550 {
551 G4cerr << __METHOD_NAME__ << "Caution: Transform3d: " << component->GetName() << G4endl;
552 G4cerr << __METHOD_NAME__ << "dz = " << dz << " < 0 -> could overlap previous element" << G4endl;
553 }
554
555 // if not the first element in the beamline, get information from
556 // the end of the last element in the beamline
557 if (!empty())
558 {
559 BDSBeamlineElement* last = back();
562 }
563
564 // apply position
565 // transform the local dx,dy,dz displacement into the global frame then apply
566 G4ThreeVector delta = G4ThreeVector(dx, dy, dz).transform(*previousReferenceRotationEnd);
568
569 // apply rotation
570 G4RotationMatrix trRotInverse = component->rotationMatrix.inverse();
571 (*previousReferenceRotationEnd) *= trRotInverse;
572}
573
575{
576 if (!element)
577 {throw BDSException(__METHOD_NAME__, "invalid BDSBeamlineElement");}
578 if (!(element->GetAcceleratorComponent()))
579 {throw BDSException(__METHOD_NAME__, "invalid BDSAcceleratorComponent");}
580
581 // update world extent for this beam line
582 UpdateExtents(element);
583
584 // append it to the beam line
585 beamline.push_back(element);
586
587 // register it by name
588 RegisterElement(element);
589
590 // no need to update any internal variables - that's done by AddSingleComponent()
591}
592
594{
595 G4ThreeVector mEA;
596 for (int i=0; i<3; i++)
597 {mEA[i] = std::max(std::abs(maximumExtentPositive[i]), std::abs(maximumExtentNegative[i]));}
598 return mEA;
599}
600
601G4Transform3D BDSBeamline::GetGlobalEuclideanTransform(G4double s, G4double x, G4double y,
602 G4int* indexOfFoundElement) const
603{
604 // check if s is in the range of the beamline
605 G4double sStart = at(0)->GetSPositionStart();
606 if (s-sStart > totalArcLength) // need to offset start S position
607 {
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();
612 }
613
614 const auto element = GetElementFromGlobalS(s, indexOfFoundElement);
615
616#ifdef BDSDEBUG
617 G4cout << __METHOD_NAME__ << G4endl;
618 G4cout << "S position requested: " << s << G4endl;
619 G4cout << "Index: " << indexOfFoundElement << G4endl;
620 G4cout << "Element: " << *element << G4endl;
621#endif
622
623 G4double dx = 0;
624 // G4double dy = 0; // currently magnets can only bend in local x so avoid extra calculation
625
626 // difference from centre of element to point in local coords)
627 // difference in s from centre, normalised to arcLength and scaled to chordLength
628 // as s is really arc length, but we must place effectively in chord length coordinates
629 const BDSAcceleratorComponent* component = element->GetAcceleratorComponent();
630 G4double arcLength = component->GetArcLength();
631 G4double chordLength = component->GetChordLength();
632 G4double dS = s - element->GetSPositionMiddle();
633 G4double localZ = dS * (chordLength / arcLength);
634 G4double angle = component->GetAngle();
635 G4RotationMatrix rotation; // will be interpolated rotation
636 G4RotationMatrix* rotMiddle = element->GetReferenceRotationMiddle();
637 // find offset of point from centre of volume - 2 methods
638 if (BDS::IsFinite(angle))
639 {
640 // finite bend angle - interpolate position and angle along arc due to change in angle
641 // local unit z at start of element
642 G4ThreeVector localUnitY = G4ThreeVector(0,1,0);
643 localUnitY.transform(*(element->GetReferenceRotationStart()));
644 // linearly interpolate angle -> angle * (s from beginning into component)/arcLength
645 G4double partialAngle = angle * std::fabs(( (0.5*arcLength + dS) / arcLength));
646 rotation = G4RotationMatrix(*(element->GetReferenceRotationStart())); // start rotation
647 rotation.rotate(partialAngle, localUnitY); // rotate it by the partial angle about local Y
648 dx = localZ*tan(partialAngle); // calculate difference of arc from chord at that point
649 }
650 else
651 {rotation = G4RotationMatrix(*rotMiddle);}
652
653 // note, magnets only bend in local x so no need to add dy as always 0
654 G4ThreeVector dLocal = G4ThreeVector(x + dx, y /*+ dy*/, localZ);
655#ifdef BDSDEBUG
656 G4cout << "Local offset from middle: " << dLocal << G4endl;
657#endif
658 // note, rotation middle is also the same as the coordinate frame of the g4 solid
659 G4ThreeVector globalPos = element->GetReferencePositionMiddle() + dLocal.transform(*rotMiddle);
660 // construct transform3d from global position and rotation matrix
661 G4Transform3D result = G4Transform3D(rotation, globalPos);
662
663#ifdef BDSDEBUG
664 G4cout << "Global offset from middle: " << dLocal << G4endl;
665 G4cout << "Resultant global position: " << globalPos << G4endl;
666#endif
667 return result;
668}
669
671 G4int* indexOfFoundElement) const
672{
673 // find element that s position belongs to
674 auto lower = std::lower_bound(sEnd.begin(), sEnd.end(), S);
675 G4int index = G4int(lower - sEnd.begin()); // subtract iterators to get index
676 if (indexOfFoundElement)
677 {*indexOfFoundElement = index;}
678 return beamline.at(index);
679}
680
682{
683 auto lower = std::lower_bound(sEnd.begin(), sEnd.end(), S);
684 auto iter = begin();
685 std::advance(iter, std::distance(sEnd.begin(), lower));
686 return iter;
687}
688
690{
691 // search for element
692 auto result = find(beamline.begin(), beamline.end(), element);
693 if (result != beamline.end())
694 {// found
695 return GetPrevious(G4int(result - beamline.begin()));
696 }
697 else
698 {return nullptr;}
699}
700
701const BDSBeamlineElement* BDSBeamline::GetPrevious(G4int index) const
702{
703 if (index < 1 || index > (G4int)(beamline.size()-1))
704 {return nullptr;} // invalid index - inc beginning or end
705 else
706 {return beamline[index-1];}
707}
708
709const BDSBeamlineElement* BDSBeamline::GetNext(const BDSBeamlineElement* element) const
710{
711 // search for element
712 auto result = find(beamline.begin(), beamline.end(), element);
713 if (result != beamline.end())
714 {// found
715 return GetNext(G4int(result - beamline.begin()));
716 }
717 else
718 {return nullptr;}
719}
720
721const BDSBeamlineElement* BDSBeamline::GetNext(G4int index) const
722{
723 if (index < 0 || index > (G4int)(beamline.size()-2))
724 {return nullptr;} // invalid index - inc beginning or end
725 else
726 {return beamline[index+1];}
727}
728
730{
731 // check if base name already registered (can be single component placed multiple times)
732 const auto search = components.find(element->GetName());
733 if (search == components.end())
734 {// not registered
735 components[element->GetPlacementName()] = element;
736 }
737}
738
739const BDSBeamlineElement* BDSBeamline::GetElement(G4String acceleratorComponentName,
740 G4int i) const
741{
742 // build placement name based on acc component name and ith placement
743 // matches construction in BDSBeamlineElement
744 G4String suffix = "_" + std::to_string(i);
745 G4String placementName = acceleratorComponentName + suffix;
746 const auto search = components.find(placementName);
747 if (search == components.end())
748 {
749 // Try again but search including naming convention of uniquely built
750 // components. Sometimes we modify an element or build it uniquely for
751 // that position in the lattice, so we therefore add a suffix to the
752 // name for storing in the component registry.
753 // Naming will be NAME_MOD_MODNUMBER_PLACEMENTNUMBER
754 // Why not search registry? -> should be found from this beam line
755 // 1) search with starts with NAME
756 std::vector<const BDSBeamlineElement*> candidates;
757 std::for_each(this->begin(),
758 this->end(),
759 [&acceleratorComponentName,&candidates](const BDSBeamlineElement* el)
760 {if (BDS::StartsWith(el->GetPlacementName(), acceleratorComponentName)){candidates.push_back(el);};});
761
762 if (candidates.empty())
763 {return nullptr;} // nothing found
764 else
765 {// 2) of things that start with NAME, search for ones that end in _PLACEMENTNUMBER
766 auto foundItem = std::find_if(candidates.begin(),
767 candidates.end(),
768 [&suffix](const BDSBeamlineElement* el)
769 {return BDS::EndsWith(el->GetPlacementName(), suffix);});
770 return foundItem != candidates.end() ? *foundItem : nullptr;
771 }
772 }
773 else
774 {return search->second;}
775}
776
777G4Transform3D BDSBeamline::GetTransformForElement(const G4String& acceleratorComponentName,
778 G4int i) const
779{
780 const BDSBeamlineElement* result = GetElement(acceleratorComponentName, i);
781 if (!result)
782 {
783 G4cerr << __METHOD_NAME__ << "No element named \""
784 << acceleratorComponentName << "\" found for placement number "
785 << i << G4endl;
786 G4cout << "Note, this may be because the element is a bend and split into " << G4endl;
787 G4cout << "multiple sections with unique names." << G4endl;
788 throw BDSException(__METHOD_NAME__, "");
789 }
790 else
791 {return G4Transform3D(*(result->GetRotationMiddle()), result->GetPositionMiddle());}
792}
793
795{
796 // calculate extents for world size determination
797 // get the boundary points in global coordinates.
798 BDSExtentGlobal extG = element->GetExtentGlobal();
799 const auto boundaryPoints = extG.AllBoundaryPointsGlobal();
800
801 // expand maximums based on the boundary points.
802 for (const auto& point : boundaryPoints)
803 {
804 for (int i = 0; i < 3; ++i)
805 {
806 if (point[i] > maximumExtentPositive[i])
807 {maximumExtentPositive[i] = point[i];}
808 if (point[i] < maximumExtentNegative[i])
809 {maximumExtentNegative[i] = point[i];}
810 }
811 }
812#ifdef BDSDEBUG
813 G4cout << "new global extent +ve: " << maximumExtentPositive << G4endl;
814 G4cout << "new global extent -ve: " << maximumExtentNegative << G4endl;
815#endif
816}
817
818BDSBeamlineElement* BDSBeamline::ProvideEndPieceElementBefore(BDSSimpleComponent* endPiece,
819 G4int index) const
820{
821 if (!IndexOK(index))
822 {return nullptr;}
823
824 const G4double pl = BDSGlobalConstants::Instance()->LengthSafetyLarge(); // shortcut - 'padding length'
825 G4double endPieceLength = endPiece->GetChordLength();
826 BDSBeamlineElement* element = beamline[index];
827 G4RotationMatrix* elRotStart = element->GetRotationMiddle();
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);
831 G4double elSPosStart = element->GetSPositionStart();
832 BDSTiltOffset* elTiltOffset = element->GetTiltOffset();
833 BDSTiltOffset* forEndPiece = nullptr;
834 if (elTiltOffset)
835 {forEndPiece = new BDSTiltOffset(*elTiltOffset);}
836 BDSBeamlineElement* result = new BDSBeamlineElement(endPiece,
837 positionStart,
838 positionMiddle,
839 elPosStart,
840 new G4RotationMatrix(*elRotStart),
841 new G4RotationMatrix(*elRotStart),
842 new G4RotationMatrix(*elRotStart),
843 positionStart,// for now the same - ie no tilt offset
844 positionMiddle,
845 elPosStart,
846 new G4RotationMatrix(*elRotStart),
847 new G4RotationMatrix(*elRotStart),
848 new G4RotationMatrix(*elRotStart),
849 elSPosStart - endPieceLength,
850 elSPosStart - 0.5*endPieceLength,
851 elSPosStart,
852 forEndPiece);
853 return result;
854}
855
857 G4int index,
858 G4bool flip) const
859{
860 if (!IndexOK(index))
861 {return nullptr;}
862
863 const G4double pl = paddingLength; // shortcut
864 G4double endPieceLength = endPiece->GetChordLength();
865 BDSBeamlineElement* element = beamline[index];
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);
870 if (flip)
871 {// rotate about local unit Y direction
872 G4ThreeVector localUnitY = G4ThreeVector(0,1,0).transform(*elRotEnd);
873 elRotEnd->rotate(CLHEP::pi, localUnitY);
874 }
875 G4double elSPosEnd = element->GetSPositionEnd();
876 BDSTiltOffset* elTiltOffset = element->GetTiltOffset();
877 BDSTiltOffset* forEndPiece = nullptr;
878 if (elTiltOffset)
879 {forEndPiece = new BDSTiltOffset(*elTiltOffset);}
880 BDSBeamlineElement* result = new BDSBeamlineElement(endPiece,
881 elPosEnd,
882 positionMiddle,
883 positionEnd,
884 new G4RotationMatrix(*elRotEnd),
885 new G4RotationMatrix(*elRotEnd),
886 new G4RotationMatrix(*elRotEnd),
887 elPosEnd,
888 positionMiddle,
889 positionEnd,
890 new G4RotationMatrix(*elRotEnd),
891 new G4RotationMatrix(*elRotEnd),
892 new G4RotationMatrix(*elRotEnd),
893 elSPosEnd,
894 elSPosEnd + 0.5*endPieceLength,
895 elSPosEnd + endPieceLength,
896 forEndPiece);
897 delete elRotEnd;
898 return result;
899}
900
901G4bool BDSBeamline::IndexOK(G4int index) const
902{
903 if (index < 0 || index > (G4int)(beamline.size()-1))
904 {return false;}
905 else
906 {return true;}
907}
908
909std::vector<G4double> BDSBeamline::GetEdgeSPositions()const
910{
911 std::vector<G4double> sPos;
912 sPos.reserve(beamline.size()+1);
913 // add start position
914 sPos.push_back(0.0);
915 for (auto element : beamline)
916 {sPos.push_back(element->GetSPositionEnd()/CLHEP::m);}
917 if (sPos.size() == 1)
918 {sPos.push_back(1*CLHEP::m);}
919 return sPos;
920}
921
923{
924 return (std::abs(totalAngle) > 0.99 * 2.0 * CLHEP::pi) and (std::abs(totalAngle) < 1.01 * 2.0 * CLHEP::pi);
925}
926
928{
930 BDSExtentGlobal extG = BDSExtentGlobal(ext, G4Transform3D());
931 return extG;
932}
933
934std::vector<G4int> BDSBeamline::GetIndicesOfElementsOfType(const G4String& type) const
935{
936 std::vector<G4int> result;
937 for (auto element : beamline)
938 {
939 if (element->GetType() == type)
940 {result.push_back(element->GetIndex());}
941 }
942 return result;
943}
944
945std::vector<G4int> BDSBeamline::GetIndicesOfElementsOfType(const std::set<G4String>& types) const
946{
947 std::vector<G4int> result;
948 for (auto element : beamline)
949 {
950 G4String type = element->GetType();
951 if (types.find(type) != types.end())
952 {result.push_back(element->GetIndex());}
953 }
954 return result;
955}
956
957std::vector<G4int> BDSBeamline::GetIndicesOfCollimators() const
958{
959 std::set<G4String> collimatorTypes = {"ecol", "rcol", "jcol", "crystalcol", "element-collimator"};
960 return GetIndicesOfElementsOfType(collimatorTypes);
961}
Abstract class that represents a component of an accelerator.
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.
Definition: BDSBeamline.hh:61
G4double totalChordLength
Sum of all chord lengths.
Definition: BDSBeamline.hh:263
std::vector< G4double > GetEdgeSPositions() const
Definition: BDSBeamline.cc:909
G4ThreeVector previousReferencePositionEnd
Current reference position at the end of the previous element.
Definition: BDSBeamline.hh:276
const BDSBeamlineElement * GetPrevious(const BDSBeamlineElement *element) const
Definition: BDSBeamline.cc:689
void PrintAllComponents(std::ostream &out) const
Definition: BDSBeamline.cc:102
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.
Definition: BDSBeamline.cc:739
BDSBeamlineElement * back() const
Return a reference to the last element.
Definition: BDSBeamline.hh:214
G4ThreeVector maximumExtentPositive
maximum extent in the positive coordinates in each dimension
Definition: BDSBeamline.hh:269
G4double previousSPositionEnd
Current s coordinate at the end of the previous element.
Definition: BDSBeamline.hh:279
G4bool ElementAnglesSumToCircle() const
Check if the sum of the angle of all elements is two pi.
Definition: BDSBeamline.cc:922
G4double totalAngle
Sum of all angles.
Definition: BDSBeamline.hh:267
BDSBeamlineElement * ProvideEndPieceElementAfter(BDSSimpleComponent *endPiece, G4int index, G4bool flip=false) const
Definition: BDSBeamline.cc:856
G4RotationMatrix * previousReferenceRotationEnd
Current reference rotation at the end of the previous element.
Definition: BDSBeamline.hh:273
G4bool transformHasJustBeenApplied
Definition: BDSBeamline.hh:283
void RegisterElement(BDSBeamlineElement *element)
Definition: BDSBeamline.cc:729
void UpdateExtents(BDSBeamlineElement *element)
Definition: BDSBeamline.cc:794
G4ThreeVector GetMaximumExtentAbsolute() const
Get the maximum extent absolute in each dimension.
Definition: BDSBeamline.cc:593
void PrintMemoryConsumption() const
Definition: BDSBeamline.cc:108
void ApplyTransform3D(BDSTransform3D *component)
Definition: BDSBeamline.cc:542
static G4double paddingLength
The gap added for padding between each component.
Definition: BDSBeamline.hh:286
const BDSBeamlineElement * at(int iElement) const
Return a reference to the element at i.
Definition: BDSBeamline.hh:120
G4bool IndexOK(G4int index) const
Whether the supplied index will lie within the beam line vector.
Definition: BDSBeamline.cc:901
BeamlineVector::size_type size() const
Get the number of elements.
Definition: BDSBeamline.hh:148
void AddBeamlineElement(BDSBeamlineElement *element)
Definition: BDSBeamline.cc:574
G4Transform3D GetTransformForElement(const G4String &acceleratorComponentName, G4int i=0) const
Definition: BDSBeamline.cc:777
const_iterator FindFromS(G4double S) const
Returns an iterator to the beamline element at s.
Definition: BDSBeamline.cc:681
G4double totalArcLength
Sum of all arc lengths.
Definition: BDSBeamline.hh:265
BeamlineVector beamline
Vector of beam line elements - the data.
Definition: BDSBeamline.hh:67
const BDSBeamlineElement * GetElementFromGlobalS(G4double S, G4int *indexOfFoundElement=nullptr) const
Return the element in this beam line according to a given s coordinate.
Definition: BDSBeamline.cc:670
G4bool empty() const
Iterator mechanics.
Definition: BDSBeamline.hh:175
void AddSingleComponent(BDSAcceleratorComponent *component, BDSTiltOffset *tiltOffset=nullptr, BDSSamplerInfo *samplerInfo=nullptr)
Definition: BDSBeamline.cc:209
G4Transform3D GetGlobalEuclideanTransform(G4double s, G4double x=0, G4double y=0, G4int *indexOfFoundElement=nullptr) const
Definition: BDSBeamline.cc:601
BDSExtentGlobal GetExtentGlobal() const
Get the global extents for this beamline.
Definition: BDSBeamline.cc:927
iterator begin()
Iterator mechanics.
Definition: BDSBeamline.hh:167
std::map< G4String, BDSBeamlineElement * > components
Map of objects by placement name stored in this beam line.
Definition: BDSBeamline.hh:289
void AddComponent(BDSAcceleratorComponent *component, BDSTiltOffset *tiltOffset=nullptr, BDSSamplerInfo *samplerInfo=nullptr)
Definition: BDSBeamline.cc:124
iterator end()
Iterator mechanics.
Definition: BDSBeamline.hh:168
std::vector< G4int > GetIndicesOfElementsOfType(const G4String &type) const
Return vector of indices for this beam line where element of type name 'type' is found.
Definition: BDSBeamline.cc:934
std::vector< G4double > sEnd
Definition: BDSBeamline.hh:295
BeamlineVector::const_iterator const_iterator
Iterator mechanics.
Definition: BDSBeamline.hh:164
std::vector< G4int > GetIndicesOfCollimators() const
Return indices in order of ecol, rcol, jcol and crystalcol elements.
Definition: BDSBeamline.cc:957
G4ThreeVector maximumExtentNegative
maximum extent in the negative coordinates in each dimension
Definition: BDSBeamline.hh:270
General exception with possible name of object and message.
Definition: BDSException.hh:35
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.
Definition: BDSExtent.hh:39
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.
Definition: BDSLine.hh:38
static void PrintProtectedNames(std::ostream &out)
Feedback for protected names.
Definition: BDSOutput.cc:398
static G4bool InvalidSamplerName(const G4String &samplerName)
Test whether a sampler name is invalid or not.
Definition: BDSOutput.cc:393
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.
Transform in beam line that takes up no space.
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())