| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4VTwistSurface implementation <<
27 // 26 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu << 27 // $Id: G4VTwistSurface.cc,v 1.5 2006/06/29 18:49:34 gunter Exp $
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), << 28 // GEANT4 tag $Name: geant4-08-01-patch-01 $
30 // from original version in Jupi << 29 //
>> 30 //
>> 31 // --------------------------------------------------------------------
>> 32 // GEANT 4 class source file
>> 33 //
>> 34 //
>> 35 // G4VTwistSurface.cc
>> 36 //
>> 37 // Author:
>> 38 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
>> 39 //
>> 40 // History:
>> 41 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
>> 42 // from original version in Jupiter-2.5.02 application.
31 // ------------------------------------------- 43 // --------------------------------------------------------------------
32 <<
33 #include <iomanip> <<
34 44
35 #include "G4VTwistSurface.hh" 45 #include "G4VTwistSurface.hh"
36 #include "G4GeometryTolerance.hh" << 46 #include <iomanip>
37 47
38 const G4int G4VTwistSurface::sOutside 48 const G4int G4VTwistSurface::sOutside = 0x00000000;
39 const G4int G4VTwistSurface::sInside 49 const G4int G4VTwistSurface::sInside = 0x10000000;
40 const G4int G4VTwistSurface::sBoundary 50 const G4int G4VTwistSurface::sBoundary = 0x20000000;
41 const G4int G4VTwistSurface::sCorner 51 const G4int G4VTwistSurface::sCorner = 0x40000000;
42 const G4int G4VTwistSurface::sC0Min1Min 52 const G4int G4VTwistSurface::sC0Min1Min = 0x40000101;
43 const G4int G4VTwistSurface::sC0Max1Min 53 const G4int G4VTwistSurface::sC0Max1Min = 0x40000201;
44 const G4int G4VTwistSurface::sC0Max1Max 54 const G4int G4VTwistSurface::sC0Max1Max = 0x40000202;
45 const G4int G4VTwistSurface::sC0Min1Max 55 const G4int G4VTwistSurface::sC0Min1Max = 0x40000102;
46 const G4int G4VTwistSurface::sAxisMin 56 const G4int G4VTwistSurface::sAxisMin = 0x00000101;
47 const G4int G4VTwistSurface::sAxisMax 57 const G4int G4VTwistSurface::sAxisMax = 0x00000202;
48 const G4int G4VTwistSurface::sAxisX 58 const G4int G4VTwistSurface::sAxisX = 0x00000404;
49 const G4int G4VTwistSurface::sAxisY 59 const G4int G4VTwistSurface::sAxisY = 0x00000808;
50 const G4int G4VTwistSurface::sAxisZ 60 const G4int G4VTwistSurface::sAxisZ = 0x00000C0C;
51 const G4int G4VTwistSurface::sAxisRho 61 const G4int G4VTwistSurface::sAxisRho = 0x00001010;
52 const G4int G4VTwistSurface::sAxisPhi 62 const G4int G4VTwistSurface::sAxisPhi = 0x00001414;
53 63
54 // mask 64 // mask
55 const G4int G4VTwistSurface::sAxis0 65 const G4int G4VTwistSurface::sAxis0 = 0x0000FF00;
56 const G4int G4VTwistSurface::sAxis1 66 const G4int G4VTwistSurface::sAxis1 = 0x000000FF;
57 const G4int G4VTwistSurface::sSizeMask 67 const G4int G4VTwistSurface::sSizeMask = 0x00000303;
58 const G4int G4VTwistSurface::sAxisMask 68 const G4int G4VTwistSurface::sAxisMask = 0x0000FCFC;
59 const G4int G4VTwistSurface::sAreaMask 69 const G4int G4VTwistSurface::sAreaMask = 0XF0000000;
60 70
61 //============================================ 71 //=====================================================================
62 //* constructors ----------------------------- 72 //* constructors ------------------------------------------------------
63 73
64 G4VTwistSurface::G4VTwistSurface(const G4Strin 74 G4VTwistSurface::G4VTwistSurface(const G4String &name)
65 : fIsValidNorm(false), fName(name) << 75 : fName(name)
66 { 76 {
67 77
68 fAxis[0] = kUndefined; 78 fAxis[0] = kUndefined;
69 fAxis[1] = kUndefined; 79 fAxis[1] = kUndefined;
70 fAxisMin[0] = kInfinity; 80 fAxisMin[0] = kInfinity;
71 fAxisMin[1] = kInfinity; 81 fAxisMin[1] = kInfinity;
72 fAxisMax[0] = kInfinity; 82 fAxisMax[0] = kInfinity;
73 fAxisMax[1] = kInfinity; 83 fAxisMax[1] = kInfinity;
74 fHandedness = 1; 84 fHandedness = 1;
75 85
76 for (auto i=0; i<4; ++i) << 86 G4int i;
77 { << 87 for (i=0; i<4; i++) {
78 fCorners[i].set(kInfinity, kInfinity, kI 88 fCorners[i].set(kInfinity, kInfinity, kInfinity);
79 fNeighbours[i] = nullptr; << 89 fNeighbours[i] = 0;
80 } 90 }
81 91
82 fCurrentNormal.p.set(kInfinity, kInfinity, 92 fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity);
83 93
84 fAmIOnLeftSide.me.set(kInfinity, kInfinity, 94 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
85 fAmIOnLeftSide.vec.set(kInfinity, kInfinity 95 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
86 kCarTolerance = G4GeometryTolerance::GetIns <<
87 } 96 }
88 97
89 G4VTwistSurface::G4VTwistSurface(const G4Strin << 98 G4VTwistSurface::G4VTwistSurface(const G4String &name,
90 const G4RotationMatrix& << 99 const G4RotationMatrix &rot,
91 const G4ThreeVector& << 100 const G4ThreeVector &tlate,
92 G4int 101 G4int handedness,
93 const EAxis 102 const EAxis axis0 ,
94 const EAxis 103 const EAxis axis1 ,
95 G4double 104 G4double axis0min,
96 G4double 105 G4double axis1min,
97 G4double 106 G4double axis0max,
98 G4double 107 G4double axis1max )
99 : fIsValidNorm(false), fName(name) << 108 : fName(name)
100 { 109 {
101 fAxis[0] = axis0; 110 fAxis[0] = axis0;
102 fAxis[1] = axis1; 111 fAxis[1] = axis1;
103 fAxisMin[0] = axis0min; 112 fAxisMin[0] = axis0min;
104 fAxisMin[1] = axis1min; 113 fAxisMin[1] = axis1min;
105 fAxisMax[0] = axis0max; 114 fAxisMax[0] = axis0max;
106 fAxisMax[1] = axis1max; 115 fAxisMax[1] = axis1max;
107 fHandedness = handedness; 116 fHandedness = handedness;
108 fRot = rot; 117 fRot = rot;
109 fTrans = tlate; 118 fTrans = tlate;
110 119
111 for (auto i=0; i<4; ++i) << 120 G4int i;
112 { << 121 for (i=0; i<4; i++) {
113 fCorners[i].set(kInfinity, kInfinity, kI 122 fCorners[i].set(kInfinity, kInfinity, kInfinity);
114 fNeighbours[i] = nullptr; << 123 fNeighbours[i] = 0;
115 } 124 }
116 125
117 fCurrentNormal.p.set(kInfinity, kInfinity, 126 fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity);
118 127
119 fAmIOnLeftSide.me.set(kInfinity, kInfinity, 128 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
120 fAmIOnLeftSide.vec.set(kInfinity, kInfinity 129 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
121 kCarTolerance = G4GeometryTolerance::GetIns <<
122 } 130 }
123 131
124 //============================================ 132 //=====================================================================
125 //* Fake default constructor ----------------- 133 //* Fake default constructor ------------------------------------------
126 134
127 G4VTwistSurface::G4VTwistSurface( __void__& ) 135 G4VTwistSurface::G4VTwistSurface( __void__& )
128 : fHandedness(0), fIsValidNorm(false), kCarT << 136 : fName("")
129 fName("") << 137 {
>> 138 }
>> 139
>> 140 //=====================================================================
>> 141 //* destructor --------------------------------------------------------
>> 142
>> 143 G4VTwistSurface::~G4VTwistSurface()
130 { 144 {
131 fAxis[0] = fAxis[1] = kXAxis; <<
132 fAxisMin[0] = fAxisMin[1] = 0.; <<
133 fAxisMax[0] = fAxisMax[1] = 0.; <<
134 fNeighbours[0] = fNeighbours[1] = fNeighbou <<
135 } 145 }
136 146
137 //============================================ 147 //=====================================================================
138 //* AmIOnLeftSide ---------------------------- 148 //* AmIOnLeftSide -----------------------------------------------------
139 149
140 G4int G4VTwistSurface::AmIOnLeftSide(const G4T << 150 G4int G4VTwistSurface::AmIOnLeftSide(const G4ThreeVector &me,
141 const G4T << 151 const G4ThreeVector &vec,
142 G4b << 152 G4bool withtol)
143 { 153 {
144 // AmIOnLeftSide returns phi-location of "m 154 // AmIOnLeftSide returns phi-location of "me"
145 // (phi relation between me and vec project 155 // (phi relation between me and vec projected on z=0 plane).
146 // If "me" is on -ve-phi-side of "vec", it 156 // If "me" is on -ve-phi-side of "vec", it returns 1.
147 // On the other hand, if "me" is on +ve-phi 157 // On the other hand, if "me" is on +ve-phi-side of "vec",
148 // it returns -1. 158 // it returns -1.
149 // (The return value represents z-coordinat 159 // (The return value represents z-coordinate of normal vector
150 // of me.cross(vec).) 160 // of me.cross(vec).)
151 // If me is on boundary of vec, return 0. 161 // If me is on boundary of vec, return 0.
152 162
153 const G4double kAngTolerance << 163 static G4RotationMatrix unitrot; // unit matrix
154 = G4GeometryTolerance::GetInstance()->Get << 164 static const G4RotationMatrix rottol = unitrot.rotateZ(0.5*kAngTolerance);
155 << 165 static const G4RotationMatrix invrottol = unitrot.rotateZ(-1.*kAngTolerance);
156 G4RotationMatrix unitrot; <<
157 const G4RotationMatrix rottol = unitrot. <<
158 const G4RotationMatrix invrottol = unitrot. <<
159 166
160 if (fAmIOnLeftSide.me == me 167 if (fAmIOnLeftSide.me == me
161 && fAmIOnLeftSide.vec == vec 168 && fAmIOnLeftSide.vec == vec
162 && fAmIOnLeftSide.withTol == withtol) << 169 && fAmIOnLeftSide.withTol == withtol) {
163 { <<
164 return fAmIOnLeftSide.amIOnLeftSide; 170 return fAmIOnLeftSide.amIOnLeftSide;
165 } 171 }
166 172
167 fAmIOnLeftSide.me = me; 173 fAmIOnLeftSide.me = me;
168 fAmIOnLeftSide.vec = vec; 174 fAmIOnLeftSide.vec = vec;
169 fAmIOnLeftSide.withTol = withtol; 175 fAmIOnLeftSide.withTol = withtol;
170 176
171 G4ThreeVector met = (G4ThreeVector(me.x() 177 G4ThreeVector met = (G4ThreeVector(me.x(), me.y(), 0.)).unit();
172 G4ThreeVector vect = (G4ThreeVector(vec.x( 178 G4ThreeVector vect = (G4ThreeVector(vec.x(), vec.y(), 0.)).unit();
173 179
174 G4ThreeVector ivect = invrottol * vect; 180 G4ThreeVector ivect = invrottol * vect;
175 G4ThreeVector rvect = rottol * vect; 181 G4ThreeVector rvect = rottol * vect;
176 182
177 G4double metcrossvect = met.x() * vect.y() 183 G4double metcrossvect = met.x() * vect.y() - met.y() * vect.x();
178 184
179 if (withtol) << 185 if (withtol) {
180 { <<
181 if (met.x() * ivect.y() - met.y() * ivec 186 if (met.x() * ivect.y() - met.y() * ivect.x() > 0 &&
182 metcrossvect >= 0) { 187 metcrossvect >= 0) {
183 fAmIOnLeftSide.amIOnLeftSide = 1; 188 fAmIOnLeftSide.amIOnLeftSide = 1;
184 } else if (met.x() * rvect.y() - met.y() 189 } else if (met.x() * rvect.y() - met.y() * rvect.x() < 0 &&
185 metcrossvect <= 0) { 190 metcrossvect <= 0) {
186 fAmIOnLeftSide.amIOnLeftSide = -1; 191 fAmIOnLeftSide.amIOnLeftSide = -1;
187 } else { 192 } else {
188 fAmIOnLeftSide.amIOnLeftSide = 0; 193 fAmIOnLeftSide.amIOnLeftSide = 0;
189 } 194 }
190 } << 195 } else {
191 else <<
192 { <<
193 if (metcrossvect > 0) { 196 if (metcrossvect > 0) {
194 fAmIOnLeftSide.amIOnLeftSide = 1; 197 fAmIOnLeftSide.amIOnLeftSide = 1;
195 } else if (metcrossvect < 0 ) { 198 } else if (metcrossvect < 0 ) {
196 fAmIOnLeftSide.amIOnLeftSide = -1; 199 fAmIOnLeftSide.amIOnLeftSide = -1;
197 } else { 200 } else {
198 fAmIOnLeftSide.amIOnLeftSide = 0; 201 fAmIOnLeftSide.amIOnLeftSide = 0;
199 } 202 }
200 } 203 }
201 204
202 #ifdef G4TWISTDEBUG << 205 #ifdef G4SPECSDEBUG
203 G4cout << " === G4VTwistSurface::Am 206 G4cout << " === G4VTwistSurface::AmIOnLeftSide() =============="
204 << G4endl; 207 << G4endl;
205 G4cout << " Name , returncode 208 G4cout << " Name , returncode : " << fName << " "
206 << fAmIOnLeftSide.amIOn 209 << fAmIOnLeftSide.amIOnLeftSide << G4endl;
207 G4cout << " me, vec : " << s 210 G4cout << " me, vec : " << std::setprecision(14) << me
208 << " 211 << " " << vec << G4endl;
209 G4cout << " met, vect : " << m 212 G4cout << " met, vect : " << met << " " << vect << G4endl;
210 G4cout << " ivec, rvec : " << i 213 G4cout << " ivec, rvec : " << ivect << " " << rvect << G4endl;
211 G4cout << " met x vect : " << m 214 G4cout << " met x vect : " << metcrossvect << G4endl;
212 G4cout << " met x ivec : " << m 215 G4cout << " met x ivec : " << met.cross(ivect) << G4endl;
213 G4cout << " met x rvec : " << m 216 G4cout << " met x rvec : " << met.cross(rvect) << G4endl;
214 G4cout << " ======================= 217 G4cout << " =============================================="
215 << G4endl; 218 << G4endl;
216 #endif 219 #endif
217 220
218 return fAmIOnLeftSide.amIOnLeftSide; 221 return fAmIOnLeftSide.amIOnLeftSide;
219 } 222 }
220 223
221 //============================================ 224 //=====================================================================
222 //* DistanceToBoundary ----------------------- 225 //* DistanceToBoundary ------------------------------------------------
223 226
224 G4double G4VTwistSurface::DistanceToBoundary(G 227 G4double G4VTwistSurface::DistanceToBoundary(G4int areacode,
225 G << 228 G4ThreeVector &xx,
226 const G << 229 const G4ThreeVector &p)
227 { 230 {
228 // DistanceToBoundary 231 // DistanceToBoundary
229 // 232 //
230 // return distance to nearest boundary from 233 // return distance to nearest boundary from arbitrary point p
231 // in local coodinate. 234 // in local coodinate.
232 // Argument areacode must be one of them: 235 // Argument areacode must be one of them:
233 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 236 // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
234 // sAxis1 & sAxisMin, sAxis1 & sAxisMax. 237 // sAxis1 & sAxisMin, sAxis1 & sAxisMax.
235 // 238 //
236 239
237 G4ThreeVector d; // direction vector of 240 G4ThreeVector d; // direction vector of the boundary
238 G4ThreeVector x0; // reference point of t 241 G4ThreeVector x0; // reference point of the boundary
239 G4double dist = kInfinity; 242 G4double dist = kInfinity;
240 G4int boundarytype; 243 G4int boundarytype;
241 244
242 if (IsAxis0(areacode) && IsAxis1(areacode)) << 245 if (IsAxis0(areacode) && IsAxis1(areacode)) {
243 { << 246 G4cerr << "ERROR - G4VTwistSurface::DistanceToBoundary()" << G4endl
244 std::ostringstream message; << 247 << " Point is in the corner area. This function returns"
245 message << "Point is in the corner area. << 248 << G4endl
246 << " Point is in the corn << 249 << " a direction vector of a boundary line." << G4endl
247 << G4endl << 250 << " areacode = " << areacode << G4endl;
248 << " a direction vector o << 251 G4Exception("G4VTwistSurface::DistanceToBoundary()", "InvalidSetup",
249 << " areacode = " << area << 252 FatalException, "Point is in the corner area.");
250 G4Exception("G4VTwistSurface::DistanceTo << 253 } else if (IsAxis0(areacode) || IsAxis1(areacode)) {
251 FatalException, message); <<
252 } <<
253 else if (IsAxis0(areacode) || IsAxis1(areac <<
254 { <<
255 GetBoundaryParameters(areacode, d, x0, b 254 GetBoundaryParameters(areacode, d, x0, boundarytype);
256 if (boundarytype == sAxisPhi) << 255 if (boundarytype == sAxisPhi) {
257 { <<
258 G4double t = x0.getRho() / p.getRho() 256 G4double t = x0.getRho() / p.getRho();
259 xx.set(t*p.x(), t*p.y(), x0.z()); 257 xx.set(t*p.x(), t*p.y(), x0.z());
260 dist = (xx - p).mag(); 258 dist = (xx - p).mag();
261 } << 259 } else {
262 else <<
263 { <<
264 // linear boundary 260 // linear boundary
265 // sAxisX, sAxisY, sAxisZ, sAxisRho 261 // sAxisX, sAxisY, sAxisZ, sAxisRho
266 dist = DistanceToLine(p, x0, d, xx); 262 dist = DistanceToLine(p, x0, d, xx);
267 } 263 }
268 } << 264 } else {
269 else << 265 G4cerr << "ERROR - G4VTwistSurface::DistanceToBoundary()" << G4endl
270 { << 266 << " areacode = " << areacode << G4endl;
271 std::ostringstream message; << 267 G4Exception("G4VTwistSurface::DistanceToBoundary()", "InvalidSetup",
272 message << "Bad areacode of boundary." < << 268 FatalException, "Bad areacode of boundary.");
273 << " areacode = " << area <<
274 G4Exception("G4VTwistSurface::DistanceTo <<
275 FatalException, message); <<
276 } 269 }
277 return dist; 270 return dist;
278 } 271 }
279 272
280 //============================================ 273 //=====================================================================
281 //* DistanceToIn ----------------------------- 274 //* DistanceToIn ------------------------------------------------------
282 275
283 G4double G4VTwistSurface::DistanceToIn(const G << 276 G4double G4VTwistSurface::DistanceToIn(const G4ThreeVector &gp,
284 const G << 277 const G4ThreeVector &gv,
285 G << 278 G4ThreeVector &gxxbest)
286 { 279 {
287 #ifdef G4TWISTDEBUG << 280 #ifdef G4SPECSDEBUG
288 G4cout << " ~~~~ G4VTwistSurface::DistanceT << 281 G4cout << " ~~~~~ G4VTwistSurface::DistanceToIn(p,v) - Start ~~~~~" << G4endl;
289 G4cout << " Name : " << fName << G4end 282 G4cout << " Name : " << fName << G4endl;
290 G4cout << " gp : " << gp << G4endl; 283 G4cout << " gp : " << gp << G4endl;
291 G4cout << " gv : " << gv << G4endl; 284 G4cout << " gv : " << gv << G4endl;
292 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 285 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
293 #endif 286 #endif
294 287
295 G4ThreeVector gxx[G4VSURFACENXX]; 288 G4ThreeVector gxx[G4VSURFACENXX];
296 G4double distance[G4VSURFACENXX] ; 289 G4double distance[G4VSURFACENXX] ;
297 G4int areacode[G4VSURFACENXX] ; 290 G4int areacode[G4VSURFACENXX] ;
298 G4bool isvalid[G4VSURFACENXX] ; 291 G4bool isvalid[G4VSURFACENXX] ;
299 292
300 for (G4int i = 0 ; i<G4VSURFACENXX ; ++i ) << 293 for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) {
301 { <<
302 distance[i] = kInfinity ; 294 distance[i] = kInfinity ;
303 areacode[i] = sOutside ; 295 areacode[i] = sOutside ;
304 isvalid[i] = false ; 296 isvalid[i] = false ;
305 } 297 }
306 298
307 G4double bestdistance = kInfinity; 299 G4double bestdistance = kInfinity;
308 #ifdef G4TWISTDEBUG <<
309 G4int besti = -1; 300 G4int besti = -1;
310 #endif <<
311 G4ThreeVector bestgxx(kInfinity, kInfinity, 301 G4ThreeVector bestgxx(kInfinity, kInfinity, kInfinity);
312 302
313 G4int nxx = DistanceToSurface(gp, g 303 G4int nxx = DistanceToSurface(gp, gv, gxx, distance, areacode,
314 isval 304 isvalid, kValidateWithTol);
315 305
316 for (G4int i=0; i<nxx; ++i) << 306 for (G4int i=0; i< nxx; i++) {
317 { <<
318 307
319 // skip this intersection if: 308 // skip this intersection if:
320 // - invalid intersection 309 // - invalid intersection
321 // - particle goes outword the surface 310 // - particle goes outword the surface
322 311
323 if (!isvalid[i]) << 312 if (!isvalid[i]) {
324 { <<
325 // xx[i] is sOutside or distance[i] < 313 // xx[i] is sOutside or distance[i] < 0
326 continue; 314 continue;
327 } 315 }
328 316
329 G4ThreeVector normal = GetNormal(gxx[i], 317 G4ThreeVector normal = GetNormal(gxx[i], true);
330 318
331 if ((normal * gv) >= 0) << 319 if ((normal * gv) >= 0) {
332 { <<
333 320
334 #ifdef G4TWISTDEBUG << 321 #ifdef G4SPECSDEBUG
335 G4cout << " G4VTwistSurface::Distan 322 G4cout << " G4VTwistSurface::DistanceToIn(p,v): "
336 << "particle goes outword the 323 << "particle goes outword the surface." << G4endl;
337 #endif 324 #endif
338 continue; 325 continue;
339 } 326 }
340 327
341 // 328 //
342 // accept this intersection if the inter 329 // accept this intersection if the intersection is inside.
343 // 330 //
344 331
345 if (IsInside(areacode[i])) << 332 if (IsInside(areacode[i])) {
346 { << 333 if (distance[i] < bestdistance) {
347 if (distance[i] < bestdistance) <<
348 { <<
349 bestdistance = distance[i]; 334 bestdistance = distance[i];
350 bestgxx = gxx[i]; 335 bestgxx = gxx[i];
351 #ifdef G4TWISTDEBUG <<
352 besti = i; 336 besti = i;
>> 337
>> 338 #ifdef G4SPECSDEBUG
353 G4cout << " G4VTwistSurface::Dis 339 G4cout << " G4VTwistSurface::DistanceToIn(p,v): "
354 << " areacode sInside name, 340 << " areacode sInside name, distance = "
355 << fName << " "<< bestdist 341 << fName << " "<< bestdistance << G4endl;
356 #endif 342 #endif
357 } 343 }
358 344
359 // 345 //
360 // else, the intersection is on boundary 346 // else, the intersection is on boundary or corner.
361 // 347 //
362 348
363 } << 349 } else {
364 else << 350
365 { << 351 G4VTwistSurface *neighbours[2];
366 G4VTwistSurface* neighbours[2]; <<
367 G4bool isaccepted[2] = {false, f 352 G4bool isaccepted[2] = {false, false};
368 G4int nneighbours = GetNeighb 353 G4int nneighbours = GetNeighbours(areacode[i], neighbours);
369 354
370 for (G4int j=0; j<nneighbours; ++j) << 355 for (G4int j=0; j< nneighbours; j++) {
371 { <<
372 // if on corner, nneighbours = 2. 356 // if on corner, nneighbours = 2.
373 // if on boundary, nneighbours = 1 357 // if on boundary, nneighbours = 1.
374 358
375 G4ThreeVector tmpgxx[G4VSURFACENXX 359 G4ThreeVector tmpgxx[G4VSURFACENXX];
376 G4double tmpdist[G4VSURFACENX 360 G4double tmpdist[G4VSURFACENXX] ;
377 G4int tmpareacode[G4VSURFA 361 G4int tmpareacode[G4VSURFACENXX] ;
378 G4bool tmpisvalid[G4VSURFAC 362 G4bool tmpisvalid[G4VSURFACENXX] ;
379 363
380 for (G4int l = 0 ; l<G4VSURFACENXX << 364 for (G4int l = 0 ; l<G4VSURFACENXX ; l++ ) {
381 { <<
382 tmpdist[l] = kInfinity ; 365 tmpdist[l] = kInfinity ;
383 tmpareacode[l] = sOutside ; 366 tmpareacode[l] = sOutside ;
384 tmpisvalid[l] = false ; 367 tmpisvalid[l] = false ;
385 } 368 }
386 369
387 G4int tmpnxx = neighbours[j]->Dist 370 G4int tmpnxx = neighbours[j]->DistanceToSurface(
388 gp, 371 gp, gv, tmpgxx, tmpdist,
389 tmpa 372 tmpareacode, tmpisvalid,
390 kVal 373 kValidateWithTol);
391 G4ThreeVector neighbournormal; 374 G4ThreeVector neighbournormal;
392 375
393 for (G4int k=0; k< tmpnxx; ++k) << 376 for (G4int k=0; k< tmpnxx; k++) {
394 { << 377
395 // 378 //
396 // if tmpxx[k] is valid && sIns 379 // if tmpxx[k] is valid && sInside, the final winner must
397 // be neighbour surface. return 380 // be neighbour surface. return kInfinity.
398 // else , choose tmpxx on same 381 // else , choose tmpxx on same boundary of xx, then check normal
399 // 382 //
400 383
401 if (IsInside(tmpareacode[k])) << 384 if (IsInside(tmpareacode[k])) {
402 { << 385
403 #ifdef G4TWISTDEBUG << 386 #ifdef G4SPECSDEBUG
404 G4cout << " G4VTwistSurfac 387 G4cout << " G4VTwistSurface:DistanceToIn(p,v): "
405 << " intersection "<< 388 << " intersection "<< tmpgxx[k] << G4endl
406 << " is inside of n 389 << " is inside of neighbour surface of " << fName
407 << " . returning kInf 390 << " . returning kInfinity." << G4endl;
408 G4cout << "~~ G4VTwistSurfac << 391 G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~"
409 << G4endl; 392 << G4endl;
410 G4cout << " No intersec 393 G4cout << " No intersections " << G4endl;
411 G4cout << " Name : " << 394 G4cout << " Name : " << fName << G4endl;
412 G4cout << "~~~~~~~~~~~~~~~~~ 395 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
413 << G4endl; 396 << G4endl;
414 #endif 397 #endif
415 if (tmpisvalid[k]) return k 398 if (tmpisvalid[k]) return kInfinity;
416 continue; 399 continue;
417 400
418 // 401 //
419 // if tmpxx[k] is valid && sIns 402 // if tmpxx[k] is valid && sInside, the final winner must
420 // be neighbour surface. return 403 // be neighbour surface. return .
421 // 404 //
422 405
423 } << 406 } else if (IsSameBoundary(this,areacode[i],
424 else if (IsSameBoundary(this,ar << 407 neighbours[j], tmpareacode[k])) {
425 neighbo <<
426 { <<
427 // tmpxx[k] is same boundary 408 // tmpxx[k] is same boundary (or corner) of xx.
428 409
429 neighbournormal = neighbours 410 neighbournormal = neighbours[j]->GetNormal(tmpgxx[k], true);
430 if (neighbournormal * gv < 0 411 if (neighbournormal * gv < 0) isaccepted[j] = true;
431 } 412 }
432 } 413 }
433 414
434 // if nneighbours = 1, chabge isac 415 // if nneighbours = 1, chabge isaccepted[1] before
435 // exiting neighboursurface loop. 416 // exiting neighboursurface loop.
436 417
437 if (nneighbours == 1) isaccepted[1 418 if (nneighbours == 1) isaccepted[1] = true;
438 419
439 } // neighboursurface loop end 420 } // neighboursurface loop end
440 421
441 // now, we can accept xx intersection 422 // now, we can accept xx intersection
442 423
443 if (isaccepted[0] && isaccepted[1]) << 424 if (isaccepted[0] == true && isaccepted[1] == true) {
444 { << 425 if (distance[i] < bestdistance) {
445 if (distance[i] < bestdistance) << 426 bestdistance = distance[i];
446 { << 427 gxxbest = gxx[i];
447 bestdistance = distance[i]; << 428 besti = i;
448 gxxbest = gxx[i]; << 429 #ifdef G4SPECSDEBUG
449 #ifdef G4TWISTDEBUG <<
450 besti = i; <<
451 G4cout << " G4VTwistSurface:: 430 G4cout << " G4VTwistSurface::DistanceToIn(p,v): "
452 << " areacode sBoundary 431 << " areacode sBoundary & sBoundary distance = "
453 << fName << " " << dist 432 << fName << " " << distance[i] << G4endl;
454 #endif 433 #endif
455 } 434 }
456 } 435 }
>> 436
457 } // else end 437 } // else end
458 } // intersection loop end 438 } // intersection loop end
459 439
460 gxxbest = bestgxx; 440 gxxbest = bestgxx;
461 441
462 #ifdef G4TWISTDEBUG << 442 #ifdef G4SPECSDEBUG
463 if (besti < 0) << 443 if (besti < 0) {
464 { << 444 G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~" << G4endl;
465 G4cout << "~~~ G4VTwistSurface::Distance <<
466 G4cout << " No intersections " << G 445 G4cout << " No intersections " << G4endl;
467 G4cout << " Name : " << fName << G4 446 G4cout << " Name : " << fName << G4endl;
468 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 447 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
469 } << 448 } else {
470 else << 449 G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) : return ~~~~" << G4endl;
471 { <<
472 G4cout << "~~~ G4VTwistSurface::Distance <<
473 G4cout << " Name, i : " << fName < 450 G4cout << " Name, i : " << fName << " , " << besti << G4endl;
474 G4cout << " gxx[i] : " << gxxbest 451 G4cout << " gxx[i] : " << gxxbest << G4endl;
475 G4cout << " bestdist : " << bestdis 452 G4cout << " bestdist : " << bestdistance << G4endl;
476 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 453 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
477 } 454 }
478 455
479 #endif 456 #endif
480 457
481 return bestdistance; 458 return bestdistance;
482 } 459 }
483 460
484 //============================================ 461 //=====================================================================
485 //* DistanceToOut(p, v) ---------------------- 462 //* DistanceToOut(p, v) -----------------------------------------------
486 463
487 G4double G4VTwistSurface::DistanceToOut(const << 464 G4double G4VTwistSurface::DistanceToOut(const G4ThreeVector &gp,
488 const << 465 const G4ThreeVector &gv,
489 << 466 G4ThreeVector &gxxbest)
490 { 467 {
491 #ifdef G4TWISTDEBUG << 468 #ifdef G4SPECSDEBUG
492 G4cout << "~~~~~ G4VTwistSurface::DistanceT 469 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - Start ~~~~" << G4endl;
493 G4cout << " Name : " << fName << G4end 470 G4cout << " Name : " << fName << G4endl;
494 G4cout << " gp : " << gp << G4endl; 471 G4cout << " gp : " << gp << G4endl;
495 G4cout << " gv : " << gv << G4endl; 472 G4cout << " gv : " << gv << G4endl;
496 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 473 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
497 #endif 474 #endif
498 475
499 G4ThreeVector gxx[G4VSURFACENXX]; 476 G4ThreeVector gxx[G4VSURFACENXX];
500 G4double distance[G4VSURFACENXX]; << 477 G4double distance[G4VSURFACENXX] ;
501 G4int areacode[G4VSURFACENXX]; << 478 G4int areacode[G4VSURFACENXX] ;
502 G4bool isvalid[G4VSURFACENXX]; << 479 G4bool isvalid[G4VSURFACENXX] ;
>> 480 G4int i;
503 481
504 for ( G4int i = 0 ; i<G4VSURFACENXX ; ++i ) << 482 for ( i = 0 ; i<G4VSURFACENXX ; i++ )
505 { 483 {
506 distance[i] = kInfinity ; 484 distance[i] = kInfinity ;
507 areacode[i] = sOutside ; 485 areacode[i] = sOutside ;
508 isvalid[i] = false ; 486 isvalid[i] = false ;
509 } 487 }
510 488
511 G4int nxx; 489 G4int nxx;
512 G4double bestdistance = kInfinity; 490 G4double bestdistance = kInfinity;
>> 491 G4int besti = -1;
513 492
514 nxx = DistanceToSurface(gp, gv, gxx, distan 493 nxx = DistanceToSurface(gp, gv, gxx, distance, areacode,
515 isvalid, kValidateW 494 isvalid, kValidateWithTol);
516 495
517 for (G4int i=0; i<nxx; ++i) << 496 for (i=0; i<nxx; i++) {
518 { << 497 if (!(isvalid[i])) {
519 if (!(isvalid[i])) <<
520 { <<
521 continue; 498 continue;
522 } 499 }
523 500
524 G4ThreeVector normal = GetNormal(gxx[i], 501 G4ThreeVector normal = GetNormal(gxx[i], true);
525 if (normal * gv <= 0) << 502 if (normal * gv <= 0) {
526 { <<
527 // particle goes toword inside of sol 503 // particle goes toword inside of solid, return kInfinity
528 #ifdef G4TWISTDEBUG << 504 #ifdef G4SPECSDEBUG
529 G4cout << " G4VTwistSurface::Dista << 505 G4cout << " G4VTwistSurface::DistanceToOut(p,v): normal*gv < 0, normal "
530 << fName << " " << normal 506 << fName << " " << normal
531 << G4endl; 507 << G4endl;
532 #endif 508 #endif
533 } << 509 } else {
534 else <<
535 { <<
536 // gxx[i] is accepted. 510 // gxx[i] is accepted.
537 if (distance[i] < bestdistance) << 511 if (distance[i] < bestdistance) {
538 { <<
539 bestdistance = distance[i]; 512 bestdistance = distance[i];
540 gxxbest = gxx[i]; 513 gxxbest = gxx[i];
>> 514 besti = i;
541 } 515 }
542 } 516 }
543 } 517 }
544 518
545 #ifdef G4TWISTDEBUG << 519 #ifdef G4SPECSDEBUG
546 if (besti < 0) << 520 if (besti < 0) {
547 { << 521 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - return ~~~" << G4endl;
548 G4cout << "~~ G4VTwistSurface::DistanceT <<
549 G4cout << " No intersections " << 522 G4cout << " No intersections " << G4endl;
550 G4cout << " Name : " << fName < 523 G4cout << " Name : " << fName << G4endl;
551 G4cout << " bestdist : " << bestdis 524 G4cout << " bestdist : " << bestdistance << G4endl;
552 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 525 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
553 } << 526 } else {
554 else << 527 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) : return ~~~" << G4endl;
555 { <<
556 G4cout << "~~ G4VTwistSurface::DistanceT <<
557 G4cout << " Name, i : " << fName < 528 G4cout << " Name, i : " << fName << " , " << i << G4endl;
558 G4cout << " gxx[i] : " << gxxbest 529 G4cout << " gxx[i] : " << gxxbest << G4endl;
559 G4cout << " bestdist : " << bestdis 530 G4cout << " bestdist : " << bestdistance << G4endl;
560 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 531 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
561 } 532 }
562 #endif 533 #endif
563 534
564 return bestdistance; 535 return bestdistance;
565 } 536 }
566 537
567 //============================================ 538 //=====================================================================
568 //* DistanceTo(p) ---------------------------- 539 //* DistanceTo(p) -----------------------------------------------------
569 540
570 G4double G4VTwistSurface::DistanceTo(const G4T << 541 G4double G4VTwistSurface::DistanceTo(const G4ThreeVector &gp,
571 G4T << 542 G4ThreeVector &gxxbest)
572 { 543 {
573 #ifdef G4TWISTDEBUG << 544 #ifdef G4SPECSDEBUG
574 G4cout << "~~~~~ G4VTwistSurface::DistanceT 545 G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - Start ~~~~~~~~~" << G4endl;
575 G4cout << " Name : " << fName << G4end 546 G4cout << " Name : " << fName << G4endl;
576 G4cout << " gp : " << gp << G4endl; 547 G4cout << " gp : " << gp << G4endl;
577 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 548 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
578 #endif 549 #endif
579 550
580 551
581 G4ThreeVector gxx[G4VSURFACENXX]; 552 G4ThreeVector gxx[G4VSURFACENXX];
582 G4double distance[G4VSURFACENXX] ; 553 G4double distance[G4VSURFACENXX] ;
583 G4int areacode[G4VSURFACENXX] ; 554 G4int areacode[G4VSURFACENXX] ;
584 555
585 for (G4int i = 0 ; i<G4VSURFACENXX ; ++i ) << 556 for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) {
586 { <<
587 distance[i] = kInfinity ; 557 distance[i] = kInfinity ;
588 areacode[i] = sOutside ; 558 areacode[i] = sOutside ;
589 } 559 }
590 560
591 DistanceToSurface(gp, gxx, distance, areaco << 561 G4int nxx;
>> 562
>> 563 nxx = DistanceToSurface(gp, gxx, distance, areacode);
592 gxxbest = gxx[0]; 564 gxxbest = gxx[0];
593 565
594 #ifdef G4TWISTDEBUG << 566 #ifdef G4SPECSDEBUG
595 G4cout << "~~~~~ G4VTwistSurface::DistanceT 567 G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - return ~~~~~~~~" << G4endl;
596 G4cout << " Name : " << fName << G 568 G4cout << " Name : " << fName << G4endl;
597 G4cout << " gxx : " << gxxbest << 569 G4cout << " gxx : " << gxxbest << G4endl;
598 G4cout << " bestdist : " << distance[0 570 G4cout << " bestdist : " << distance[0] << G4endl;
599 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 571 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
600 #endif 572 #endif
601 573
602 return distance[0]; 574 return distance[0];
603 } 575 }
604 576
605 //============================================ 577 //=====================================================================
606 //* IsSameBoundary --------------------------- 578 //* IsSameBoundary ----------------------------------------------------
607 579
608 G4bool << 580 G4bool G4VTwistSurface::IsSameBoundary(G4VTwistSurface *surface1, G4int areacode1,
609 G4VTwistSurface::IsSameBoundary(G4VTwistSurfac << 581 G4VTwistSurface *surface2, G4int areacode2 ) const
610 G4VTwistSurfac <<
611 { 582 {
612 // 583 //
613 // IsSameBoundary 584 // IsSameBoundary
614 // 585 //
615 // checking tool whether two boundaries on 586 // checking tool whether two boundaries on different surfaces are same or not.
616 // 587 //
617 588
618 G4bool testbitmode = true; 589 G4bool testbitmode = true;
619 G4bool iscorner[2] = {IsCorner(areacode1, t 590 G4bool iscorner[2] = {IsCorner(areacode1, testbitmode),
620 IsCorner(areacode2, t 591 IsCorner(areacode2, testbitmode)};
621 592
622 if (iscorner[0] && iscorner[1]) << 593 if (iscorner[0] && iscorner[1]) {
623 { <<
624 // on corner 594 // on corner
625 G4ThreeVector corner1 = 595 G4ThreeVector corner1 =
626 surf1->ComputeGlobalPoint(surf1->Ge << 596 surface1->ComputeGlobalPoint(surface1->GetCorner(areacode1));
627 G4ThreeVector corner2 = 597 G4ThreeVector corner2 =
628 surf2->ComputeGlobalPoint(surf2->Ge << 598 surface2->ComputeGlobalPoint(surface2->GetCorner(areacode2));
629 599
630 return (corner1 - corner2).mag() < kCarT << 600 if ((corner1 - corner2).mag() < kCarTolerance) {
631 } << 601 return true;
632 else if ((IsBoundary(areacode1, testbitmode << 602 } else {
633 (IsBoundary(areacode2, testbitmode << 603 return false;
634 { << 604 }
>> 605
>> 606 } else if ((IsBoundary(areacode1, testbitmode) && (!iscorner[0])) &&
>> 607 (IsBoundary(areacode2, testbitmode) && (!iscorner[1]))) {
635 // on boundary 608 // on boundary
636 G4ThreeVector d1, d2, ld1, ld2; 609 G4ThreeVector d1, d2, ld1, ld2;
637 G4ThreeVector x01, x02, lx01, lx02; 610 G4ThreeVector x01, x02, lx01, lx02;
638 G4int type1, type2; 611 G4int type1, type2;
639 surf1->GetBoundaryParameters(areacode1, << 612 surface1->GetBoundaryParameters(areacode1, ld1, lx01, type1);
640 surf2->GetBoundaryParameters(areacode2, << 613 surface2->GetBoundaryParameters(areacode2, ld2, lx02, type2);
641 614
642 x01 = surf1->ComputeGlobalPoint(lx01); << 615 x01 = surface1->ComputeGlobalPoint(lx01);
643 x02 = surf2->ComputeGlobalPoint(lx02); << 616 x02 = surface2->ComputeGlobalPoint(lx02);
644 d1 = surf1->ComputeGlobalDirection(ld1) << 617 d1 = surface1->ComputeGlobalDirection(ld1);
645 d2 = surf2->ComputeGlobalDirection(ld2) << 618 d2 = surface2->ComputeGlobalDirection(ld2);
>> 619
>> 620 if ((x01 - x02).mag() < kCarTolerance &&
>> 621 (d1 - d2).mag() < kCarTolerance) {
>> 622 return true;
>> 623 } else {
>> 624 return false;
>> 625 }
646 626
647 return (x01 - x02).mag() < kCarTolerance << 627 } else {
648 && (d1 - d2).mag() < kCarTolerance <<
649 } <<
650 else <<
651 { <<
652 return false; 628 return false;
653 } 629 }
654 } 630 }
655 631
656 //============================================ 632 //=====================================================================
657 //* GetBoundaryParameters -------------------- 633 //* GetBoundaryParameters ---------------------------------------------
658 634
659 void G4VTwistSurface::GetBoundaryParameters(co << 635 void G4VTwistSurface::GetBoundaryParameters(const G4int &areacode,
660 G << 636 G4ThreeVector &d,
661 G << 637 G4ThreeVector &x0,
662 G << 638 G4int &boundarytype) const
663 { 639 {
664 // areacode must be one of them: 640 // areacode must be one of them:
665 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 641 // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
666 // sAxis1 & sAxisMin, sAxis1 & sAxisMax. 642 // sAxis1 & sAxisMin, sAxis1 & sAxisMax.
667 643
668 for (const auto & boundary : fBoundaries) << 644 G4int i;
669 { << 645 for (i=0; i<4; i++) {
670 if (boundary.GetBoundaryParameters(areac << 646 if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0,
671 { << 647 boundarytype)) {
672 return; 648 return;
673 } 649 }
674 } 650 }
675 651
676 std::ostringstream message; << 652 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryParameters()" << G4endl
677 message << "Not registered boundary." << G4 << 653 << " Boundary at areacode " << std::hex << areacode
678 << " Boundary at areacode " << 654 << std::dec << G4endl
679 << std::dec << G4endl << 655 << " is not be registered." << G4endl;
680 << " is not registered."; << 656 G4Exception("G4VTwistSurface::GetBoundaryParameters()", "InvalidSetup",
681 G4Exception("G4VTwistSurface::GetBoundaryPa << 657 FatalException, "Not registered boundary.");
682 FatalException, message); <<
683 } 658 }
684 659
685 //============================================ 660 //=====================================================================
686 //* GetBoundaryAtPZ -------------------------- 661 //* GetBoundaryAtPZ ---------------------------------------------------
687 662
688 G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ 663 G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ(G4int areacode,
689 << 664 const G4ThreeVector &p) const
690 { 665 {
691 // areacode must be one of them: 666 // areacode must be one of them:
692 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 667 // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
693 // sAxis1 & sAxisMin, sAxis1 & sAxisMax. 668 // sAxis1 & sAxisMin, sAxis1 & sAxisMax.
694 669
695 if (((areacode & sAxis0) != 0) && ((areacod << 670 if (areacode & sAxis0 && areacode & sAxis1) {
696 { << 671 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAtPZ()" << G4endl
697 std::ostringstream message; << 672 << " Point is in the corner area. This function returns"
698 message << "Point is in the corner area." << 673 << G4endl
699 << " This function returns << 674 << " a direction vector of a boundary line." << G4endl
700 << "a direction vector of a bound << 675 << " areacode = " << areacode << G4endl;
701 << " areacode = " << areac << 676 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "InvalidCondition",
702 G4Exception("G4VTwistSurface::GetBoundary << 677 FatalException, "Point is in the corner area.");
703 FatalException, message); <<
704 } 678 }
705 679
706 G4ThreeVector d; 680 G4ThreeVector d;
707 G4ThreeVector x0; 681 G4ThreeVector x0;
708 G4int boundarytype = 0; << 682 G4int boundarytype;
709 G4bool found = false; 683 G4bool found = false;
710 684
711 for (const auto & boundary : fBoundaries) << 685 for (G4int i=0; i<4; i++) {
712 { << 686 if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0,
713 if (boundary.GetBoundaryParameters(areac << 687 boundarytype)){
714 { <<
715 found = true; 688 found = true;
716 continue; 689 continue;
717 } 690 }
718 } 691 }
719 692
720 if (!found) << 693 if (!found) {
721 { << 694 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAtPZ()" << G4endl
722 std::ostringstream message; << 695 << " Boundary at areacode " << areacode << G4endl
723 message << "Not registered boundary." << << 696 << " is not be registered." << G4endl;
724 << " Boundary at areacode << 697 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "InvalidSetup",
725 << " is not registered."; << 698 FatalException, "Not registered boundary.");
726 G4Exception("G4VTwistSurface::GetBoundary <<
727 FatalException, message); <<
728 } 699 }
729 700
730 if (((boundarytype & sAxisPhi) == sAxisPhi) 701 if (((boundarytype & sAxisPhi) == sAxisPhi) ||
731 ((boundarytype & sAxisRho) == sAxisRho) << 702 ((boundarytype & sAxisRho) == sAxisRho)) {
732 { << 703 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAtPZ()" << G4endl
733 std::ostringstream message; << 704 << " Boundary at areacode " << areacode << G4endl
734 message << "Not a z-depended line boundar << 705 << " is not a z-depended line." << G4endl;
735 << " Boundary at areacode << 706 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "InvalidSetup",
736 << " is not a z-depended l << 707 FatalException, "Not a z-depended line boundary.");
737 G4Exception("G4VTwistSurface::GetBoundary <<
738 FatalException, message); <<
739 } 708 }
740 return ((p.z() - x0.z()) / d.z()) * d + x0; 709 return ((p.z() - x0.z()) / d.z()) * d + x0;
741 } 710 }
742 711
743 //============================================ 712 //=====================================================================
744 //* SetCorner -------------------------------- 713 //* SetCorner ---------------------------------------------------------
745 714
746 void G4VTwistSurface::SetCorner(G4int areacode << 715 void G4VTwistSurface::SetCorner(G4int areacode, G4double x, G4double y, G4double z)
747 G4double x, G4 <<
748 { 716 {
749 if ((areacode & sCorner) != sCorner) << 717 if ((areacode & sCorner) != sCorner){
750 { << 718 G4cerr << "ERROR - G4VTwistSurface::SetCorner()" << G4endl
751 std::ostringstream message; << 719 << " areacode " << areacode << G4endl;
752 message << "Area code must represents cor << 720 G4Exception("G4VTwistSurface::SetCorner()", "InvalidSetup",
753 << " areacode " << areacod << 721 FatalException, "Area code must represents corner.");
754 G4Exception("G4VTwistSurface::SetCorner() <<
755 FatalException, message); <<
756 } 722 }
757 723
758 if ((areacode & sC0Min1Min) == sC0Min1Min) 724 if ((areacode & sC0Min1Min) == sC0Min1Min) {
759 fCorners[0].set(x, y, z); 725 fCorners[0].set(x, y, z);
760 } else if ((areacode & sC0Max1Min) == sC0Ma 726 } else if ((areacode & sC0Max1Min) == sC0Max1Min) {
761 fCorners[1].set(x, y, z); 727 fCorners[1].set(x, y, z);
762 } else if ((areacode & sC0Max1Max) == sC0Ma 728 } else if ((areacode & sC0Max1Max) == sC0Max1Max) {
763 fCorners[2].set(x, y, z); 729 fCorners[2].set(x, y, z);
764 } else if ((areacode & sC0Min1Max) == sC0Mi 730 } else if ((areacode & sC0Min1Max) == sC0Min1Max) {
765 fCorners[3].set(x, y, z); 731 fCorners[3].set(x, y, z);
766 } 732 }
767 } 733 }
768 734
769 //============================================ 735 //=====================================================================
770 //* SetBoundaryAxis -------------------------- 736 //* SetBoundaryAxis ---------------------------------------------------
771 737
772 void G4VTwistSurface::GetBoundaryAxis(G4int ar 738 void G4VTwistSurface::GetBoundaryAxis(G4int areacode, EAxis axis[]) const
773 { 739 {
774 if ((areacode & sBoundary) != sBoundary) { 740 if ((areacode & sBoundary) != sBoundary) {
775 G4Exception("G4VTwistSurface::GetBoundary << 741 G4Exception("G4VTwistSurface::GetBoundaryAxis()", "InvalidCondition",
776 FatalException, "Not located 742 FatalException, "Not located on a boundary!");
777 } 743 }
778 for (G4int i=0; i<2; ++i) << 744 G4int i;
779 { << 745 for (i=0; i<2; i++) {
>> 746
780 G4int whichaxis = 0 ; 747 G4int whichaxis = 0 ;
781 if (i == 0) { 748 if (i == 0) {
782 whichaxis = sAxis0; 749 whichaxis = sAxis0;
783 } else if (i == 1) { 750 } else if (i == 1) {
784 whichaxis = sAxis1; 751 whichaxis = sAxis1;
785 } 752 }
786 753
787 // extracted axiscode of whichaxis 754 // extracted axiscode of whichaxis
788 G4int axiscode = whichaxis & sAxisMask & 755 G4int axiscode = whichaxis & sAxisMask & areacode ;
789 if (axiscode != 0) { << 756 if (axiscode) {
790 if (axiscode == (whichaxis & sAxisX)) 757 if (axiscode == (whichaxis & sAxisX)) {
791 axis[i] = kXAxis; 758 axis[i] = kXAxis;
792 } else if (axiscode == (whichaxis & s 759 } else if (axiscode == (whichaxis & sAxisY)) {
793 axis[i] = kYAxis; 760 axis[i] = kYAxis;
794 } else if (axiscode == (whichaxis & s 761 } else if (axiscode == (whichaxis & sAxisZ)) {
795 axis[i] = kZAxis; 762 axis[i] = kZAxis;
796 } else if (axiscode == (whichaxis & s 763 } else if (axiscode == (whichaxis & sAxisRho)) {
797 axis[i] = kRho; 764 axis[i] = kRho;
798 } else if (axiscode == (whichaxis & s 765 } else if (axiscode == (whichaxis & sAxisPhi)) {
799 axis[i] = kPhi; 766 axis[i] = kPhi;
800 } else { 767 } else {
801 std::ostringstream message; << 768 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAxis()" << G4endl
802 message << "Not supported areacode. << 769 << " areacode " << areacode << G4endl;
803 << " areacode " << a << 770 G4Exception("G4VTwistSurface::GetBoundaryAxis()", "InvalidSetup",
804 G4Exception("G4VTwistSurface::GetBo << 771 FatalException, "Not supported areacode.");
805 FatalException, message <<
806 } 772 }
807 } 773 }
808 } 774 }
809 } 775 }
810 776
811 //============================================ 777 //=====================================================================
812 //* SetBoundaryLimit ------------------------- 778 //* SetBoundaryLimit --------------------------------------------------
813 779
814 void G4VTwistSurface::GetBoundaryLimit(G4int a 780 void G4VTwistSurface::GetBoundaryLimit(G4int areacode, G4double limit[]) const
815 { 781 {
816 if ((areacode & sCorner) != 0) { << 782 if (areacode & sCorner) {
817 if ((areacode & sC0Min1Min) != 0) { << 783 if (areacode & sC0Min1Max) {
818 limit[0] = fAxisMin[0]; 784 limit[0] = fAxisMin[0];
819 limit[1] = fAxisMin[1]; 785 limit[1] = fAxisMin[1];
820 } else if ((areacode & sC0Max1Min) != 0) << 786 } else if (areacode & sC0Max1Min) {
821 limit[0] = fAxisMax[0]; 787 limit[0] = fAxisMax[0];
822 limit[1] = fAxisMin[1]; 788 limit[1] = fAxisMin[1];
823 } else if ((areacode & sC0Max1Max) != 0) << 789 } else if (areacode & sC0Max1Max) {
824 limit[0] = fAxisMax[0]; 790 limit[0] = fAxisMax[0];
825 limit[1] = fAxisMax[1]; 791 limit[1] = fAxisMax[1];
826 } else if ((areacode & sC0Min1Max) != 0) << 792 } else if (areacode & sC0Min1Max) {
827 limit[0] = fAxisMin[0]; 793 limit[0] = fAxisMin[0];
828 limit[1] = fAxisMax[1]; 794 limit[1] = fAxisMax[1];
829 } 795 }
830 } else if ((areacode & sBoundary) != 0) { << 796 } else if (areacode & sBoundary) {
831 if ((areacode & (sAxis0 | sAxisMin)) != << 797 if (areacode & (sAxis0 | sAxisMin)) {
832 limit[0] = fAxisMin[0]; 798 limit[0] = fAxisMin[0];
833 } else if ((areacode & (sAxis1 | sAxisMi << 799 } else if (areacode & (sAxis1 | sAxisMin)) {
834 limit[0] = fAxisMin[1]; 800 limit[0] = fAxisMin[1];
835 } else if ((areacode & (sAxis0 | sAxisMa << 801 } else if (areacode & (sAxis0 | sAxisMax)) {
836 limit[0] = fAxisMax[0]; 802 limit[0] = fAxisMax[0];
837 } else if ((areacode & (sAxis1 | sAxisMa << 803 } else if (areacode & (sAxis1 | sAxisMax)) {
838 limit[0] = fAxisMax[1]; 804 limit[0] = fAxisMax[1];
839 } 805 }
840 } else { 806 } else {
841 std::ostringstream message; << 807 G4cerr << "WARNING - G4VTwistSurface::GetBoundaryAxis()" << G4endl
842 message << "Not located on a boundary!" < << 808 << " areacode " << areacode << G4endl;
843 << " areacode " << areac << 809 G4Exception("G4VTwistSurface::GetBoundaryLimit()", "InvalidCondition",
844 G4Exception("G4VTwistSurface::GetBoundary << 810 JustWarning, "Not located on a boundary!");
845 JustWarning, message); <<
846 } 811 }
847 } 812 }
848 813
849 //============================================ 814 //=====================================================================
850 //* SetBoundary ------------------------------ 815 //* SetBoundary -------------------------------------------------------
851 816
852 void G4VTwistSurface::SetBoundary(const G4int& << 817 void G4VTwistSurface::SetBoundary(const G4int &axiscode,
853 const G4Thre << 818 const G4ThreeVector &direction,
854 const G4Thre << 819 const G4ThreeVector &x0,
855 const G4int& << 820 const G4int &boundarytype)
856 { 821 {
857 G4int code = (~sAxisMask) & axiscode; 822 G4int code = (~sAxisMask) & axiscode;
858 if ((code == (sAxis0 & sAxisMin)) || 823 if ((code == (sAxis0 & sAxisMin)) ||
859 (code == (sAxis0 & sAxisMax)) || 824 (code == (sAxis0 & sAxisMax)) ||
860 (code == (sAxis1 & sAxisMin)) || 825 (code == (sAxis1 & sAxisMin)) ||
861 (code == (sAxis1 & sAxisMax))) << 826 (code == (sAxis1 & sAxisMax))) {
862 { << 827
>> 828 G4int i;
863 G4bool done = false; 829 G4bool done = false;
864 for (auto & boundary : fBoundaries) << 830 for (i=0; i<4; i++) {
865 { << 831 if (fBoundaries[i].IsEmpty()) {
866 if (boundary.IsEmpty()) << 832 fBoundaries[i].SetFields(axiscode, direction,
867 { << 833 x0, boundarytype);
868 boundary.SetFields(axiscode, direc <<
869 done = true; 834 done = true;
870 break; 835 break;
871 } 836 }
872 } 837 }
873 838
874 if (!done) << 839 if (!done) {
875 { << 840 G4Exception("G4VTwistSurface::SetBoundary()", "InvalidCondition",
876 G4Exception("G4VTwistSurface::SetBoun <<
877 FatalException, "Number 841 FatalException, "Number of boundary exceeding 4!");
878 } 842 }
879 } << 843 } else {
880 else << 844 G4cerr << "ERROR - G4VTwistSurface::SetBoundary()" << G4endl
881 { << 845 << " invalid axiscode. axiscode = "
882 std::ostringstream message; << 846 << std::hex << axiscode << std::dec << G4endl;
883 message << "Invalid axis-code." << G4end << 847 G4Exception("G4VTwistSurface::SetBoundary()", "InvalidCondition",
884 << " axiscode = " << 848 FatalException, "Invalid axis-code.");
885 << std::hex << axiscode << std:: <<
886 G4Exception("G4VTwistSurface::SetBoundar <<
887 FatalException, message); <<
888 } 849 }
889 } 850 }
890 851
891 //============================================ 852 //=====================================================================
892 //* GetFace ---------------------------------- 853 //* GetFace -----------------------------------------------------------
893 854
894 G4int G4VTwistSurface::GetFace( G4int i, G4int << 855 G4int G4VTwistSurface::GetFace( G4int i, G4int j, G4int m,
895 G4int n, G4int 856 G4int n, G4int iside )
896 { 857 {
897 // this is the face mapping function 858 // this is the face mapping function
898 // (i,j) -> face number 859 // (i,j) -> face number
899 860
900 if ( iside == 0 ) { 861 if ( iside == 0 ) {
901 return i * ( k - 1 ) + j ; << 862 return i * ( m - 1 ) + j ;
902 } 863 }
903 864
904 else if ( iside == 1 ) { 865 else if ( iside == 1 ) {
905 return (k-1)*(k-1) + i*(k-1) + j ; << 866 return (m-1)*(m-1) + i*(m-1) + j ;
906 } 867 }
907 868
908 else if ( iside == 2 ) { 869 else if ( iside == 2 ) {
909 return 2*(k-1)*(k-1) + i*(k-1) + j ; << 870 return 2*(m-1)*(m-1) + i*(m-1) + j ;
910 } 871 }
911 872
912 else if ( iside == 3 ) { 873 else if ( iside == 3 ) {
913 return 2*(k-1)*(k-1) + (n-1)*(k-1) + i*(k- << 874 return 2*(m-1)*(m-1) + (n-1)*(m-1) + i*(m-1) + j ;
914 } 875 }
915 876
916 else if ( iside == 4 ) { 877 else if ( iside == 4 ) {
917 return 2*(k-1)*(k-1) + 2*(n-1)*(k-1) + i*( << 878 return 2*(m-1)*(m-1) + 2*(n-1)*(m-1) + i*(m-1) + j ;
918 } 879 }
919 880
920 else if ( iside == 5 ) { 881 else if ( iside == 5 ) {
921 return 2*(k-1)*(k-1) + 3*(n-1)*(k-1) + i*( << 882 return 2*(m-1)*(m-1) + 3*(n-1)*(m-1) + i*(m-1) + j ;
922 } 883 }
923 884
924 else { 885 else {
925 std::ostringstream message; << 886
926 message << "Not correct side number: " << 887 G4cerr << "ERROR - G4VTwistSurface::GetFace(): "
927 << GetName() << G4endl << 888 << GetName() << G4endl
928 << "iside is " << iside << " but s << 889 << " Not correct side number! - " << G4endl
929 << "0,1,2,3,4 or 5" << "."; << 890 << "iside is " << iside << " but should be "
930 G4Exception("G4TwistSurface::G4GetFace()", << 891 << "0,1,2,3,4 or 5" << "." << G4endl ;
931 FatalException, message); << 892 G4Exception("G4TwistSurface::G4GetFace()", "InvalidSetup",
>> 893 FatalException, "Not correct side number.");
>> 894
>> 895
932 } 896 }
933 897
934 return -1 ; // wrong face 898 return -1 ; // wrong face
935 } 899 }
936 900
937 //============================================ 901 //=====================================================================
938 //* GetNode ---------------------------------- 902 //* GetNode -----------------------------------------------------------
939 903
940 G4int G4VTwistSurface::GetNode( G4int i, G4int << 904 G4int G4VTwistSurface::GetNode( G4int i, G4int j, G4int m,
941 G4int n, G4int 905 G4int n, G4int iside )
942 { 906 {
943 // this is the node mapping function 907 // this is the node mapping function
944 // (i,j) -> node number 908 // (i,j) -> node number
945 // Depends on the side iside and the used me 909 // Depends on the side iside and the used meshing of the surface
946 910
947 if ( iside == 0 ) << 911 if ( iside == 0 ) {
948 { << 912 // lower endcap is mxm squared.
949 // lower endcap is kxk squared. << 913 // n = m
950 // n = k << 914 return i * m + j ;
951 return i * k + j ; <<
952 } 915 }
953 916
954 if ( iside == 1 ) << 917 if ( iside == 1 ) {
955 { << 918 // upper endcap is mxm squared. Shift by m*m
956 // upper endcap is kxk squared. Shift by k << 919 // n = m
957 // n = k << 920 return m*m + i*m + j ;
958 return k*k + i*k + j ; <<
959 } 921 }
960 922
961 else if ( iside == 2 ) << 923 else if ( iside == 2 ) {
962 { <<
963 // front side. 924 // front side.
964 if ( i == 0 ) { return j ; } << 925 if ( i == 0 ) { return j ; }
965 else if ( i == n-1 ) { return k*k + j ; } << 926 else if ( i == n-1 ) { return m*m + j ; }
966 else { return 2*k*k + 4*(i << 927 else { return 2*m*m + 4*(i-1)*(m-1) + j ; }
967 } 928 }
968 929
969 else if ( iside == 3 ) << 930 else if ( iside == 3 ) {
970 { <<
971 // right side 931 // right side
972 if ( i == 0 ) { return (j+1)* << 932 if ( i == 0 ) { return (j+1)*m - 1 ; }
973 else if ( i == n-1 ) { return k*k + (j+1)* << 933 else if ( i == n-1 ) { return m*m + (j+1)*m - 1 ; }
974 else << 934 else { return 2*m*m + 4*(i-1)*(m-1) + (m-1) + j ; }
975 { <<
976 return 2*k*k + 4*(i-1)*(k-1) + (k-1) + j <<
977 } <<
978 } 935 }
979 else if ( iside == 4 ) << 936 else if ( iside == 4 ) {
980 { << 937 // back side.
981 // back side << 938 if ( i == 0 ) { return m*m - 1 - j ; } // reversed order
982 if ( i == 0 ) { return k*k - 1 - << 939 else if ( i == n-1 ) { return 2*m*m - 1 - j ; } // reversed order
983 else if ( i == n-1 ) { return 2*k*k - 1 - << 940 else { return 2*m*m + 4*(i-1)*(m-1) + 2*(m-1) + j ; // normal order
984 else <<
985 { <<
986 return 2*k*k + 4*(i-1)*(k-1) + 2*(k-1) + <<
987 } 941 }
988 } 942 }
989 else if ( iside == 5 ) << 943 else if ( iside == 5 ) {
990 { <<
991 // left side 944 // left side
992 if ( i == 0 ) { return k*k - (j+1 << 945 if ( i == 0 ) { return m*m - (j+1)*m ; } // reversed order
993 else if ( i == n-1) { return 2*k*k - (j+1 << 946 else if ( i == n-1) { return 2*m*m - (j+1)*m ; } // reverded order
994 else << 947 else {
995 { << 948 if ( j == m-1 ) { return 2*m*m + 4*(i-1)*(m-1) ; } // special case
996 if ( j == k-1 ) { return 2*k*k + 4*(i-1) << 949 else { return 2*m*m + 4*(i-1)*(m-1) + 3*(m-1) + j ; } // normal order
997 else << 950 }
998 { <<
999 return 2*k*k + 4*(i-1)*(k-1) + 3*(k-1) <<
1000 } <<
1001 } <<
1002 } 951 }
1003 else << 952
1004 { << 953 else {
1005 std::ostringstream message; << 954
1006 message << "Not correct side number: " << 955 G4cerr << "ERROR - G4VTwistSurface::GetNode(): "
1007 << GetName() << G4endl << 956 << GetName() << G4endl
1008 << "iside is " << iside << " but << 957 << " Not correct side number! - " << G4endl
1009 << "0,1,2,3,4 or 5" << "."; << 958 << "iside is " << iside << " but should be "
1010 G4Exception("G4TwistSurface::G4GetNode()" << 959 << "0,1,2,3,4 or 5" << "." << G4endl ;
1011 FatalException, message); << 960 G4Exception("G4TwistSurface::G4GetNode()", "InvalidSetup",
>> 961 FatalException, "Not correct side number.");
1012 } 962 }
1013 return -1 ; // wrong node 963 return -1 ; // wrong node
1014 } 964 }
1015 965
1016 //=========================================== 966 //=====================================================================
1017 //* GetEdgeVisiblility ---------------------- 967 //* GetEdgeVisiblility ------------------------------------------------
1018 968
1019 G4int G4VTwistSurface::GetEdgeVisibility( G4i << 969 G4int G4VTwistSurface::GetEdgeVisibility( G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
1020 G4i <<
1021 { 970 {
>> 971
1022 // clockwise filling -> positive or 972 // clockwise filling -> positive orientation
1023 // counter clockwise filling -> negative or 973 // counter clockwise filling -> negative orientation
1024 974
1025 // 975 //
1026 // d C c 976 // d C c
1027 // +------+ 977 // +------+
1028 // | | 978 // | |
1029 // | | 979 // | |
1030 // | | 980 // | |
1031 // D | |B 981 // D | |B
1032 // | | 982 // | |
1033 // | | 983 // | |
1034 // | | 984 // | |
1035 // +------+ 985 // +------+
1036 // a A b 986 // a A b
1037 // 987 //
1038 // a = +--+ A = ---+ 988 // a = +--+ A = ---+
1039 // b = --++ B = --+- 989 // b = --++ B = --+-
1040 // c = -++- C = -+-- 990 // c = -++- C = -+--
1041 // d = ++-- D = +--- 991 // d = ++-- D = +---
1042 992
1043 993
1044 // check first invisible faces 994 // check first invisible faces
1045 995
1046 if ( ( i>0 && i<n-2 ) && ( j>0 && j<k-2 ) ) << 996 if ( ( i>0 && i<n-2 ) && ( j>0 && j<m-2 ) ) {
1047 { <<
1048 return -1 ; // always invisible, signs: 997 return -1 ; // always invisible, signs: ----
1049 } 998 }
1050 999
1051 // change first the vertex number (depends 1000 // change first the vertex number (depends on the orientation)
1052 // 0,1,2,3 -> 3,2,1,0 1001 // 0,1,2,3 -> 3,2,1,0
1053 if ( orientation < 0 ) { number = ( 3 - num 1002 if ( orientation < 0 ) { number = ( 3 - number ) ; }
1054 1003
1055 // check true edges 1004 // check true edges
1056 if ( ( j>=1 && j<=k-3 ) ) << 1005 if ( ( j>=1 && j<=m-3 ) ) {
1057 { << 1006
1058 if ( i == 0 ) { // signs (A): -- 1007 if ( i == 0 ) { // signs (A): ---+
1059 return ( number == 3 ) ? 1 : -1 ; 1008 return ( number == 3 ) ? 1 : -1 ;
1060 } 1009 }
1061 1010
1062 else if ( i == n-2 ) { // signs (C): -+ 1011 else if ( i == n-2 ) { // signs (C): -+--
1063 return ( number == 1 ) ? 1 : -1 ; 1012 return ( number == 1 ) ? 1 : -1 ;
1064 } 1013 }
1065 1014
1066 else << 1015 else {
1067 { << 1016 G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): "
1068 std::ostringstream message; << 1017 << GetName() << G4endl
1069 message << "Not correct face number: " << 1018 << " Not correct face number! - " << G4endl ;
1070 G4Exception("G4TwistSurface::G4GetEdgeV << 1019 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup",
1071 "GeomSolids0003", FatalExce << 1020 FatalException, "Not correct face number.");
1072 } 1021 }
1073 } 1022 }
1074 1023
1075 if ( ( i>=1 && i<=n-3 ) ) << 1024 if ( ( i>=1 && i<=n-3 ) ) {
1076 { << 1025
1077 if ( j == 0 ) { // signs (D): +- 1026 if ( j == 0 ) { // signs (D): +---
1078 return ( number == 0 ) ? 1 : -1 ; 1027 return ( number == 0 ) ? 1 : -1 ;
1079 } 1028 }
1080 1029
1081 else if ( j == k-2 ) { // signs (B): -- << 1030 else if ( j == m-2 ) { // signs (B): --+-
1082 return ( number == 2 ) ? 1 : -1 ; 1031 return ( number == 2 ) ? 1 : -1 ;
1083 } 1032 }
1084 1033
1085 else << 1034 else {
1086 { << 1035 G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): "
1087 std::ostringstream message; << 1036 << GetName() << G4endl
1088 message << "Not correct face number: " << 1037 << " Not correct face number! - " << G4endl ;
1089 G4Exception("G4TwistSurface::G4GetEdgeV << 1038 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup",
1090 "GeomSolids0003", FatalExce << 1039 FatalException, "Not correct face number.");
1091 } 1040 }
1092 } 1041 }
1093 1042
1094 // now the corners 1043 // now the corners
1095 if ( i == 0 && j == 0 ) { // signs 1044 if ( i == 0 && j == 0 ) { // signs (a) : +--+
1096 return ( number == 0 || number == 3 ) ? 1 1045 return ( number == 0 || number == 3 ) ? 1 : -1 ;
1097 } 1046 }
1098 else if ( i == 0 && j == k-2 ) { // signs << 1047 else if ( i == 0 && j == m-2 ) { // signs (b) : --++
1099 return ( number == 2 || number == 3 ) ? 1 1048 return ( number == 2 || number == 3 ) ? 1 : -1 ;
1100 } 1049 }
1101 else if ( i == n-2 && j == k-2 ) { // signs << 1050 else if ( i == n-2 && j == m-2 ) { // signs (c) : -++-
1102 return ( number == 1 || number == 2 ) ? 1 1051 return ( number == 1 || number == 2 ) ? 1 : -1 ;
1103 } 1052 }
1104 else if ( i == n-2 && j == 0 ) { // signs << 1053 else if ( i == n-2 && j == j ) { // signs (d) : ++--
1105 return ( number == 0 || number == 1 ) ? 1 1054 return ( number == 0 || number == 1 ) ? 1 : -1 ;
1106 } 1055 }
1107 else << 1056 else {
1108 { << 1057 G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): "
1109 std::ostringstream message; << 1058 << GetName() << G4endl
1110 message << "Not correct face number: " << << 1059 << " Not correct face number! - " << G4endl ;
1111 G4Exception("G4TwistSurface::G4GetEdgeVis << 1060 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup",
1112 "GeomSolids0003", FatalExcept << 1061 FatalException, "Not correct face number.");
1113 } 1062 }
1114 1063
1115 std::ostringstream message; << 1064 G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): "
1116 message << "Not correct face number: " << G << 1065 << GetName() << G4endl
1117 G4Exception("G4TwistSurface::G4GetEdgeVisib << 1066 << " Not correct face number! - " << G4endl ;
1118 FatalException, message); << 1067 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup",
1119 << 1068 FatalException, "Not correct face number.");
>> 1069
1120 return 0 ; 1070 return 0 ;
>> 1071
1121 } 1072 }
1122 1073
1123 1074
1124 //=========================================== 1075 //=====================================================================
1125 //* DebugPrint ------------------------------ 1076 //* DebugPrint --------------------------------------------------------
1126 1077
1127 void G4VTwistSurface::DebugPrint() const 1078 void G4VTwistSurface::DebugPrint() const
1128 { 1079 {
1129 G4ThreeVector A = fRot * GetCorner(sC0Min1 1080 G4ThreeVector A = fRot * GetCorner(sC0Min1Min) + fTrans;
1130 G4ThreeVector B = fRot * GetCorner(sC0Max1 1081 G4ThreeVector B = fRot * GetCorner(sC0Max1Min) + fTrans;
1131 G4ThreeVector C = fRot * GetCorner(sC0Max1 1082 G4ThreeVector C = fRot * GetCorner(sC0Max1Max) + fTrans;
1132 G4ThreeVector D = fRot * GetCorner(sC0Min1 1083 G4ThreeVector D = fRot * GetCorner(sC0Min1Max) + fTrans;
1133 1084
1134 G4cout << "/* G4VTwistSurface::DebugPrint( << 1085 G4cout << "/* G4VTwistSurface::DebugPrint():-------------------------------"
1135 << G4endl; 1086 << G4endl;
1136 G4cout << "/* Name = " << fName << G4endl; 1087 G4cout << "/* Name = " << fName << G4endl;
1137 G4cout << "/* Axis = " << std::hex << fAxi 1088 G4cout << "/* Axis = " << std::hex << fAxis[0] << " "
1138 << std::hex << fAxis[1] 1089 << std::hex << fAxis[1]
1139 << " (0,1,2,3,5 = kXAxis,kYAxis,kZA 1090 << " (0,1,2,3,5 = kXAxis,kYAxis,kZAxis,kRho,kPhi)"
1140 << std::dec << G4endl; 1091 << std::dec << G4endl;
1141 G4cout << "/* BoundaryLimit(in local) fAxi 1092 G4cout << "/* BoundaryLimit(in local) fAxis0(min, max) = ("<<fAxisMin[0]
1142 << ", " << fAxisMax[0] << ")" << G4 1093 << ", " << fAxisMax[0] << ")" << G4endl;
1143 G4cout << "/* BoundaryLimit(in local) fAxi 1094 G4cout << "/* BoundaryLimit(in local) fAxis1(min, max) = ("<<fAxisMin[1]
1144 << ", " << fAxisMax[1] << ")" << G4 1095 << ", " << fAxisMax[1] << ")" << G4endl;
1145 G4cout << "/* Cornar point sC0Min1Min = " 1096 G4cout << "/* Cornar point sC0Min1Min = " << A << G4endl;
1146 G4cout << "/* Cornar point sC0Max1Min = " 1097 G4cout << "/* Cornar point sC0Max1Min = " << B << G4endl;
1147 G4cout << "/* Cornar point sC0Max1Max = " 1098 G4cout << "/* Cornar point sC0Max1Max = " << C << G4endl;
1148 G4cout << "/* Cornar point sC0Min1Max = " 1099 G4cout << "/* Cornar point sC0Min1Max = " << D << G4endl;
1149 G4cout << "/*----------------------------- 1100 G4cout << "/*---------------------------------------------------------"
1150 << G4endl; 1101 << G4endl;
1151 } 1102 }
1152 1103
1153 //=========================================== 1104 //=====================================================================
1154 // G4VTwistSurface::CurrentStatus class 1105 // G4VTwistSurface::CurrentStatus class
1155 //=========================================== 1106 //=====================================================================
1156 1107
1157 //=========================================== 1108 //=====================================================================
1158 //* CurrentStatus::CurrentStatus ------------ 1109 //* CurrentStatus::CurrentStatus --------------------------------------
1159 1110
1160 G4VTwistSurface::CurrentStatus::CurrentStatus 1111 G4VTwistSurface::CurrentStatus::CurrentStatus()
1161 { 1112 {
1162 for (size_t i=0; i<G4VSURFACENXX; ++i) << 1113 for (size_t i=0; i<G4VSURFACENXX; i++)
1163 { 1114 {
1164 fDistance[i] = kInfinity; 1115 fDistance[i] = kInfinity;
1165 fAreacode[i] = sOutside; 1116 fAreacode[i] = sOutside;
1166 fIsValid[i] = false; 1117 fIsValid[i] = false;
1167 fXX[i].set(kInfinity, kInfinity, kInfinit 1118 fXX[i].set(kInfinity, kInfinity, kInfinity);
1168 } 1119 }
1169 fNXX = 0; 1120 fNXX = 0;
1170 fLastp.set(kInfinity, kInfinity, kInfinity) 1121 fLastp.set(kInfinity, kInfinity, kInfinity);
1171 fLastv.set(kInfinity, kInfinity, kInfinity) 1122 fLastv.set(kInfinity, kInfinity, kInfinity);
1172 fLastValidate = kUninitialized; 1123 fLastValidate = kUninitialized;
1173 fDone = false; 1124 fDone = false;
1174 } 1125 }
1175 1126
1176 //=========================================== 1127 //=====================================================================
1177 //* CurrentStatus::~CurrentStatus ----------- 1128 //* CurrentStatus::~CurrentStatus -------------------------------------
1178 1129
1179 G4VTwistSurface::CurrentStatus::~CurrentStatu 1130 G4VTwistSurface::CurrentStatus::~CurrentStatus()
1180 = default; << 1131 {
>> 1132 }
1181 1133
1182 //=========================================== 1134 //=====================================================================
1183 //* CurrentStatus::SetCurrentStatus --------- 1135 //* CurrentStatus::SetCurrentStatus -----------------------------------
1184 1136
1185 void 1137 void
1186 G4VTwistSurface::CurrentStatus::SetCurrentSta << 1138 G4VTwistSurface::CurrentStatus::SetCurrentStatus(G4int i,
1187 << 1139 G4ThreeVector &xx,
1188 << 1140 G4double &dist,
1189 << 1141 G4int &areacode,
1190 << 1142 G4bool &isvalid,
1191 << 1143 G4int nxx,
1192 << 1144 EValidate validate,
1193 co << 1145 const G4ThreeVector *p,
1194 co << 1146 const G4ThreeVector *v)
1195 { 1147 {
1196 fDistance[i] = dist; 1148 fDistance[i] = dist;
1197 fAreacode[i] = areacode; 1149 fAreacode[i] = areacode;
1198 fIsValid[i] = isvalid; 1150 fIsValid[i] = isvalid;
1199 fXX[i] = xx; 1151 fXX[i] = xx;
1200 fNXX = nxx; 1152 fNXX = nxx;
1201 fLastValidate = validate; 1153 fLastValidate = validate;
1202 if (p != nullptr) << 1154 if (p)
1203 { 1155 {
1204 fLastp = *p; 1156 fLastp = *p;
1205 } 1157 }
1206 else 1158 else
1207 { 1159 {
1208 G4Exception("G4VTwistSurface::CurrentStat << 1160 G4Exception("G4VTwistSurface::CurrentStatus::CurrentStatus()",
1209 "GeomSolids0003", FatalExcept << 1161 "InvalidCondition", FatalException,
>> 1162 "SetCurrentStatus: p = 0!");
1210 } 1163 }
1211 if (v != nullptr) << 1164 if (v)
1212 { 1165 {
1213 fLastv = *v; 1166 fLastv = *v;
1214 } 1167 }
1215 else 1168 else
1216 { 1169 {
1217 fLastv.set(kInfinity, kInfinity, kInfinit 1170 fLastv.set(kInfinity, kInfinity, kInfinity);
1218 } 1171 }
1219 fDone = true; 1172 fDone = true;
1220 } 1173 }
1221 1174
1222 //=========================================== 1175 //=====================================================================
1223 //* CurrentStatus::ResetfDone --------------- 1176 //* CurrentStatus::ResetfDone -----------------------------------------
1224 1177
1225 void 1178 void
1226 G4VTwistSurface::CurrentStatus::ResetfDone(EV << 1179 G4VTwistSurface::CurrentStatus::ResetfDone(EValidate validate,
1227 const G4 << 1180 const G4ThreeVector *p,
1228 const G4 << 1181 const G4ThreeVector *v)
1229 1182
1230 { 1183 {
1231 if (validate == fLastValidate && p != nullp << 1184 if (validate == fLastValidate && p && *p == fLastp)
1232 { 1185 {
1233 if (v == nullptr || (*v == fLastv)) retu << 1186 if (!v || (v && *v == fLastv)) return;
1234 } 1187 }
1235 G4ThreeVector xx(kInfinity, kInfinity, kInf 1188 G4ThreeVector xx(kInfinity, kInfinity, kInfinity);
1236 for (size_t i=0; i<G4VSURFACENXX; ++i) << 1189 for (size_t i=0; i<G4VSURFACENXX; i++)
1237 { 1190 {
1238 fDistance[i] = kInfinity; 1191 fDistance[i] = kInfinity;
1239 fAreacode[i] = sOutside; 1192 fAreacode[i] = sOutside;
1240 fIsValid[i] = false; 1193 fIsValid[i] = false;
1241 fXX[i] = xx; // bug in old code ( was f 1194 fXX[i] = xx; // bug in old code ( was fXX[i] = xx[i] )
1242 } 1195 }
1243 fNXX = 0; << 1196 fNXX = 0;
1244 fLastp.set(kInfinity, kInfinity, kInfinity) 1197 fLastp.set(kInfinity, kInfinity, kInfinity);
1245 fLastv.set(kInfinity, kInfinity, kInfinity) 1198 fLastv.set(kInfinity, kInfinity, kInfinity);
1246 fLastValidate = kUninitialized; 1199 fLastValidate = kUninitialized;
1247 fDone = false; 1200 fDone = false;
1248 } 1201 }
1249 1202
1250 //=========================================== 1203 //=====================================================================
1251 //* CurrentStatus::DebugPrint --------------- 1204 //* CurrentStatus::DebugPrint -----------------------------------------
1252 1205
1253 void 1206 void
1254 G4VTwistSurface::CurrentStatus::DebugPrint() 1207 G4VTwistSurface::CurrentStatus::DebugPrint() const
1255 { 1208 {
1256 G4cout << "CurrentStatus::Dist0,1= " << fDi 1209 G4cout << "CurrentStatus::Dist0,1= " << fDistance[0]
1257 << " " << fDistance[1] << " areacode 1210 << " " << fDistance[1] << " areacode = " << fAreacode[0]
1258 << " " << fAreacode[1] << G4endl; 1211 << " " << fAreacode[1] << G4endl;
1259 } 1212 }
1260 1213
1261 //=========================================== 1214 //=====================================================================
1262 // G4VTwistSurface::Boundary class 1215 // G4VTwistSurface::Boundary class
1263 //=========================================== 1216 //=====================================================================
1264 1217
1265 //=========================================== 1218 //=====================================================================
>> 1219 //* Boundary::Boundary ------------------------------------------------
>> 1220
>> 1221 G4VTwistSurface::Boundary::Boundary()
>> 1222 : fBoundaryAcode(-1), fBoundaryType(0)
>> 1223 {
>> 1224 }
>> 1225
>> 1226 //=====================================================================
>> 1227 //* Boundary::~Boundary -----------------------------------------------
>> 1228
>> 1229 G4VTwistSurface::Boundary::~Boundary()
>> 1230 {
>> 1231 }
>> 1232
>> 1233 //=====================================================================
1266 //* Boundary::SetFields --------------------- 1234 //* Boundary::SetFields -----------------------------------------------
1267 1235
1268 void 1236 void
1269 G4VTwistSurface::Boundary::SetFields(const G4 << 1237 G4VTwistSurface::Boundary::SetFields(const G4int &areacode,
1270 const G4Three << 1238 const G4ThreeVector &d,
1271 const G4Three << 1239 const G4ThreeVector &x0,
1272 const G4int& << 1240 const G4int &boundarytype)
1273 { 1241 {
1274 fBoundaryAcode = areacode; 1242 fBoundaryAcode = areacode;
1275 fBoundaryDirection = d; 1243 fBoundaryDirection = d;
1276 fBoundaryX0 = x0; 1244 fBoundaryX0 = x0;
1277 fBoundaryType = boundarytype; 1245 fBoundaryType = boundarytype;
1278 } 1246 }
1279 1247
1280 //=========================================== 1248 //=====================================================================
1281 //* Boundary::IsEmpty ----------------------- 1249 //* Boundary::IsEmpty -------------------------------------------------
1282 1250
1283 G4bool G4VTwistSurface::Boundary::IsEmpty() c 1251 G4bool G4VTwistSurface::Boundary::IsEmpty() const
1284 { 1252 {
1285 return fBoundaryAcode == -1; << 1253 if (fBoundaryAcode == -1) return true;
>> 1254 return false;
1286 } 1255 }
1287 1256
1288 //=========================================== 1257 //=====================================================================
1289 //* Boundary::GetBoundaryParameters --------- 1258 //* Boundary::GetBoundaryParameters -----------------------------------
1290 1259
1291 G4bool 1260 G4bool
1292 G4VTwistSurface::Boundary::GetBoundaryParamet << 1261 G4VTwistSurface::Boundary::GetBoundaryParameters(const G4int &areacode,
1293 << 1262 G4ThreeVector &d,
1294 << 1263 G4ThreeVector &x0,
1295 << 1264 G4int &boundarytype) const
1296 { 1265 {
1297 // areacode must be one of them: 1266 // areacode must be one of them:
1298 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 1267 // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
1299 // sAxis1 & sAxisMin, sAxis1 & sAxisMax 1268 // sAxis1 & sAxisMin, sAxis1 & sAxisMax
1300 // << 1269 if ((areacode & sAxis0) && (areacode & sAxis1))
1301 if (((areacode & sAxis0) != 0) && ((areacod <<
1302 { 1270 {
1303 std::ostringstream message; << 1271 G4cerr << "ERROR - G4VTwistSurface::Boundary::GetBoundaryParameters()"
1304 message << "Located in the corner area." << 1272 << G4endl
1305 << " This function returns << 1273 << " Located in the corner area. This function"
1306 << "a boundary line." << G4endl << 1274 << " returns a direction vector of a boundary line."
1307 << " areacode = " << areac << 1275 << " areacode = " << areacode << G4endl;
1308 G4Exception("G4VTwistSurface::Boundary::G 1276 G4Exception("G4VTwistSurface::Boundary::GetBoundaryParameters()",
1309 "GeomSolids0003", FatalExcept << 1277 "InvalidCondition", FatalException,
>> 1278 "Located in the corner area.");
1310 } 1279 }
1311 if ((areacode & sSizeMask) != (fBoundaryAco 1280 if ((areacode & sSizeMask) != (fBoundaryAcode & sSizeMask))
1312 { 1281 {
1313 return false; 1282 return false;
1314 } 1283 }
1315 d = fBoundaryDirection; 1284 d = fBoundaryDirection;
1316 x0 = fBoundaryX0; 1285 x0 = fBoundaryX0;
1317 boundarytype = fBoundaryType; 1286 boundarytype = fBoundaryType;
1318 return true; 1287 return true;
1319 } 1288 }
1320 1289