BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSBeamPipeFactoryPointsFile.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 "BDSAperturePointsLoader.hh"
20#include "BDSBeamPipeFactoryPointsFile.hh"
21
22#include "G4Types.hh"
23#include "G4TwoVector.hh"
24
25#include <algorithm>
26#include <vector>
27
28BDSBeamPipeFactoryPointsFile::BDSBeamPipeFactoryPointsFile():
29 maximumRadius(0)
30{;}
31
32void BDSBeamPipeFactoryPointsFile::GeneratePoints(G4double, G4double, G4double, G4double,
33 G4double beamPipeThickness, G4int)
34{
35 maximumRadius = 0;
36
37 // the product will not own this points vector pointer - the cache does inside this function
38 std::vector<G4TwoVector>* rawPoints = BDS::LoadAperturePoints(pointsFile, pointsUnit);
39
40 // create a local copy
41 vacuumEdge = std::vector<G4TwoVector>(rawPoints->begin(), rawPoints->end());
42
43 // determine the winding order and if wrong, flip it
44 // we rely on this order, so we calculate the outward facing normals correctly
45 // further down - the sign of the pi/2
46 auto detV1 = vacuumEdge[1] - vacuumEdge[0];
47 auto detV2 = vacuumEdge[2] - vacuumEdge[1];
48 G4double det = Determinant(detV1, detV2);
49
50 if (det < 0)
51 {std::reverse(vacuumEdge.begin(), vacuumEdge.end());}
52
53 // calculate the unit vector for each point - so a direction from the 0,0 point to that point
54 std::vector<G4TwoVector> vertexNormals;
55 vertexNormals.reserve(vacuumEdge.size());
56
57 // Loop round the polygon and calculate the normal of each segment
58 // sum the normals of two segments and normalise to give an outwards
59 // pointing unit vector at each vertex. We have a bit of if else for
60 // the boundaries of the loop for back / front.
61 G4TwoVector v1;
62 G4TwoVector v2;
63 G4TwoVector v3;
64 for (G4int i = 0; i < (G4int)vacuumEdge.size(); i++)
65 {
66 if (i == 0)
67 {
68 v1 = vacuumEdge.back();
69 v3 = vacuumEdge[i+1];
70 }
71 else if (i == (G4int)vacuumEdge.size() - 1)
72 {
73 v1 = vacuumEdge[i-1];
74 v3 = vacuumEdge[0];
75 }
76 else
77 {
78 v1 = vacuumEdge[i-1];
79 v3 = vacuumEdge[i+1];
80 }
81 v2 = vacuumEdge[i];
82
83 auto normal1 = v2 - v1;
84 auto normal2 = v3 - v2;
85 normal1.rotate(-CLHEP::halfpi); // now it's a normal
86 normal2.rotate(-CLHEP::halfpi);
87
88 auto vertexNormal = (normal1 + normal2).unit();
89
90 vertexNormals.push_back(vertexNormal); // outward pointing direction for expansion
91 }
92
93 for (G4int i = 0; i < (G4int)vacuumEdge.size(); i++)
94 {
95 auto ur = vertexNormals[i];
96 beamPipeInnerEdge.emplace_back( vacuumEdge[i] + ur*lengthSafetyLarge);
97 beamPipeOuterEdge.emplace_back( vacuumEdge[i] + ur*(lengthSafetyLarge + beamPipeThickness));
98 containerEdge.emplace_back( vacuumEdge[i] + ur*(5*lengthSafetyLarge + beamPipeThickness));
99 containerSubtractionEdge.emplace_back(vacuumEdge[i] + ur*(7*lengthSafetyLarge + beamPipeThickness));
100 maximumRadius = std::max(maximumRadius, containerSubtractionEdge.back().mag());
101 }
102}
103
104G4double BDSBeamPipeFactoryPointsFile::CalculateIntersectionRadius(G4double, G4double, G4double, G4double, G4double)
105{
106 return maximumRadius*1.01; // 1% margin on absolute size
107}
108
109G4double BDSBeamPipeFactoryPointsFile::Determinant(const G4TwoVector& v1,
110 const G4TwoVector& v2) const
111{
112 G4double det = v1.x()*v2.y() - v2.x()*v1.y();
113 return det;
114}
G4double Determinant(const G4TwoVector &v1, const G4TwoVector &v2) const
Calculate the determinant of 2x G4TwoVectors.
virtual G4double CalculateIntersectionRadius(G4double aper1, G4double aper2, G4double aper3, G4double aper4, G4double beamPipeThickness)
Calculate the radius of the solid used for intersection for angled faces.
virtual void GeneratePoints(G4double, G4double, G4double, G4double, G4double beamPipeThickness, G4int=40)
Purely to fulfill interface - should not be used!
std::vector< G4TwoVector > containerSubtractionEdge
Vector of x,y coordinates for the container subtraction volume.
std::vector< G4TwoVector > beamPipeInnerEdge
Vector of x,y coordinates for beam pipe inner edge.
std::vector< G4TwoVector > vacuumEdge
Vector of x,y coordinates for vacuum extruded solid edge.
G4String pointsFile
Temporary copy of input variables.
std::vector< G4TwoVector > containerEdge
Vector of x,y coordinates for the container volume.
std::vector< G4TwoVector > beamPipeOuterEdge
Vector of x,y coordinates for beam pipe outer edge.
G4String pointsUnit
Temporary copy of input variables.
G4double lengthSafetyLarge
Cache of global constants variable.
std::vector< G4TwoVector > * LoadAperturePoints(const G4String &fileName, const G4String &unit="")